home *** CD-ROM | disk | FTP | other *** search
/ The World of Computer Software / World_Of_Computer_Software-02-385-Vol-1of3.iso / i / iritsm3s.zip / cagd_lib / bspcoxdb.c < prev    next >
C/C++ Source or Header  |  1991-05-18  |  5KB  |  118 lines
  1. /******************************************************************************
  2. * BspCoxDB.c - Bspline evaluation using Cox - de Boor recursive algorithm.    *
  3. *******************************************************************************
  4. * Written by Gershon Elber, Aug. 90.                          *
  5. ******************************************************************************/
  6.  
  7. #ifdef __MSDOS__
  8. #include <stdlib.h>
  9. #endif /* __MSDOS__ */
  10.  
  11. #include <ctype.h>
  12. #include <stdio.h>
  13. #include <string.h>
  14. #include "cagd_loc.h"
  15.  
  16. /******************************************************************************
  17. * Returns a pointer to a static data, holding the value of the curve at given *
  18. * parametric location t. The curve is assumed to be non uniform spline.          *
  19. * Uses the recursive Cox de Boor algorithm, to evaluate the spline.          *
  20. ******************************************************************************/
  21. CagdRType *BspCrvEvalCoxDeBoor(CagdCrvStruct *Crv, CagdRType t)
  22. {
  23.     static CagdRType Pt[CAGD_MAX_PT_COORD];
  24.     CagdBType
  25.     IsNotRational = !CAGD_IS_RATIONAL_CRV(Crv);
  26.     CagdRType *p, *BasisFunc;
  27.     int i, j, l, IndexFirst,
  28.     k = Crv -> Order,
  29.     MaxCoord = CAGD_NUM_OF_PT_COORD(Crv -> PType);
  30.  
  31.     BasisFunc = BspCrvCoxDeBoorBasis(Crv -> KnotVector, k, Crv -> Length,
  32.                                      t, &IndexFirst);
  33.  
  34.     /* And finally multiply the basis functions with the control polygon. */
  35.     for (i = IsNotRational; i <= MaxCoord; i++) {
  36.     Pt[i] = 0;
  37.     p = Crv -> Points[i];
  38.     for (j = IndexFirst, l = 0; l < k; )
  39.         Pt[i] += p[j++] * BasisFunc[l++];
  40.     }
  41.  
  42.     return Pt;
  43. }
  44.  
  45. /******************************************************************************
  46. * Returns a pointer to a vector of size order, holding values of the non zero *
  47. * basis functions of a given curve at given parametric location t.          *
  48. * This vector SHOULD NOT BE FREED. Although it is dynamically allocated, the  *
  49. * returned pointer does not point to the beginning of this memory and it will *
  50. * be maintained by this routine (i.e. it will be freed next time this routine *
  51. * is called).                                      *
  52. * IndexFirst returns the index of first non zero basis function for given t.  *
  53. * The curve is assumed to be non uniform bspline.                  *
  54. * Uses the recursive Cox de Boor algorithm, to evaluate the spline.          *
  55. * Algorithm:                                      *
  56. * Use the following recursion relation with B(i,0) == 1.                      *
  57. *                                          *
  58. *          t     -    t(i)            t(i+k)    -   t                         *
  59. * B(i,k) = --------------- B(i,k-1) + --------------- B(i+1,k-1)              *
  60. *          t(i+k-1) - t(i)            t(i+k) - t(i+1)                         *
  61. *                                          *
  62. * Starting with constant spline (k == 1) only one basis func. is non zero and *
  63. * equal to one. This is the constant spline spanning interval t(i)..t(i+1)    *
  64. * such that t(i) <= t and t(i+1) > t. We then raise this constant spline      *
  65. * to the Crv order and finding in this process all the basis functions that   *
  66. * are non zero in t for order Order. Sound limple hah!?                  *
  67. ******************************************************************************/
  68. CagdRType *BspCrvCoxDeBoorBasis(CagdRType *KnotVector, int Order, int Len,
  69.                         CagdRType t, int *IndexFirst)
  70. {
  71.     static CagdRType *B = NULL;
  72.     CagdRType s1, s2, *BasisFunc;
  73.     int i, l,
  74.     KVLen = Order + Len,
  75.     Index = BspKnotLastIndexLE(KnotVector, KVLen, t);
  76.  
  77.     if (!BspKnotParamInDomain(KnotVector, Len, Order, t))
  78.     FATAL_ERROR(CAGD_ERR_T_NOT_IN_CRV);
  79.  
  80.     /* Starting the recursion from constant splines - one spline is non     */
  81.     /* zero and is equal to one. This is the spline that starts at Index.   */
  82.     /* As We are going to reference index -1 we increment the buffer by one */
  83.     /* and save 0.0 at index -1. We then initialize the constant spline     */
  84.     /* values - all are zero but the one from t(i) to t(i+1).               */
  85.     if (B != NULL) CagdFree((VoidPtr) B);
  86.     BasisFunc = B = (CagdRType *) CagdMalloc(sizeof(CagdRType) * (1 + Order));
  87.     *BasisFunc++ = 0.0;
  88.     if (Index >= Len + Order - 1) {
  89.     /* We are at the end of the parametric domain and this is open      */
  90.     /* end comdition - simply return last point.                */
  91.     for (i = 0; i < Order; i++)
  92.         BasisFunc[i] = (CagdRType) (i == Order - 1);
  93.  
  94.     *IndexFirst = Len - Order;
  95.     return BasisFunc;
  96.     }
  97.     else
  98.     for (i = 0; i < Order; i++) BasisFunc[i] = (CagdRType) (i == 0);
  99.  
  100.     /* Here is the tricky part. we raise these constant splines to the      */
  101.     /* required order of the curve Crv for the given parameter t. There are */
  102.     /* at most order non zero function at param. value t. These functions   */
  103.     /* start at Index-order+1 up to Index (order functions overwhole).      */
  104.     for (i = 2; i <= Order; i++) {            /* Goes through all orders... */
  105.     for (l = i - 1; l >= 0; l--) {  /* And all basis funcs. of order i. */
  106.         s1 = (KnotVector[Index + l] - KnotVector[Index + l - i + 1]);
  107.         s1 = APX_EQ(s1, 0.0) ? 0.0 : (t - KnotVector[Index + l - i + 1]) / s1;
  108.         s2 = (KnotVector[Index + l + 1] - KnotVector[Index + l - i + 2]);
  109.         s2 = APX_EQ(s2, 0.0) ? 0.0 : (KnotVector[Index + l + 1] - t) / s2;
  110.  
  111.         BasisFunc[l] = s1 * BasisFunc[l - 1] + s2 * BasisFunc[l];
  112.     }
  113.     }
  114.  
  115.     *IndexFirst = Index - Order + 1;
  116.     return BasisFunc;
  117. }
  118.