/******************************************************************************
* SBsp-Aux.c - Bspline surface auxilary routines.			      *
*******************************************************************************
* Written by Gershon Elber, July. 90.					      *
******************************************************************************/

#ifdef __MSDOS__
#include <stdlib.h>
#endif /* __MSDOS__ */

#include <ctype.h>
#include <stdio.h>
#include <string.h>
#include "cagd_loc.h"

/* Define some marcos to make some of the routines below look better. They  */
/* calculate the index of the U, V point of the control mesh in Points.	    */
#define DERIVED_SRF(U, V)	CAGD_MESH_UV(DerivedSrf, U, V)
#define RAISED_SRF(U, V)	CAGD_MESH_UV(RaisedSrf, U, V)
#define SRF(U, V)		CAGD_MESH_UV(Srf, U, V)

/******************************************************************************
* Given a bspline surface - subdivide it into two at given parametric value.  *
* Returns pointer to first surface in a list of two srfs (subdivided ones).   *
* The subdivision is exact result of evaluating the surface int a curve at t  *
* using the recursive algorithm - the left resulting points is left surface,  *
* and the right resulting points is right surface (left is below t).	      *
******************************************************************************/
CagdSrfStruct *BspSrfSubdivAtParam(CagdSrfStruct *Srf, CagdRType t,
							CagdSrfDirType Dir)
{
    int Row, Col,
        ULength1 = 0,
	ULength2 = 0,
	VLength1 = 0,
	VLength2 = 0,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder;
    CagdCrvStruct *Crv, *LCrv, *RCrv;
    CagdSrfStruct
	*RSrf = NULL,
	*LSrf = NULL;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    for (Row = 0; Row < VLength; Row++) {
		Crv = BspSrfCrvFromMesh(Srf, Row, CAGD_CONST_V_DIR);
		LCrv = BspCrvSubdivAtParam(Crv, t);
		RCrv = LCrv -> Pnext;

		if (Row == 0) {    /* Figure out surfaces size and allocate. */
		    ULength1 = LCrv -> Length;
		    ULength2 = RCrv -> Length;
		    LSrf = BspSrfNew(ULength1, VLength, UOrder, VOrder,
								Srf -> PType);
		    RSrf = BspSrfNew(ULength2, VLength, UOrder, VOrder,
								Srf -> PType);
		    CagdFree((VoidPtr) RSrf -> UKnotVector);
		    CagdFree((VoidPtr) RSrf -> VKnotVector);
		    CagdFree((VoidPtr) LSrf -> UKnotVector);
		    CagdFree((VoidPtr) LSrf -> VKnotVector);
		    RSrf -> UKnotVector = BspKnotCopy(RCrv -> KnotVector,
					      RCrv -> Length + RCrv -> Order);
		    RSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
					      Srf -> VLength + Srf -> VOrder);
		    LSrf -> UKnotVector = BspKnotCopy(LCrv -> KnotVector,
					      LCrv -> Length + LCrv -> Order);
		    LSrf -> VKnotVector = BspKnotCopy(Srf -> VKnotVector,
					      Srf -> VLength + Srf -> VOrder);
		}

		CagdCrvToMesh(LCrv, Row, CAGD_CONST_V_DIR, LSrf);
		CagdCrvToMesh(RCrv, Row, CAGD_CONST_V_DIR, RSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(LCrv);
		CagdCrvFree(RCrv);
	    }
	    break;
	case CAGD_CONST_V_DIR:
	    for (Col = 0; Col < ULength; Col++) {
		Crv = BspSrfCrvFromMesh(Srf, Col, CAGD_CONST_U_DIR);
		LCrv = BspCrvSubdivAtParam(Crv, t);
		RCrv = LCrv -> Pnext;

		if (Col == 0) {    /* Figure out surfaces size and allocate. */
		    VLength1 = LCrv -> Length;
		    VLength2 = RCrv -> Length;
		    LSrf = BspSrfNew(ULength, VLength1, UOrder, VOrder,
								Srf -> PType);
		    RSrf = BspSrfNew(ULength, VLength2, UOrder, VOrder,
								Srf -> PType);
		    CagdFree((VoidPtr) RSrf -> UKnotVector);
		    CagdFree((VoidPtr) RSrf -> VKnotVector);
		    CagdFree((VoidPtr) LSrf -> UKnotVector);
		    CagdFree((VoidPtr) LSrf -> VKnotVector);
		    RSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
					      Srf -> ULength + Srf -> UOrder);
		    RSrf -> VKnotVector = BspKnotCopy(RCrv -> KnotVector,
					      RCrv -> Length + RCrv -> Order);
		    LSrf -> UKnotVector = BspKnotCopy(Srf -> UKnotVector,
					      Srf -> ULength + Srf -> UOrder);
		    LSrf -> VKnotVector = BspKnotCopy(LCrv -> KnotVector,
					      LCrv -> Length + LCrv -> Order);
		}

		CagdCrvToMesh(LCrv, Col, CAGD_CONST_U_DIR, LSrf);
		CagdCrvToMesh(RCrv, Col, CAGD_CONST_U_DIR, RSrf);

		CagdCrvFree(Crv);
		CagdCrvFree(LCrv);
		CagdCrvFree(RCrv);
	    }
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    LSrf -> Pnext = RSrf;
    return LSrf;
}

/******************************************************************************
* Return a new surface, identical to the original but with one degree higher  *
* in the given direction.						      *
******************************************************************************/
CagdSrfStruct *BspSrfDegreeRaise(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    FATAL_ERROR(CAGD_ERR_NOT_IMPLEMENTED);

    return NULL;
}

/******************************************************************************
* Return a new surface equal to the derived surface in the direction Dir.     *
* A test is made to make sure the original surface is not rational.	      *
* Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:    *
* Q(i) = (k - 1) * P(i+1) - P(i), i = 0 to k-2.				      *
* This is applied to all rows/cols of the surface.			      *
******************************************************************************/
CagdSrfStruct *BspSrfDerive(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, Row, Col,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength,
	UOrder = Srf -> UOrder,
	VOrder = Srf -> VOrder,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct
        *DerivedSrf = NULL;


    switch (Dir) {
	case CAGD_CONST_V_DIR:
	    if (!IsNotRational || UOrder < 3)
		FATAL_ERROR(CAGD_ERR_RAT_LIN_NO_SUPPORT);

	    DerivedSrf = BspSrfNew(ULength - 1, VLength,
				   UOrder - 1, VOrder, Srf -> PType);
	    GEN_COPY(DerivedSrf -> UKnotVector,
		     &Srf -> UKnotVector[1],
		     sizeof(CagdRType) * (ULength + UOrder - 2));
	    GEN_COPY(DerivedSrf -> VKnotVector,
		     Srf -> VKnotVector,
		     sizeof(CagdRType) * (VLength + VOrder));

       	    for (Row = 0; Row < VLength; Row++)
		for (i = 0; i < ULength - 1; i++)
		    for (j = IsNotRational; j <= MaxCoord; j++)
			DerivedSrf -> Points[j][DERIVED_SRF(i, Row)] =
                            (ULength - 1) *
			    (Srf -> Points[j][SRF(i + 1, Row)] -
			     Srf -> Points[j][SRF(i, Row)]);
	    break;
	case CAGD_CONST_U_DIR:
	    if (!IsNotRational || VOrder < 3)
		FATAL_ERROR(CAGD_ERR_RAT_LIN_NO_SUPPORT);

	    DerivedSrf = BspSrfNew(ULength, VLength - 1,
	    			   UOrder, VOrder - 1, Srf -> PType);
	    GEN_COPY(DerivedSrf -> UKnotVector,
		     Srf -> UKnotVector,
		     sizeof(CagdRType) * (ULength + UOrder));
	    GEN_COPY(DerivedSrf -> VKnotVector,
		     &Srf -> VKnotVector[1],
		     sizeof(CagdRType) * (VLength + VOrder - 2));

	    for (Col = 0; Col < ULength; Col++)
		for (i = 0; i < VLength - 1; i++)
		    for (j = IsNotRational; j <= MaxCoord; j++)
			DerivedSrf -> Points[j][DERIVED_SRF(Col, i)] =
                            (VLength - 1) *
			    (Srf -> Points[j][SRF(Col, i + 1)] -
			     Srf -> Points[j][SRF(Col, i)]);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return DerivedSrf;
}

/******************************************************************************
* Evaluate the tangent to a surface at a given point and given direction.     *
******************************************************************************/
CagdVecStruct *BspSrfTangent(CagdSrfStruct *Srf, CagdRType u, CagdRType v,
							 CagdSrfDirType Dir)
{
    CagdVecStruct
	*Tangent = NULL;
    CagdCrvStruct *Crv;

    switch (Dir) {
	case CAGD_CONST_V_DIR:
	    Crv = BspSrfCrvFromSrf(Srf, v, Dir);
	    Tangent = BspCrvTangent(Crv, u);
	    CagdCrvFree(Crv);
	    break;
	case CAGD_CONST_U_DIR:
	    Crv = BspSrfCrvFromSrf(Srf, u, Dir);
	    Tangent = BspCrvTangent(Crv, v);
	    CagdCrvFree(Crv);
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return Tangent;
}

/******************************************************************************
* Evaluate the normal of a surface at a given point.			      *
******************************************************************************/
CagdVecStruct *BspSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
{
    static CagdVecStruct Normal;
    CagdVecStruct *V, T1, T2;

    V = BspSrfTangent(Srf, u, v, CAGD_CONST_U_DIR);
    CAGD_COPY_VECTOR(T1, *V);

    V = BspSrfTangent(Srf, u, v, CAGD_CONST_V_DIR);
    CAGD_COPY_VECTOR(T2, *V);

    /* The normal is the cross product of T1 and T2: */
    Normal.Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1];
    Normal.Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2];
    Normal.Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0];

    CAGD_NORMALIZE_VECTOR(Normal);		    /* Normalize the vector. */

    return &Normal;
}

/******************************************************************************
* Evaluate the normals of a surface at a mesh defined by subdividing the      *
* parametric space int a grid of size UFineNess by VFineNess.		      *
* The normals are saved in a linear CagdVecStruct vector which is allocated   *
* dynamically. Data is saved u inc. first.				      *
* This routine much faster than evaluating normal per each point of such mesh.*
******************************************************************************/
CagdVecStruct *BspSrfMeshNormals(CagdSrfStruct *Srf, int UFineNess,
								int VFineNess)
{
    int i, j;
    CagdRType u, v, UMin, UMax, VMin, VMax;
    CagdVecStruct *Normals, *NPtr, *T, T1, T2;
    CagdCrvStruct **UCurves, **VCurves, *UCrv, *VCrv;

    BspSrfDomain(Srf, &UMin, &UMax, &VMin, &VMax);
    Normals = (CagdVecStruct *) CagdMalloc(sizeof(CagdVecStruct) * UFineNess *
								   VFineNess);

    UCurves = (CagdCrvStruct **) CagdMalloc(sizeof(CagdCrvStruct *) *
								   UFineNess);
    VCurves = (CagdCrvStruct **) CagdMalloc(sizeof(CagdCrvStruct *) *
								   VFineNess);

    UFineNess--;
    VFineNess--;
    for (i = 0; i <= UFineNess; i++)		    /* Prepare Iso U curves. */
	UCurves[i] = BspSrfCrvFromSrf(Srf,
    			UMin + (UMax - UMin) * ((CagdRType) i) / UFineNess,
			CAGD_CONST_U_DIR);
    for (j = 0; j <= VFineNess; j++)	            /* Prepare Iso V curves. */
	VCurves[j] = BspSrfCrvFromSrf(Srf,
			VMin + (VMax - VMin) * ((CagdRType) j) / VFineNess,
			CAGD_CONST_V_DIR);

    NPtr = Normals;
    for (i = 0; i <= UFineNess; i++) {
	UCrv = UCurves[i];
	u = UMin + (UMax - UMin) * ((CagdRType) i) / UFineNess;

	for (j = 0; j <= VFineNess; j++) {
	    VCrv = VCurves[j];
	    v = VMin + (VMax - VMin) * ((CagdRType) j) / VFineNess;

	    /* We need to copy the tangents as BspCrvTangent save it in as   */
	    /* static so the second call will overwrite first call value.    */
	    T = BspCrvTangent(UCrv, v);
	    CAGD_COPY_VECTOR(T1, *T);
	    T = BspCrvTangent(VCrv, u);
	    CAGD_COPY_VECTOR(T2, *T);

	    /* The normal is the cross product of T1 and T2: */
	    NPtr -> Vec[0] = T1.Vec[1] * T2.Vec[2] - T1.Vec[2] * T2.Vec[1];
	    NPtr -> Vec[1] = T1.Vec[2] * T2.Vec[0] - T1.Vec[0] * T2.Vec[2];
	    NPtr -> Vec[2] = T1.Vec[0] * T2.Vec[1] - T1.Vec[1] * T2.Vec[0];

	    CAGD_NORMALIZE_VECTOR(*NPtr);           /* Normalize the vector. */
	    NPtr++;
	}
    }

    for (i = 0; i <= UFineNess; i++) CagdCrvFree(UCurves[i]);
    CagdFree(UCurves);
    for (j = 0; j <= VFineNess; j++) CagdCrvFree(VCurves[j]);
    CagdFree(VCurves);

    return Normals;
}