ProfitPress Mega CDROM2 …eeware (MSDOS)(1992)(Eng)/ProfitPress-MegaCDROM2.B6I/GRAPHICS/PLOT/IRITSM3S.ZIP/CAGD_LIB/BSPCOXDB.C

C/C++ Source or Header | 1991-05-18 | 5.3 KB | 118 lines

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

#ifdef __MSDOS__
#include <stdlib.h>
#endif /* __MSDOS__ */

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

/******************************************************************************
* Returns a pointer to a static data, holding the value of the curve at given *
* parametric location t. The curve is assumed to be non uniform spline.          *
* Uses the recursive Cox de Boor algorithm, to evaluate the spline.          *
******************************************************************************/
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,
    MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);

    BasisFunc = BspCrvCoxDeBoorBasis(Crv -> KnotVector, k, Crv -> Length,
                                     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++] * BasisFunc[l++];
    }

    return Pt;
}

/******************************************************************************
* Returns a pointer to a vector of size order, holding values of the non zero *
* basis functions of a given curve at given parametric location t.          *
* This vector SHOULD NOT BE FREED. Although it is dynamically allocated, the  *
* returned pointer does not point to the beginning of this memory and it will *
* be maintained by this routine (i.e. it will be freed next time this routine *
* is called).                                      *
* IndexFirst returns the index of first non zero basis function for given t.  *
* The curve is assumed to be non uniform bspline.                  *
* Uses the recursive Cox de Boor algorithm, to evaluate the spline.          *
* Algorithm:                                      *
* Use the following recursion relation with B(i,0) == 1.                      *
*                                          *
*          t     -    t(i)            t(i+k)    -   t                         *
* B(i,k) = --------------- B(i,k-1) + --------------- B(i+1,k-1)              *
*          t(i+k-1) - t(i)            t(i+k) - t(i+1)                         *
*                                          *
* Starting with constant spline (k == 1) only one basis func. is non zero and *
* equal to one. This is the constant spline spanning interval t(i)..t(i+1)    *
* such that t(i) <= t and t(i+1) > t. We then raise this constant spline      *
* to the Crv order and finding in this process all the basis functions that   *
* are non zero in t for order Order. Sound limple hah!?                  *
******************************************************************************/
CagdRType *BspCrvCoxDeBoorBasis(CagdRType *KnotVector, int Order, int Len,
                        CagdRType t, int *IndexFirst)
{
    static CagdRType *B = NULL;
    CagdRType s1, s2, *BasisFunc;
    int i, l,
    KVLen = Order + Len,
    Index = BspKnotLastIndexLE(KnotVector, KVLen, t);

    if (!BspKnotParamInDomain(KnotVector, Len, Order, t))
    FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);

    /* 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) CagdFree((VoidPtr) B);
    BasisFunc = B = (CagdRType *) CagdMalloc(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 comdition - 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;
}