/******************************************************************************
* BspCoxDB.c - Bspline evaluation using Cox - de Boor recursive algorithm.    *
*******************************************************************************
* Written by Gershon Elber, Aug. 90.					      *
******************************************************************************/

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

/*****************************************************************************
* DESCRIPTION:                                                               M
* Returns a pointer to a static data, holding the value of the curve at      M
* the prescribed parametric location t.                                      M
*   Uses the recursive Cox de-Boor algorithm, to evaluate the spline, which  M
* is not very efficient if many evaluations of the same curve are necessary  M
*   Use knot insertion when multiple evaluations are to be performed.        M
*                                                                            *
* PARAMETERS:                                                                M
*   Crv:      To evaluate at the given parametric location t.                M
*   t:        The parameter value at which the curve Crv is to be evaluated. M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:  A vector holding all the coefficients of all components    M
*                 of curve Crv's point type. If for example the curve's      M
*                 point type is P2, the W, X, and Y will be saved in the     M
*                 first three locations of the returned vector. The first    M
*                 location (index 0) of the returned vector is reserved for  M
*                 the rational coefficient W and XYZ always starts at second M
*                 location of the returned vector (index 1).                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvEvalCoxDeBoor, evaluation, Bsplines                                M
*****************************************************************************/
CagdRType *BspCrvEvalCoxDeBoor(CagdCrvStruct *Crv, CagdRType t)
{
    static CagdRType Pt[CAGD_MAX_PT_COORD];
    CagdBType
	IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
    CagdRType *p, *BasisFunc;
    int i, j, l, IndexFirst,
	k = Crv -> Order,
	Length = Crv -> Length,
	MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);

    BasisFunc = BspCrvCoxDeBoorBasis(Crv -> KnotVector, k,
				     CAGD_CRV_PT_LST_LEN(Crv),
				     t, &IndexFirst);

    /* And finally multiply the basis functions with the control polygon. */
    for (i = IsNotRational; i <= MaxCoord; i++) {
	Pt[i] = 0;
	p = Crv -> Points[i];
	for (j = IndexFirst, l = 0; l < k; )
	    Pt[i] += p[j++ % Length] * BasisFunc[l++];
    }

    return Pt;
}

/******************************************************************************
* DESCRIPTION:                                                               M
* Returns a pointer to a vector of size Order, holding values of the non     M
* zero basis functions of a given curve at given parametric location t.	     M
*   This vector SHOULD NOT BE FREED. Although it is dynamically allocated,   M
* the returned pointer does not point to the beginning of this memory and it M
* it be maintained by this routine (i.e. it might be freed next time this    M
* routine is being called).						     M
*   IndexFirst returns the index of first non zero basis function for the    M
* given parameter value t.						     M
*   Uses the recursive Cox de Boor algorithm, to evaluate the Bspline basis  M
* functions.								     M
*   Algorithm:								     M
* Use the following recursion relation with B(i,0) == 1.                     M
*									     M
*          t     -    t(i)            t(i+k)    -   t                        V
* B(i,k) = --------------- B(i,k-1) + --------------- B(i+1,k-1)             V
*          t(i+k-1) - t(i)            t(i+k) - t(i+1)                        V
*									     M
*   Starting with constant Bspline (k == 0) only one basis function is non   M
* zero and is equal to one. This is the constant Bspline spanning interval   M
* t(i)...t(i+1) such that t(i) <= t < t(i+1). We then raise this constant    M
* Bspline to the prescribed Order and find in this process all the basis     M
* functions that are non zero in t for order Order. Sound simple hah!?	     M
*                                                                            *
* PARAMETERS:                                                                M
*   KnotVector:    To evaluate the Bspline Basis functions for this space.   M
*   Order:         Of the geometry.                                          M
*   Len:           Number of control points in the geometry. The length of   M
*                  KnotVector is equal to Len + Order.                       M
*   t:             At which the Bspline basis functions are to be evaluated. M
*   IndexFirst:    Index of the first Bspline basis function that might be   M
*                  non zero.                                                 M
*                                                                            *
* RETURN VALUE:                                                              M
*   CagdRType *:   A vector of length Order thats holds the values of the    M
*                  Bspline basis functions for the given t. A Bspline of     M
*                  order Order might have at most Order non zero basis       M
*                  functions that will hence start at IndexFirst and upto    M
*                  (*IndexFirst + Order - 1).                                M
*                                                                            *
* KEYWORDS:                                                                  M
*   BspCrvCoxDeBoorBasis, evaluation, Bsplines                               M
*****************************************************************************/
CagdRType *BspCrvCoxDeBoorBasis(CagdRType *KnotVector,
				int Order,
				int Len,
				CagdRType t,
				int *IndexFirst)
{
    static CagdRType
	*B = NULL;
    CagdRType s1, s2, *BasisFunc;
    int i, l, Index,
	KVLen = Order + Len;

    if (!BspKnotParamInDomain(KnotVector, Len, Order, FALSE, t))
	CAGD_FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
    if (t == KnotVector[Len])
	t -= IRIT_EPSILON;
    Index = BspKnotLastIndexLE(KnotVector, KVLen, t);

    /* Starting the recursion from constant splines - one spline is non     */
    /* zero and is equal to one. This is the spline that starts at Index.   */
    /* As We are going to reference index -1 we increment the buffer by one */
    /* and save 0.0 at index -1. We then initialize the constant spline     */
    /* values - all are zero but the one from t(i) to t(i+1).               */
    if (B != NULL)
	IritFree((VoidPtr) B);
    BasisFunc = B = (CagdRType *) IritMalloc(sizeof(CagdRType) * (1 + Order));
    *BasisFunc++ = 0.0;
    if (Index >= Len + Order - 1) {
	/* We are at the end of the parametric domain and this is open      */
	/* end condition - simply return last point.			    */
	for (i = 0; i < Order; i++)
	    BasisFunc[i] = (CagdRType) (i == Order - 1);

	*IndexFirst = Len - Order;
	return BasisFunc;
    }
    else
	for (i = 0; i < Order; i++)
	    BasisFunc[i] = (CagdRType) (i == 0);

    /* Here is the tricky part. we raise these constant splines to the      */
    /* required order of the curve Crv for the given parameter t. There are */
    /* at most order non zero function at param. value t. These functions   */
    /* start at Index-order+1 up to Index (order functions overwhole).      */
    for (i = 2; i <= Order; i++) {            /* Goes through all orders... */
	for (l = i - 1; l >= 0; l--) {  /* And all basis funcs. of order i. */
	    s1 = (KnotVector[Index + l] - KnotVector[Index + l - i + 1]);
	    s1 = APX_EQ(s1, 0.0) ? 0.0 : (t - KnotVector[Index + l - i + 1]) / s1;
	    s2 = (KnotVector[Index + l + 1] - KnotVector[Index + l - i + 2]);
	    s2 = APX_EQ(s2, 0.0) ? 0.0 : (KnotVector[Index + l + 1] - t) / s2;

	    BasisFunc[l] = s1 * BasisFunc[l - 1] + s2 * BasisFunc[l];
	}
    }

    *IndexFirst = Index - Order + 1;
    return BasisFunc;
}
