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

#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)

static CagdSrfStruct *CnvrtBspline2BezierSrfAux(CagdSrfStruct *Srf);

/******************************************************************************
* 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.     *
* 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;

    if (!IsNotRational)
	return BzrSrfDeriveRational(Srf, Dir);

    switch (Dir) {
	case CAGD_CONST_U_DIR:
	    if (VLength < 2)
		FATAL_ERROR(CAGD_ERR_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 (ULength < 2)
		FATAL_ERROR(CAGD_ERR_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.			      *
* If we fail to compute the normal at given location we try by moving a tad.  *
******************************************************************************/
CagdVecStruct *BzrSrfNormal(CagdSrfStruct *Srf, CagdRType u, CagdRType v)
{
    static CagdVecStruct Normal;
    CagdVecStruct *V, T1, T2;

    V = BzrSrfTangent(Srf, u, v, CAGD_CONST_U_DIR);
    if (CAGD_LEN_VECTOR(*V) < EPSILON)
	V = BzrSrfTangent(Srf,
			  u > 0.5 ? u - EPSILON : u + EPSILON,
			  v > 0.5 ? v - EPSILON : v + EPSILON,
			  CAGD_CONST_U_DIR);
    CAGD_COPY_VECTOR(T1, *V);

    V = BzrSrfTangent(Srf, u, v, CAGD_CONST_V_DIR);
    if (CAGD_LEN_VECTOR(*V) < EPSILON)
	V = BzrSrfTangent(Srf,
			  u > 0.5 ? u - EPSILON : u + EPSILON,
			  v > 0.5 ? v - EPSILON : v + EPSILON,
			  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;

    if (Srf -> GType != CAGD_SBEZIER_TYPE) {
	FATAL_ERROR(CAGD_ERR_WRONG_SRF);
	return NULL;
    }

    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;
}

/******************************************************************************
* Convert a Bspline srf into a set of Bezier srfs by subdiving the Bspline    *
* surface at all its internal knots.					      *
*   Returned is a list of Bezier surface.				      *
******************************************************************************/
CagdSrfStruct *CnvrtBspline2BezierSrf(CagdSrfStruct *Srf)
{
    int i,
	UOrder = Srf -> UOrder,
	ULength = Srf -> ULength;
    CagdRType LastT,
	*UKnotVector = Srf -> UKnotVector;
    CagdSrfStruct *BezierSrfsAux, *TSrf,
	*BezierSrfs = NULL,
	*OrigSrf = Srf;

    if (Srf -> GType != CAGD_SBSPLINE_TYPE) {
	FATAL_ERROR(CAGD_ERR_WRONG_SRF);
	return NULL;
    }

    for (i = ULength - 1, LastT = UKnotVector[ULength]; i >= UOrder; i--) {
    	CagdRType
    	    t = UKnotVector[i];
    	    
	if (!APX_EQ(LastT, t)) {
    	    CagdSrfStruct
    		*Srfs = BspSrfSubdivAtParam(Srf, t, CAGD_CONST_U_DIR);

    	    if (Srf != OrigSrf)
    	        CagdSrfFree(Srf);

	    BezierSrfsAux = CnvrtBspline2BezierSrfAux(Srfs -> Pnext);
	    for (TSrf = BezierSrfsAux; TSrf -> Pnext != NULL; TSrf = TSrf -> Pnext);
	    TSrf -> Pnext = BezierSrfs;
	    BezierSrfs = BezierSrfsAux;
	    CagdSrfFree(Srfs -> Pnext);

    	    Srf = Srfs;
    	    Srf -> Pnext = NULL;
    	}
    }

    if (Srf == OrigSrf)
	BezierSrfs = CnvrtBspline2BezierSrfAux(Srf);
    else {
	BezierSrfsAux = CnvrtBspline2BezierSrfAux(Srf);
	for (TSrf = BezierSrfsAux; TSrf -> Pnext != NULL; TSrf = TSrf -> Pnext);
	TSrf -> Pnext = BezierSrfs;
	BezierSrfs = BezierSrfsAux;

	CagdSrfFree(Srf);
    }

    return BezierSrfs;
}


/******************************************************************************
* Aux. function to CnvrtBspline2BezierSrf. Does the other subdivision.	      *
******************************************************************************/
static CagdSrfStruct *CnvrtBspline2BezierSrfAux(CagdSrfStruct *Srf)
{
    int i,
	VOrder = Srf -> VOrder,
	VLength = Srf -> VLength;
    CagdRType LastT,
	*VKnotVector = Srf -> VKnotVector;
    CagdSrfStruct
	*BezierSrfs = NULL,
	*OrigSrf = Srf;

    for (i = VLength - 1, LastT = VKnotVector[VLength]; i >= VOrder; i--) {
    	CagdRType
    	    t = VKnotVector[i];
    	    
	if (!APX_EQ(LastT, t)) {
    	    CagdSrfStruct
    		*Srfs = BspSrfSubdivAtParam(Srf, t, CAGD_CONST_V_DIR);

    	    if (Srf != OrigSrf)
    	        CagdSrfFree(Srf);

    	    Srfs -> Pnext -> Pnext = BezierSrfs;
	    BezierSrfs = Srfs -> Pnext;

    	    Srf = Srfs;
    	    Srf -> Pnext = NULL;
    	}
    }

    if (Srf == OrigSrf) {
	/* No interior knots in this surface - just copy it: */
        BezierSrfs = CagdSrfCopy(Srf);
    }
    else {
    	Srf -> Pnext = BezierSrfs;
    	BezierSrfs = Srf;
    }

    for (Srf = BezierSrfs; Srf != NULL; Srf = Srf -> Pnext) {
    	Srf -> GType = CAGD_SBEZIER_TYPE;
	IritFree((VoidPtr) Srf -> UKnotVector);
	IritFree((VoidPtr) Srf -> VKnotVector);
	Srf -> UKnotVector = NULL;
	Srf -> VKnotVector = NULL;
    }
    
    return BezierSrfs;
}