/******************************************************************************
* CBsp-Aux.c - Bspline curve auxilary routines.				      *
*******************************************************************************
* Written by Gershon Elber, Aug. 90.					      *
******************************************************************************/

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

/*****************************************************************************
* DESCRIPTION:                                                               M
* Given a Bspline 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
*   The subdivision is achieved by inserting (order-1) knot at the given     M
* parameter value t and spliting the control polygon and knot vector at that M
* location.                                                                  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
*   BspCrvSubdivAtParam, subdivision, refinement                             M
*****************************************************************************/
CagdCrvStruct *BspCrvSubdivAtParam(CagdCrvStruct *Crv, CagdRType t)
{
    CagdBType
	NewCrv = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, Len, KVLen, Index1, Index2, Mult,
	k = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType TMin, TMax, *LOnePts, *ROnePts, *OnePts, *NewKV,
	**Pts, **LPts, **RPts;
    CagdCrvStruct *LCrv, *RCrv;
    BspKnotAlphaCoeffType *A;

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
	NewCrv = TRUE;
	Crv = CnvrtPeriodic2FloatCrv(Crv);
    }

    Len = Crv -> Length;
    KVLen = k + Len;
    Index1 = BspKnotLastIndexL(Crv -> KnotVector, KVLen, t);
    if (Index1 + 1 < k)
	Index1 = k - 1;
    Index2 = BspKnotFirstIndexG(Crv -> KnotVector, KVLen, t);
    if (Index2 > Len)
	Index2 = Len;
    Mult = k - 1 - (Index2 - Index1 - 1);

    CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, Crv, TMin, TMax);

    LCrv = BspCrvNew(Index1 + 1, k, Crv -> PType);
    RCrv = BspCrvNew(Len - Index2 + k, k, Crv -> PType);

    /* Update the new knot vectors. */
    CAGD_GEN_COPY(LCrv -> KnotVector,
	     Crv -> KnotVector,
	     sizeof(CagdRType) * (Index1 + 1));
    /* Close the knot vector with multiplicity Order: */
    for (j = Index1 + 1; j <= Index1 + k; j++)
	LCrv -> KnotVector[j] = t;
    CAGD_GEN_COPY(&RCrv -> KnotVector[k],
		  &Crv -> KnotVector[Index2],
		  sizeof(CagdRType) * (Len + k - Index2));
    /* Make sure knot vector starts with multiplicity Order: */
    for (j = 0; j < k; j++)
	RCrv -> KnotVector[j] = t;

    /* Now handle the control polygon refinement. */
    Pts = Crv -> Points,
    LPts = LCrv -> Points,
    RPts = RCrv -> Points;

    if (Mult > 0) {
	NewKV = IritMalloc(sizeof(CagdRType) * Mult);
	for (i = 0; i < Mult; i++)
	    NewKV[i] = (t == TMax ? t - CAGD_DOMAIN_EPSILON : t);
	A = BspKnotEvalAlphaCoefMerge(k, Crv -> KnotVector, Len, NewKV,
				      Mult);
	IritFree((VoidPtr) NewKV);
    }
    else {
	A = BspKnotEvalAlphaCoef(k, Crv -> KnotVector, Len,
				    Crv -> KnotVector, Len);
    }

    /* Note that Mult can be negative in cases where original       */
    /* multiplicity was order or more and we need to compensate     */
    /* here, since Alpha matrix will be just a unit matrix then.    */
    Mult = Mult >= 0 ? 0 : -Mult;

    /* Blend Crv into LCrv. */
    for (j = IsNotRational; j <= MaxCoord; j++) {
	LOnePts = LPts[j];
	OnePts = Pts[j];
	for (i = 0; i < LCrv -> Length; i++, LOnePts++)
	    CAGD_ALPHA_BLEND(A, i, OnePts, Crv -> Length, LOnePts);
    }

    /* Blend Crv into RCrv. */
    for (j = IsNotRational; j <= MaxCoord; j++) {
	ROnePts = RPts[j];
	OnePts = Pts[j];
	for (i = LCrv -> Length - 1 + Mult;
	     i < LCrv -> Length + RCrv -> Length - 1 + Mult;
	     i++, ROnePts++)
	    CAGD_ALPHA_BLEND(A, i, OnePts, Crv -> Length, ROnePts);
    }

    BspKnotFreeAlphaCoef(A);

    BspKnotMakeRobustKV(RCrv -> KnotVector,
			RCrv -> Order + RCrv -> Length);
    BspKnotMakeRobustKV(LCrv -> KnotVector,
			LCrv -> Order + LCrv -> Length);

    LCrv -> Pnext = RCrv;

    if (NewCrv)
	CagdCrvFree(Crv);

    return LCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Inserts n knot, all with the value t. In no case will the multiplicity of M
* a knot be greater or equal to the curve order.			     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To refine by insertion (upto) n knot of value t.             M
*   t:          Parameter value of new knot to insert.                       M
*   n:          Maximum number of times t should be inserted.                M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Refined Crv with n knots of value t.                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvKnotInsertNSame, refinement, subdivision                           M
*****************************************************************************/
CagdCrvStruct *BspCrvKnotInsertNSame(CagdCrvStruct *Crv, CagdRType t, int n)
{
    int i,
	CrntMult = BspKnotFindMult(Crv -> KnotVector, Crv -> Order,
							   Crv -> Length, t),
	Mult = MIN(n, Crv -> Order - CrntMult - 1);
    CagdCrvStruct *RefinedCrv;

    if (Mult > 0) {
	CagdRType
	    *NewKV = (CagdRType *) IritMalloc(sizeof(CagdRType) * Mult);

	for (i = 0; i < Mult; i++)
	    NewKV[i] = t;

	RefinedCrv = BspCrvKnotInsertNDiff(Crv, FALSE, NewKV, Mult);

	IritFree((VoidPtr) NewKV);
    }
    else {
	RefinedCrv = CagdCrvCopy(Crv);
    }

    return RefinedCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Inserts n knot with different values as defined by the vector t. If,      M
* however, Replace is TRUE, the knot are simply replacing the current knot   M
* vector.								     M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To refine by insertion (upto) n knot of value t.             M
*   Replace:    if TRUE, the n knots in t should replace the knot vector     M
*               of size n of Crv. Sizes must match. If False, n new knots    M
*               as defined by t will be introduced into Crv.                 M
*   t:          New knots to introduce/replace knot vector of Crv.           M
*   n:          Size of t.                                                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Refined Crv with n new knots.                         M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvKnotInsertNDiff, refinement, subdivision                           M
*****************************************************************************/
CagdCrvStruct *BspCrvKnotInsertNDiff(CagdCrvStruct *Crv,
				     CagdBType Replace,
				     CagdRType *t, int n)
{
    CagdBType
	NewCrv = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    CagdRType *KnotVector;
    int Length,
	Order = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *RefinedCrv;

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
	NewCrv = TRUE;
	Crv = CnvrtPeriodic2FloatCrv(Crv);
    }

    Length = Crv -> Length;
    KnotVector = Crv -> KnotVector;

    if (Replace) {
	int i;

	if (Order + Length != n)
	    CAGD_FATAL_ERROR(CAGD_ERR_NUM_KNOT_MISMATCH);
	for (i = 1; i < n; i++)
	    if (t[i] < t[i - 1])
		CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);

    	RefinedCrv = CagdCrvCopy(Crv);
	for (i = 0; i < n; i++)
	    RefinedCrv -> KnotVector[i] = *t++;
    }
    else if (n == 0) {
	RefinedCrv = CagdCrvCopy(Crv);
    }
    else {
	int i, j, LengthKVt;
	BspKnotAlphaCoeffType *A;
	CagdRType *MergedKVt, TMin, TMax;

	BspCrvDomain(Crv, &TMin, &TMax);

	for (i = 1; i < n; i++)
	    if (t[i] < t[i - 1])
		CAGD_FATAL_ERROR(CAGD_ERR_KNOT_NOT_ORDERED);
	for (i = 0; i < n; i++)
	    if (t[i] >= TMax)
	        t[i] -= CAGD_DOMAIN_EPSILON;

	/* Compute the Alpha refinement matrix. */
	MergedKVt = BspKnotMergeTwo(KnotVector, Length + Order,
				    t, n, 0, &LengthKVt);
	A = BspKnotEvalAlphaCoef(Order, KnotVector, Length,
				 MergedKVt, LengthKVt - Order);

	RefinedCrv = BspCrvNew(Crv -> Length + n, Order, Crv -> PType);

	/* Update the knot vector. */
	IritFree((VoidPtr) RefinedCrv -> KnotVector);
	RefinedCrv -> KnotVector = MergedKVt;

	/* Update the control mesh */
	for (j = IsNotRational; j <= MaxCoord; j++) {
	    CagdRType
		*ROnePts = RefinedCrv -> Points[j],
		*OnePts = Crv -> Points[j];
	    for (i = 0; i < RefinedCrv -> Length; i++, ROnePts++)
		CAGD_ALPHA_BLEND(A, i, OnePts, Crv -> Length, ROnePts);
	}
	BspKnotFreeAlphaCoef(A);
    }

    BspKnotMakeRobustKV(RefinedCrv -> KnotVector,
			RefinedCrv -> Order + RefinedCrv -> Length);

    if (NewCrv)
	CagdCrvFree(Crv);

    return RefinedCrv;
}

/*****************************************************************************
* 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
*   BspCrvDegreeRaiseN, degree raising                                       M
*****************************************************************************/
CagdCrvStruct *BspCrvDegreeRaiseN(CagdCrvStruct *Crv, int NewOrder)
{
    CagdBType
	NewCrv = FALSE;
    int i, j, RaisedOrder, Length, KvLen1,
	Order = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *RaisedCrv, *UnitCrv;
    CagdRType *Kv;

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
	NewCrv = TRUE;
	Crv = CnvrtPeriodic2FloatCrv(Crv);
    }

    Length = Crv -> Length;
    KvLen1 = Order + Length - 1;
    Kv = Crv -> KnotVector;

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

    UnitCrv = BspCrvNew(RaisedOrder, RaisedOrder,
			CAGD_MAKE_PT_TYPE(FALSE, MaxCoord));
    for (i = 0; i < RaisedOrder * 2; i++)
	UnitCrv -> KnotVector[i] = i >= RaisedOrder ? Kv[KvLen1] : Kv[0];
    for (i = 1; i <= MaxCoord; i++)
	for (j = 0; j < RaisedOrder; j++)
	    UnitCrv -> Points[i][j] = 1.0;

    RaisedCrv = BspCrvMult(Crv, UnitCrv);

    CagdCrvFree(UnitCrv);

    if (NewCrv)
	CagdCrvFree(Crv);

    return RaisedCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new curve, identical to the original but with one degree higher. M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To raise it degree by one.                                   M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A curve with one degree higher representing the same   M
*                     geometry as Crv.                                       M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvDegreeRaise, degree raising                                        M
*****************************************************************************/
CagdCrvStruct *BspCrvDegreeRaise(CagdCrvStruct *Crv)
{
    CagdBType
	NewCrv = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, i2, j, RaisedLen, Length,
	Order = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *RaisedCrv;

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
	NewCrv = TRUE;
	Crv = CnvrtPeriodic2FloatCrv(Crv);
    }

    Length = Crv -> Length;

    if (Order > 2)
	return BspCrvDegreeRaiseN(Crv, Order + 1);

    /* If curve is linear, degree raising means basically to increase the    */
    /* knot multiplicity of each segment by one and add a middle point for   */
    /* each such segment.						     */
    RaisedLen = Length * 2 - 1;
    RaisedCrv = BspCrvNew(RaisedLen, Order + 1, Crv -> PType);

    /* Update the control polygon; */
    for (j = IsNotRational; j <= MaxCoord; j++)		     /* First point. */
	RaisedCrv -> Points[j][0] = Crv -> Points[j][0];

    for (i = 1, i2 = 1; i < Length; i++, i2 += 2)
	for (j = IsNotRational; j <= MaxCoord; j++) {
	    RaisedCrv -> Points[j][i2] =
		Crv -> Points[j][i-1] * 0.5 + Crv -> Points[j][i] * 0.5;
	    RaisedCrv -> Points[j][i2 + 1] = Crv -> Points[j][i];
	}

    /* Update the knot vector. */
    for (i = 0; i < 3; i++)
	RaisedCrv -> KnotVector[i] = Crv -> KnotVector[0];

    for (i = 2, j = 3; i < Length; i++, j += 2)
	RaisedCrv -> KnotVector[j] = RaisedCrv -> KnotVector[j + 1] = 
	    Crv -> KnotVector[i];

    for (i = j; i < j + 3; i++)
	RaisedCrv -> KnotVector[i] = Crv -> KnotVector[Length];

    if (NewCrv)
	CagdCrvFree(Crv);

    return RaisedCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a unit vector, equal to the tangent to Crv at parameter value t.   M
*   Algorithm: insert (order - 1) knots and return control polygon tangent.  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
*   BspCrvTangent, tangent                                                   M
*****************************************************************************/
CagdVecStruct *BspCrvTangent(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdVecStruct P2;
    CagdVecStruct P1, *T;
    CagdRType TMin, TMax;
    CagdCrvStruct
	*FCrv = CAGD_IS_PERIODIC_CRV(Crv) ? CnvrtPeriodic2FloatCrv(Crv) : Crv;
    int k = FCrv -> Order,
	Len = FCrv -> Length,
	OpenEnd = BspCrvHasOpenEC(FCrv),
	Index = BspKnotLastIndexL(FCrv -> KnotVector, k + Len, t);
    CagdPointType
	PType = FCrv -> PType;
    CagdCrvStruct *RefinedCrv;

    CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, FCrv, TMin, TMax);

    if (APX_EQ(t, TMin) && OpenEnd) {
	/* Use Crv starting tangent direction. */
	CagdCoerceToE3(P1.Vec, FCrv -> Points, 0, PType);
	CagdCoerceToE3(P2.Vec, FCrv -> Points, 1, PType);
    }
    else if (APX_EQ(t, TMax) && OpenEnd) {
	/* Use Crv ending tangent direction. */
	CagdCoerceToE3(P1.Vec, FCrv -> Points, Len - 2, PType);
	CagdCoerceToE3(P2.Vec, FCrv -> Points, Len - 1, PType);
    }
    else {
	RefinedCrv = BspCrvKnotInsertNSame(FCrv, t, k - 1);

	CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
	CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);

	CagdCrvFree(RefinedCrv);
    }

    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 = BspCrvTangent(Crv, t - TMin < TMax - t ? t + EPSILON
						       : t - EPSILON);

	    AttrFreeOneAttribute(&Crv -> Attr, "_tan");
	    if (FCrv != Crv)
		CagdCrvFree(FCrv);
	    return T;
	}
	else {
	    /* A zero length vector signals failure to compute tangent. */
	    if (FCrv != Crv)
		CagdCrvFree(FCrv);
	    return &P2;
	}
    }
    else {
	CAGD_NORMALIZE_VECTOR(P2);		    /* Normalize the vector. */
	if (FCrv != Crv)
	    CagdCrvFree(FCrv);
	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
*   BspCrvBiNormal, binormal                                                 M
*****************************************************************************/
CagdVecStruct *BspCrvBiNormal(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdVecStruct P3;
    CagdVecStruct P1, P2;
    CagdRType TMin, TMax;
    CagdCrvStruct
	*FCrv = CAGD_IS_PERIODIC_CRV(Crv) ? CnvrtPeriodic2FloatCrv(Crv) : Crv;
    int k = FCrv -> Order,
	Len = FCrv -> Length,
	OpenEnd = BspCrvHasOpenEC(FCrv),
	Index = BspKnotLastIndexL(FCrv -> KnotVector, k + Len, t);
    CagdPointType
	PType = FCrv -> PType;
    CagdCrvStruct *RefinedCrv;

    CAGD_DOMAIN_GET_AND_VERIFY_CRV(t, FCrv, TMin, TMax);

    /* 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(FCrv -> PType) == 2) {
	P3.Vec[0] = P3.Vec[1] = 0.0;
	P3.Vec[2] = 1.0;
	return &P3;
    }

    if (APX_EQ(t, TMin) && OpenEnd) {
	/* Use Crv starting tangent direction. */
	CagdCoerceToE3(P1.Vec, FCrv -> Points, 0, PType);
	CagdCoerceToE3(P2.Vec, FCrv -> Points, 1, PType);
	CagdCoerceToE3(P3.Vec, FCrv -> Points, 2, PType);
    }
    else if (APX_EQ(t, TMax) && OpenEnd) {
	/* Use Crv ending tangent direction. */
	CagdCoerceToE3(P1.Vec, FCrv -> Points, Len - 3, PType);
	CagdCoerceToE3(P2.Vec, FCrv -> Points, Len - 2, PType);
	CagdCoerceToE3(P3.Vec, FCrv -> Points, Len - 1, PType);
    }
    else {
	RefinedCrv = BspCrvKnotInsertNSame(FCrv, t, k - 1);

	CagdCoerceToE3(P1.Vec, RefinedCrv -> Points, Index, PType);
	CagdCoerceToE3(P2.Vec, RefinedCrv -> Points, Index + 1, PType);
	CagdCoerceToE3(P3.Vec, RefinedCrv -> Points, Index + 2, PType);

	CagdCrvFree(RefinedCrv);
    }

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

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

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

    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
*   BspCrvNoraml, normal                                                     M
*****************************************************************************/
CagdVecStruct *BspCrvNormal(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdVecStruct N, *T, *B;

    T = BspCrvTangent(Crv, t);
    B = BspCrvBiNormal(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)) / (Kv(i + k) - Kv(i + 1)), i = 0 to k-2.  V
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:        To differentiate.                                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:   Differentiated curve.                                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvDerive, derivatives                                                M
*****************************************************************************/
CagdCrvStruct *BspCrvDerive(CagdCrvStruct *Crv)
{
    CagdBType
	NewCrv = FALSE,
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    int i, j, Len,
	k = Crv -> Order,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdRType *Kv;
    CagdCrvStruct *DerivedCrv;

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
	NewCrv = TRUE;
	Crv = CnvrtPeriodic2FloatCrv(Crv);
    }

    if (!IsNotRational) {
	DerivedCrv = BspCrvDeriveRational(Crv);
	if (NewCrv)
	    CagdCrvFree(Crv);
	return DerivedCrv;
    }

    Len = Crv -> Length;
    Kv = Crv -> KnotVector;
    DerivedCrv = BspCrvNew(MAX(Len - 1, 1), MAX(k - 1, 1), Crv -> PType);

    if (k >= 2) {
	for (i = 0; i < Len - 1; i++) {
	    CagdRType
		Denom = Kv[i + k] - Kv[i + 1];

	    for (j = IsNotRational; j <= MaxCoord; j++)
		DerivedCrv -> Points[j][i] = (k - 1) *
		    (Crv -> Points[j][i + 1] - Crv -> Points[j][i]) / Denom;
	}
    }
    else {
	for (i = 0; i < MAX(Len - 1, 1); i++)
	    for (j = IsNotRational; j <= MaxCoord; j++)
		DerivedCrv -> Points[j][i] = 0.0;
		 
    }

    CAGD_GEN_COPY(DerivedCrv -> KnotVector,
		  &Crv -> KnotVector[k < 2 ? 0 : 1],
		  sizeof(CagdRType) * (MAX(k - 1, 1) + MAX(Len - 1, 1)));

    if (NewCrv)
	CagdCrvFree(Crv);

    return DerivedCrv;
}

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

    if (CAGD_IS_PERIODIC_CRV(Crv)) {
	NewCrv = TRUE;
	Crv = CnvrtPeriodic2FloatCrv(Crv);
    }

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

    Len = Crv -> Length;
    Kv = Crv -> KnotVector;

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

    /* Copy the knot vector and duplicate the two end knots. */
    CAGD_GEN_COPY(&IntCrv -> KnotVector[1], Kv, sizeof(CagdRType) * (Len + n));
    IntCrv -> KnotVector[0] = Kv[0];
    IntCrv -> KnotVector[Len + n + 1] = Kv[Len + n - 1];
    Kv = IntCrv -> KnotVector;

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

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

    if (NewCrv)
	CagdCrvFree(Crv);

    return IntCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a new linear Bspline curve constructed from the given polyline.    M
*                                                                            *
* PARAMETERS:                                                                M
*   Poly:     To convert to a linear bspline curve.                          M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A linear Bspline curve representing Poly.              M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtPolyline2LinBsplineCrv, linear curves, conversion                   M
*****************************************************************************/
CagdCrvStruct *CnvrtPolyline2LinBsplineCrv(CagdPolylineStruct *Poly)
{
    int i,
	Length = Poly -> Length;
    CagdCrvStruct
	*Crv = BspCrvNew(Length, 2, CAGD_PT_E3_TYPE);
    CagdRType
	**Points = Crv -> Points;
    CagdPolylnStruct
	*Pts = Poly -> Polyline;

    BspKnotUniformOpen(Length, 2, Crv -> KnotVector);

    for (i = 0; i < Length; i++, Pts++) {
	Points[1][i] = Pts -> Pt[0];
	Points[2][i] = Pts -> Pt[1];
	Points[3][i] = Pts -> Pt[2];
    }

    return Crv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Converts a Bspline curve to a Bspline curve with floating end conditions.  M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       Bspline curve to convert to floating end conditions. Assume   M
*              Crv is either periodic or has floating end condition.         M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A Bspline curve with floating end conditions,          M
*                     representing the same geometry as Crv.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtPeriodic2FloatCrv, conversion                                       M
*****************************************************************************/
CagdCrvStruct *CnvrtPeriodic2FloatCrv(CagdCrvStruct *Crv)
{
    int i,
	Order = Crv -> Order,
	Length = Crv -> Length,
	MaxAxis = CAGD_NUM_OF_PT_COORD(Crv -> PType);
    CagdCrvStruct *NewCrv;

    if (!CAGD_IS_PERIODIC_CRV(Crv)) {
	CAGD_FATAL_ERROR(CAGD_ERR_PERIODIC_EXPECTED);
	return NULL;
    }

    NewCrv = BspCrvNew(Length + Order - 1, Order, Crv -> PType);

    NewCrv -> KnotVector = BspKnotCopy(Crv -> KnotVector,
				       Length + Order + Order - 1);

    for (i = !CAGD_IS_RATIONAL_PT(Crv -> PType); i <= MaxAxis; i++) {
	NewCrv -> Points[i] = (CagdRType *) IritMalloc(sizeof(CagdRType) *
					       (Length + Order + Order - 1));
	CAGD_GEN_COPY(NewCrv -> Points[i], Crv -> Points[i],
		      sizeof(CagdRType) * Length);
	CAGD_GEN_COPY(&NewCrv -> Points[i][Length], Crv -> Points[i],
		      sizeof(CagdRType) * (Order - 1));
    }

    for (i = MaxAxis + 1; i <= CAGD_MAX_PT_COORD; i++)
	NewCrv -> Points[i] = NULL;
    
    return NewCrv;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Converts a Bspline curve to a Bspline curve with open end conditions.      M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:       Bspline curve to convert to open end conditions.		     M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdCrvStruct *:  A Bspline curve with open end conditions,	             M
*                     representing the same geometry as Crv.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   CnvrtPeriodic2FloatCrv, conversion                                       M
*****************************************************************************/
CagdCrvStruct *CnvrtFloat2OpenCrv(CagdCrvStruct *Crv)
{
    CagdRType TMin, TMax;

    if (!CAGD_IS_BSPLINE_CRV(Crv)) {
	CAGD_FATAL_ERROR(CAGD_ERR_BSP_CRV_EXPECT);
	return NULL;
    }

    CagdCrvDomain(Crv, &TMin, &TMax);
    return CagdCrvRegionFromCrv(Crv, TMin, TMax);
}

