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

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a Bezier surface - subdivides it into two sub-surfaces at the given  M
* parametric value.                                                          M
*   Returns pointer to first surface in a list of two subdivided surfaces.   M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:      To subdivide at parameter value t.                             M
*   t:        Parameter value to subdivide Srf at.                           M
*   Dir:      Direction of subdivision. Either U or V.                       M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A list of the two subdivided surfaces.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrSrfSubdivAtParam, subdivision, refinement                             M
*****************************************************************************/
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:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    LSrf -> Pnext = RSrf;
    return LSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new Bezier surface, identical to the original but with one       M
* degree higher, in the requested direction Dir.                             M
* Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then:   M
*		       i	    k-i					     V
* Q(0) = P(0), Q(i) = --- P(i-1) + (---) P(i), Q(k) = P(k-1).		     V
*		       k	     k					     V
* This is applied to all rows/cols of the surface.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To raise it degree by one.                                   M
*   Dir:        Direction to degree raise. Either U or V.                    M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A surface with one degree higher in direction Dir,     M
*                     representing the same geometry as Srf.		     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrSrfDegreeRaise, degree raising                                        M
*****************************************************************************/
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:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return RaisedSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new surface equal to the given surface, differentiated once in   M
* the direction Dir.							     M
*   Let old control polygon be P(i), i = 0 to k-1, and Q(i) be new one then: M
* Q(i) = (k - 1) * (P(i+1) - P(i)), i = 0 to k-2.			     V
* This is applied to all rows/cols of the surface.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        To differentiate.                                            M
*   Dir:        Direction of differentiation. Either U or V.                 M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:   Differentiated surface.                               M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrSrfDerive, derivatives, partial derivatives                           M
*****************************************************************************/
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:
	    DerivedSrf = BzrSrfNew(MAX(ULength - 1, 1), VLength, Srf -> PType);

	    for (Row = 0; Row < VLength; Row++)
		for (i = 0; i < MAX(ULength - 1, 1); i++)
		    for (j = IsNotRational; j <= MaxCoord; j++)
			DerivedSrf -> Points[j][DERIVED_SRF(i, Row)] =
			    ULength < 2 ? 0.0
					: (ULength - 1) *
					   (Srf -> Points[j][SRF(i + 1, Row)] -
					    Srf -> Points[j][SRF(i, Row)]);
	    break;
	case CAGD_CONST_V_DIR:
	    DerivedSrf = BzrSrfNew(ULength, MAX(VLength - 1, 1), Srf -> PType);

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

    return DerivedSrf;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluates the unit tangent to a surface at a given parametric location     M
( u, v) and given direction Dir.                                             M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Bezier surface to evaluate tangent vector for.               M
*   u, v:       Parametric location of required unit tangent.                M
*   Dir:        Direction of tangent vector. Either U or V.                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:  A pointer to a static vector holding the unit tangent  M
*                     information.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrSrfTangent, tangent                                                   M
*****************************************************************************/
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:
	    CAGD_FATAL_ERROR(CAGD_ERR_DIR_NOT_CONST_UV);
	    break;
    }

    return Tangent;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Evaluate the unit normal of a surface at a given parametric location.      M
*   If we fail to compute the normal at given location we retry by moving a  M
* tad.                                                                       M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:        Bezier surface to evaluate normal vector for.                M
*   u, v:       Parametric location of required unit normal.                 M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:  A pointer to a static vector holding the unit normal   M
*                     information.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrSrfNormal, normal                                                     M
*****************************************************************************/
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;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Converts a bezier surface into a Bspline surface by adding open end        M
* knot vector with no interior knots.                                        M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       Bezier surface to convert to a Bspline surface.               M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A Bspline surface representing same geometry as Srf.   M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtBezier2BsplineSrf, conversion                                       M
*****************************************************************************/
CagdSrfStruct *CnvrtBezier2BsplineSrf(CagdSrfStruct *Srf)
{
    CagdSrfStruct *BspSrf;

    if (Srf -> GType != CAGD_SBEZIER_TYPE) {
	CAGD_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;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Convert a Bspline surface into a set of Bezier surfaces by subdiving the   M
* Bspline surface at all its internal knots.				     M
*   Returned is a list of Bezier surface.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Srf:       Bspline surface to convert to a Bezier surface.               M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdSrfStruct *:  A list of Bezier surfaces representing same geometry   M
*                     as Srf.						     M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtBezier2BsplineSrf, conversion                                       M
*****************************************************************************/
CagdSrfStruct *CnvrtBspline2BezierSrf(CagdSrfStruct *Srf)
{
    CagdBType
	NewSrf = FALSE;
    int i, UOrder, ULength;
    CagdRType LastT, *UKnotVector;
    CagdSrfStruct *BezierSrfsAux, *TSrf, *OrigSrf,
	*BezierSrfs = NULL;

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

    if (!BspSrfHasOpenEC(Srf)) {
	Srf = BspSrfOpenEnd(Srf);
	NewSrf = TRUE;
    }
    UOrder = Srf -> UOrder;
    ULength = Srf -> ULength;
    UKnotVector = Srf -> UKnotVector;
    OrigSrf = Srf;

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

    if (NewSrf)
	CagdSrfFree(Srf);

    return BezierSrfs;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Auxiliary function of CnvrtBspline2BezierSrf. Does the other V direction's *
* subdivision.	                                                             *
*                                                                            *
* PARAMETERS:                                                                *
*   Srf:       To subdivide at all its interior knots in the V direction.    *
*                                                                            *
* RETURN VALUE:                                                              *
*   CagdSrfStruct *:  A list of Bezier surfaces representing same geometry   *
*                     as Srf.						     *
*****************************************************************************/
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;
}


