/******************************************************************************
* SBzr-Aux.c - Bezier 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 bezier 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 subdivising the surface one row/col at a *
* time, using the Bezier curve subdivision.				      *
******************************************************************************/
CagdSrfStruct *BzrSrfSubdivAtParam(CagdSrfStruct *Srf, CagdRType t,
							CagdSrfDirType Dir)
{
    int Row, Col,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength;
    CagdCrvStruct *Crv, *LCrv, *RCrv;
    CagdSrfStruct
	*RSrf = BzrSrfNew(ULength, VLength, Srf ->PType),
	*LSrf = BzrSrfNew(ULength, VLength, Srf ->PType);

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    for (Row = 0; Row < VLength; Row++) {
		Crv = BzrSrfCrvFromMesh(Srf, Row, CAGD_CONST_V_DIR);
		LCrv = BzrCrvSubdivAtParam(Crv, t);
		RCrv = LCrv -> Pnext;
		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 = BzrSrfCrvFromMesh(Srf, Col, CAGD_CONST_U_DIR);
		LCrv = BzrCrvSubdivAtParam(Crv, t);
		RCrv = LCrv -> Pnext;
		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.						      *
* Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:    *
*		       i	    k-i					      *
* Q(0) = P(0), Q(i) = --- P(i-1) + (---) P(i), Q(k) = P(k-1).		      *
*		       k	     k					      *
* This is applied to all rows/cols of the surface.			      *
******************************************************************************/
CagdSrfStruct *BzrSrfDegreeRaise(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, Row, Col,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct
	*RaisedSrf = NULL;

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    RaisedSrf = BzrSrfNew(ULength, VLength + 1, Srf -> PType);

	    for (Col = 0; Col < ULength; Col++) {
		for (j = IsNotRational; j <= MaxCoord; j++)	    /* Q(0). */
		    RaisedSrf -> Points[j][RAISED_SRF(Col, 0)] =
		       Srf -> Points[j][SRF(Col, 0)];

		for (i = 1; i < VLength; i++)			    /* Q(i). */
		    for (j = IsNotRational; j <= MaxCoord; j++)
			RaisedSrf -> Points[j][RAISED_SRF(Col, i)] =
			    Srf -> Points[j][SRF(Col, i - 1)] *
			    		         (i / ((CagdRType) VLength)) +
			    Srf -> Points[j][SRF(Col, i)] *
				     ((VLength - i) / ((CagdRType) VLength));

		for (j = IsNotRational; j <= MaxCoord; j++)	    /* Q(k). */
		    RaisedSrf -> Points[j][RAISED_SRF(Col, VLength)] =
			Srf -> Points[j][SRF(Col, VLength - 1)];
            }
	    break;
	case CAGD_CONST_V_DIR:
	    RaisedSrf = BzrSrfNew(ULength + 1, VLength, Srf -> PType);

	    for (Row = 0; Row < VLength; Row++) {
		for (j = IsNotRational; j <= MaxCoord; j++)	    /* Q(0). */
		    RaisedSrf -> Points[j][RAISED_SRF(0, Row)] =
		       Srf -> Points[j][SRF(0, Row)];

		for (i = 1; i < ULength; i++)			    /* Q(i). */
		    for (j = IsNotRational; j <= MaxCoord; j++)
			RaisedSrf -> Points[j][RAISED_SRF(i, Row)] =
			    Srf -> Points[j][SRF(i - 1, Row)] *
			    		         (i / ((CagdRType) ULength)) +
			    Srf -> Points[j][SRF(i, Row)] *
	    			     ((ULength - i) / ((CagdRType) ULength));

		for (j = IsNotRational; j <= MaxCoord; j++)	    /* Q(k). */
		    RaisedSrf -> Points[j][RAISED_SRF(ULength, Row)] =
			Srf -> Points[j][SRF(ULength - 1, Row)];
	    }
	    break;
	default:
	    FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return RaisedSrf;
}

/******************************************************************************
* 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 *BzrSrfDerive(CagdSrfStruct *Srf, CagdSrfDirType Dir)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_SRF(Srf);
    int i, j, Row, Col,
	ULength = Srf -> ULength,
	VLength = Srf -> VLength,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Srf -> PType);
    CagdSrfStruct
        *DerivedSrf = NULL;


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

	    DerivedSrf = BzrSrfNew(ULength, VLength - 1, Srf -> PType);

	    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;
	case CAGD_CONST_V_DIR:
	    if (!IsNotRational || ULength < 3)
		FATAL_ERROR(CAGD_ERR_RAT_LIN_NO_SUPPORT);

	    DerivedSrf = BzrSrfNew(ULength - 1, VLength, Srf -> PType);

	    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;
	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 *BzrSrfTangent(CagdSrfStruct *Srf, CagdRType u, CagdRType v,
							 CagdSrfDirType Dir)
{
    CagdVecStruct
	*Tangent = NULL;
    CagdCrvStruct *Crv;

    switch (Dir) {
	case CAGD_CONST_V_DIR:
	    Crv = BzrSrfCrvFromSrf(Srf, v, Dir);
	    Tangent = BzrCrvTangent(Crv, u);
	    CagdCrvFree(Crv);
	    break;
	case CAGD_CONST_U_DIR:
	    Crv = BzrSrfCrvFromSrf(Srf, u, Dir);
	    Tangent = BzrCrvTangent(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 *BzrSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
{
    static CagdVecStruct Normal;
    CagdVecStruct *V, T1, T2;

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

    V = BzrSrfTangent(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;
}

/******************************************************************************
* Convert a Bezier srf into Bspline srf by adding two open knot vectors.      *
******************************************************************************/
CagdSrfStruct *CnvrtBezier2BsplineSrf(CagdSrfStruct *Srf)
{
    CagdSrfStruct *BspSrf = CagdSrfCopy(Srf);

    BspSrf -> UOrder = BspSrf -> ULength;
    BspSrf -> VOrder = BspSrf -> VLength;
    BspSrf -> UKnotVector = BspKnotUniformOpen(BspSrf -> ULength,
						    BspSrf -> UOrder, NULL);
    BspSrf -> VKnotVector = BspKnotUniformOpen(BspSrf -> VLength,
						    BspSrf -> VOrder, NULL);
    BspSrf -> GType = CAGD_SBSPLINE_TYPE;
    return BspSrf;
}