/******************************************************************************
* CBzr-Aux.c - Bezier curve auxilary routines.				      *
*******************************************************************************
* Written by Gershon Elber, Mar. 90.					      *
******************************************************************************/

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

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a Bezier curve - subdivides it into two sub-curves at the given      M
* parametric value.                                                          M
*   Returns pointer to first curve in a list of two subdivided curves.       M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       To subdivide at parametr value t.                             M
*   t:         Parameter value to subdivide Crv at.                          M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A list of the two subdivided curves.                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvSubdivAtParam, subdivision, refinement                             M
*****************************************************************************/
CagdCrvStruct *BzrCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, l,
	k = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType
	t1 = 1.0 - t;
    CagdCrvStruct *LCrv, *RCrv;

    if (t < 0.0 || t > 1.0)
	CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);

    LCrv = BzrCrvNew(k, Crv -> PType);
    RCrv = BzrCrvNew(k, Crv -> PType);

    /* Copy Curve into RCrv, so we can apply the recursive algo. to it.      */
    for (i = 0; i < k; i++)
	for (j = IsNotRational; j <= MaxCoord; j++)
	    RCrv -> Points[j][i] = Crv -> Points[j][i];

    for (j = IsNotRational; j <= MaxCoord; j++)
	LCrv -> Points[j][0] = Crv -> Points[j][0];

    /* Apply the recursive algorithm to RCrv, and update LCrv with the	     */
    /* temporary results. Note we updated the first point of LCrv above.     */
    for (i = 1; i < k; i++) {
	for (l = 0; l < k - i; l++)
	    for (j = IsNotRational; j <= MaxCoord; j++)
		RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 +
					 RCrv -> Points[j][l+1] * t;
	/* Copy temporary result to LCrv: */
	for (j = IsNotRational; j <= MaxCoord; j++)
	    LCrv -> Points[j][i] = RCrv -> Points[j][0];
    }
    LCrv -> Pnext = RCrv;

    return LCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new curve, identical to the original but with order N.	     M
*   Degree raise is computed by multiplying by a constant 1 curve of order   M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To raise its degree to a NewOrder.                           M
*   NewOrder:   NewOrder for Crv.                                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A curve of order NewOrder representing the same        M
*                     geometry as Crv.                                       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvDegreeRaiseN, degree raising                                       M
*****************************************************************************/
CagdCrvStruct *BzrCrvDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder)
{
    int i, j, RaisedOrder,
	Order = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *RaisedCrv, *UnitCrv;

    if (NewOrder < Order) {
	CAGD_FATAL_ERROR(CAGD_ERR_WRONG_ORDER);
	return NULL;
    }
    RaisedOrder = NewOrder - Order + 1;

    UnitCrv = BzrCrvNew(RaisedOrder, CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
    for (i = 1; i <= MaxCoord; i++)
	for (j = 0; j < RaisedOrder; j++)
	    UnitCrv -> Points[i][j] = 1.0;

    RaisedCrv = BzrCrvMult(Crv, UnitCrv);

    CagdCrvFree(UnitCrv);

    return RaisedCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new curve, identical to the original but with one degree higher. 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
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To raise it degree by one.                                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A curve of one order higher representing the same      M
*                     geometry as Crv.                                       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvDegreeRaise, degree raising                                        M
*****************************************************************************/
CagdCrvStruct *BzrCrvDegreeRaise(CagdCrvStruct *Crv)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j,
	k = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct
	*RaisedCrv = BzrCrvNew(k + 1, Crv -> PType);

    for (j = IsNotRational; j <= MaxCoord; j++)			    /* Q(0). */
	RaisedCrv -> Points[j][0] = Crv -> Points[j][0];

    for (i = 1; i < k; i++)					    /* Q(i). */
	for (j = IsNotRational; j <= MaxCoord; j++)
	    RaisedCrv -> Points[j][i] =
		Crv -> Points[j][i-1] * (i / ((CagdRType) k)) +
		Crv -> Points[j][i] * ((k - i) / ((CagdRType) k));

    for (j = IsNotRational; j <= MaxCoord; j++)			    /* Q(k). */
	RaisedCrv -> Points[j][k] = Crv -> Points[j][k-1];

    return RaisedCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a unit vector, equal to the tangent to Crv at parameter value t.   M
*   Algorithm: pseudo subdivide Crv at t and using control point of          M
* subdivided curve find the tangent as the difference of the 2 end points.   M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       Crv for which to compute a unit tangent.                      M
*   t:         The parameter at which to compute the unit tangent.           M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:  A pointer to a static vector holding the tangent       M
*                     information.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvTangent, tangent                                                   M
*****************************************************************************/
CagdVecStruct *BzrCrvTangent(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdVecStruct P2;
    CagdVecStruct P1, *T;
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, l,
	k = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType
	t1 = 1.0 - t;
    CagdCrvStruct *RCrv;

    if (APX_EQ(t, 0.0)) {
	/* Use Crv starting tangent direction. */
	CagdCoerceToE3(P1.Vec, Crv -> Points, 0, Crv -> PType);
	CagdCoerceToE3(P2.Vec, Crv -> Points, 1, Crv -> PType);
    }
    else if (APX_EQ(t, 1.0)) {
	/* Use Crv ending tangent direction. */
	CagdCoerceToE3(P1.Vec, Crv -> Points, k - 2, Crv -> PType);
	CagdCoerceToE3(P2.Vec, Crv -> Points, k - 1, Crv -> PType);
    }
    else {
	if (t < 0.0 || t > 1.0)
	    CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);

	RCrv = BzrCrvNew(k, Crv -> PType);

	/* Copy Crv into RCrv, so we can apply the recursive algo. to it. */
	for (i = 0; i < k; i++)
	    for (j = IsNotRational; j <= MaxCoord; j++)
		RCrv -> Points[j][i] = Crv -> Points[j][i];

	/* Apply the recursive algorithm to RCrv. */
	for (i = 1; i < k; i++)
	    for (l = 0; l < k - i; l++)
		for (j = IsNotRational; j <= MaxCoord; j++)
		    RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 +
					   RCrv -> Points[j][l+1] * t;

	CagdCoerceToE3(P1.Vec, RCrv -> Points, 0, RCrv -> PType);
	CagdCoerceToE3(P2.Vec, RCrv -> Points, 1, RCrv -> PType);

	CagdCrvFree(RCrv);
    }

    CAGD_SUB_VECTOR(P2, P1);

    if (CAGD_LEN_VECTOR(P2) < IRIT_EPSILON) {
	if (AttrGetIntAttrib(Crv -> Attr, "_tan") != TRUE) {
	    /* Try to move a little. This location has zero speed. However,  */
	    /* do it only once since we can be here forever. The "_tan"      */
	    /* attribute guarantee we will try to move EPSILON only once.    */
	    AttrSetIntAttrib(&Crv -> Attr, "_tan", TRUE);

	    T = BzrCrvTangent(Crv, t < 0.5 ? t + EPSILON : t - EPSILON);

	    AttrFreeOneAttribute(&Crv -> Attr, "_tan");

	    return T;
	}
	else {
	    /* A zero length vector signals failure to compute tangent. */
	    return &P2;
	}
    }
    else {
	CAGD_NORMALIZE_VECTOR(P2);		    /* Normalize the vector. */

	return &P2;
    }
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a unit vector, equal to the binormal to Crv at parameter value t.  M
*   Algorithm: insert (order - 1) knots and using 3 consecutive control      M
* points at the refined location (p1, p2, p3), compute to binormal to be the M
* cross product of the two vectors (p1 - p2) and (p2 - p3).		     M
*   Since a curve may have not BiNormal at inflection points or if the 3     M
* points are colinear, NULL will be returned at such cases.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       Crv for which to compute a unit binormal.                     M
*   t:         The parameter at which to compute the unit binormal.          M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:  A pointer to a static vector holding the binormal      M
*                     information.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvBiNormal, binormal                                                 M
*****************************************************************************/
CagdVecStruct *BzrCrvBiNormal(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdVecStruct P3;
    CagdVecStruct P1, P2;
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, l,
	k = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType
	t1 = 1.0 - t;
    CagdCrvStruct *RCrv;

    /* Can not compute for linear curves. */
    if (k <= 2)
	return NULL;

    /* For planar curves, B is trivially the Z axis. */
    if (CAGD_NUM_OF_PT_COORD(Crv -> PType) == 2) {
	P3.Vec[0] = P3.Vec[1] = 0.0;
	P3.Vec[2] = 1.0;
	return &P3;
    }

    if (APX_EQ(t, 0.0)) {
	/* Use Crv starting tangent direction. */
	CagdCoerceToE3(P1.Vec, Crv -> Points, 0, Crv -> PType);
	CagdCoerceToE3(P2.Vec, Crv -> Points, 1, Crv -> PType);
	CagdCoerceToE3(P3.Vec, Crv -> Points, 2, Crv -> PType);
    }
    else if (APX_EQ(t, 1.0)) {
	/* Use Crv ending tangent direction. */
	CagdCoerceToE3(P1.Vec, Crv -> Points, k - 3, Crv -> PType);
	CagdCoerceToE3(P2.Vec, Crv -> Points, k - 2, Crv -> PType);
	CagdCoerceToE3(P3.Vec, Crv -> Points, k - 1, Crv -> PType);
    }
    else {
	if (t < 0.0 || t > 1.0)
	    CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);

	RCrv = BzrCrvNew(k, Crv -> PType);

	/* Copy Crv into RCrv, so we can apply the recursive algo. to it. */
	for (i = 0; i < k; i++)
	    for (j = IsNotRational; j <= MaxCoord; j++)
		RCrv -> Points[j][i] = Crv -> Points[j][i];

	/* Apply the recursive algorithm to RCrv. */
	for (i = 1; i < k; i++)
	    for (l = 0; l < k - i; l++)
		for (j = IsNotRational; j <= MaxCoord; j++)
		    RCrv -> Points[j][l] = RCrv -> Points[j][l] * t1 +
					   RCrv -> Points[j][l+1] * t;

	CagdCoerceToE3(P1.Vec, RCrv -> Points, 0, RCrv -> PType);
	CagdCoerceToE3(P2.Vec, RCrv -> Points, 1, RCrv -> PType);
	CagdCoerceToE3(P3.Vec, RCrv -> Points, 2, RCrv -> PType);

	CagdCrvFree(RCrv);
    }

    CAGD_SUB_VECTOR(P1, P2);
    CAGD_SUB_VECTOR(P2, P3);

    CROSS_PROD(P3.Vec, P1.Vec, P2.Vec);

    if ((t = CAGD_LEN_VECTOR(P3)) < IRIT_EPSILON)
	return NULL;
    else
	CAGD_DIV_VECTOR(P3, t);			    /* Normalize the vector. */

    return &P3;
}


/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a unit vector, equal to the normal of Crv at parameter value t.    M
*   Algorithm: returns the cross product of the curve tangent and binormal.  M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       Crv for which to compute a unit normal.                       M
*   t:         The parameter at which to compute the unit normal.            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdVecStruct *:  A pointer to a static vector holding the normal        M
*                     information.                                           M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvNoraml, normal                                                     M
*****************************************************************************/
CagdVecStruct *BzrCrvNormal(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdVecStruct N, *T, *B;

    T = BzrCrvTangent(Crv, t);
    B = BzrCrvBiNormal(Crv, t);

    if (T == NULL || B == NULL)
	return NULL;

    CROSS_PROD(N.Vec, T -> Vec, B -> Vec);

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

    return &N;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new curve, equal to the given curve, differentiated once.        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.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To differentiate.                                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Differentiated curve.                                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvDerive, derivatives                                                M
*****************************************************************************/
CagdCrvStruct *BzrCrvDerive(CagdCrvStruct *Crv)
{
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j,
	k = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *DerivedCrv;

    if (!IsNotRational)
	return BzrCrvDeriveRational(Crv);

    DerivedCrv = BzrCrvNew(MAX(1, k - 1), Crv -> PType);

    if (k >= 2) {
	for (i = 0; i < k - 1; i++)
	    for (j = IsNotRational; j <= MaxCoord; j++)
		DerivedCrv -> Points[j][i] =
		    (k - 1) * (Crv -> Points[j][i+1] - Crv -> Points[j][i]);
    }
    else {
	for (j = IsNotRational; j <= MaxCoord; j++)
	    DerivedCrv -> Points[j][0] = 0.0;
    }

    return DerivedCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new Bezier curve, equal to the integral of the given Bezier      M
* crv.		                                                             M
* The given Bezier curve should be nonrational.				     M
*									     V
* 	      n		   n		   n	   n+1			     V
*   /	     /-		   -	  /	   -   P    -			     V
*  |	    | \	    n	   \	 |  n	   \    i   \  n+1		     V
*  | C(t) = | / P  B (t) = / P   | B (t) = / -----  / B   (t) =		     V
* / 	   /  -	 i  i	   -  i /   i	   - n + 1  -  j		     V
*            i=0          i=0             i=0     j=i+1			     V
*									     V
*        n+1 j-1							     V
*         -   -								     V
*     1   \   \	    n+1							     V
* = ----- /   / P  B   (t)						     V
*   n + 1 -   -  i  j							     V
*        j=1 i=0							     V
*									     V
*									     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:         Curve to integrate.                                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Integrated curve.                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   BzrCrvIntegrate, integrals                                               M
*****************************************************************************/
CagdCrvStruct *BzrCrvIntegrate(CagdCrvStruct *Crv)
{
    int i, j, k,
	n = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *IntCrv;

    if (CAGD_IS_RATIONAL_CRV(Crv))
	CAGD_FATAL_ERROR(CAGD_ERR_RATIONAL_NO_SUPPORT);

    IntCrv = BzrCrvNew(n + 1, Crv -> PType);

    for (k = 1; k <= MaxCoord; k++) {
	CagdRType
	    *Points = Crv -> Points[k],
	    *IntPoints = IntCrv -> Points[k];

	for (j = 0; j < n + 1; j++) {
	    IntPoints[j] = 0.0;
	    for (i = 0; i < j; i++)
	        IntPoints[j] += Points[i];
	    IntPoints[j] /= n;
	}
    }

    return IntCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Converts a Bezier curve into Bspline curve by adding an open knot vector.  M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:     A Bezier curve to convert to a Bspline curve.                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A Bspline curve representing Bezier curve Crv.        M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtBezier2BsplineCrv, conversion                                       M
*****************************************************************************/
CagdCrvStruct *CnvrtBezier2BsplineCrv(CagdCrvStruct *Crv)
{
    CagdCrvStruct *BspCrv;

    if (Crv -> GType != CAGD_CBEZIER_TYPE) {
	CAGD_FATAL_ERROR(CAGD_ERR_WRONG_CRV);
	return NULL;
    }

    BspCrv = CagdCrvCopy(Crv);

    BspCrv -> Order = BspCrv -> Length;
    BspCrv -> KnotVector = BspKnotUniformOpen(BspCrv -> Length,
					           BspCrv -> Order, NULL);
    BspCrv -> GType = CAGD_CBSPLINE_TYPE;
    return BspCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Converts a Bspline curve into a set of Bezier curves by subdividing the    M
* Bspline curve at all its internal knots.				     M
*   Returned is a list of Bezier curves.				     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:     A Bspline curve to convert to a Bezier curve.                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   A list of Bezier curves representing the Bspline      M
*                      curve Crv.					     M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtBezier2BsplineCrv, conversion                                       M
*****************************************************************************/
CagdCrvStruct *CnvrtBspline2BezierCrv(CagdCrvStruct *Crv)
{
    CagdBType
	NewCrv = FALSE;
    int i, Order, Length;
    CagdRType LastT, *KnotVector;
    CagdCrvStruct *OrigCrv,
	*BezierCrvs = NULL;

    if (Crv -> GType != CAGD_CBSPLINE_TYPE) {
	CAGD_FATAL_ERROR(CAGD_ERR_WRONG_CRV);
	return NULL;
    }

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
        NewCrv = TRUE;
        Crv = CnvrtPeriodic2FloatCrv(Crv);
    }
    if (CAGD_IS_BSPLINE_CRV(Crv) && !BspCrvHasOpenEC(Crv)) {
	CagdCrvStruct
	    *TCrv = BspCrvOpenEnd(Crv);

	if (NewCrv)
	    CagdCrvFree(Crv);
	Crv = TCrv;
	NewCrv = TRUE;
    }

    Order = Crv -> Order,
    Length = Crv -> Length;
    KnotVector = Crv -> KnotVector;
    OrigCrv = Crv;

    for (i = Length - 1, LastT = KnotVector[Length]; i >= Order; i--) {
    	CagdRType
    	    t = KnotVector[i];
    	    
	if (!APX_EQ(LastT, t)) {
    	    CagdCrvStruct
		*Crvs = BspCrvSubdivAtParam(Crv, t);

    	    if (Crv != OrigCrv)
    	        CagdCrvFree(Crv);

    	    Crvs -> Pnext -> Pnext = BezierCrvs;
    	    BezierCrvs = Crvs -> Pnext;

    	    Crv = Crvs;
    	    Crv -> Pnext = NULL;

	    LastT = t;
    	}
    }

    if (Crv == OrigCrv) {
	/* No interior knots in this curve - just copy it: */
	BezierCrvs = CagdCrvCopy(Crv);
    }
    else {
    	Crv -> Pnext = BezierCrvs;
    	BezierCrvs = Crv;
    }

    for (Crv = BezierCrvs; Crv != NULL; Crv = Crv -> Pnext) {
    	Crv -> GType = CAGD_CBEZIER_TYPE;
	Crv -> Length = Crv -> Order;
	IritFree((VoidPtr) Crv -> KnotVector);
	Crv -> KnotVector = NULL;
    }

    if (NewCrv)
	CagdCrvFree(Crv);

    return BezierCrvs;
}
