home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch4_9 / lincrv.c next >
Encoding:
C/C++ Source or Header  |  1995-04-04  |  2.3 KB  |  63 lines
  1. /****** lincrv.c ******/
  2. /* Ken Shoemake, 1994 */
  3.  
  4. #include "lincrv.h"
  5.  
  6. /* Perform a generic vector unary operation. */
  7. #define V_Op(vdst,gets,op,vsrc,n) {register int V_i;\
  8.     for(V_i=(n)-1;V_i>=0;V_i--) (vdst)[V_i] gets op ((vsrc)[V_i]);}
  9.  
  10. static void lerp(Knot t, Knot a0, Knot a1, Vect p0, Vect p1, int m, Vect p)
  11. {
  12.     register Knot t0=(a1-t)/(a1-a0), t1=1-t0;
  13.     register int i;
  14.     for (i=m-1; i>=0; i--) p[i] = t0*p0[i] + t1*p1[i];
  15. }
  16.  
  17. /* DialASpline(t,a,p,m,n,work,Cn,interp,val) computes a point val at parameter
  18.     t on a spline with knot values a and control points p. The curve will have
  19.     Cn continuity, and if interp is TRUE it will interpolate the control points.
  20.     Possibilities include Langrange interpolants, Bezier curves, Catmull-Rom
  21.     interpolating splines, and B-spline curves. Points have m coordinates, and
  22.     n+1 of them are provided. The work array must have room for n+1 points.
  23.  */
  24. int DialASpline(Knot t, Knot a[], Vect p[], int m, int n, Vect work[],
  25.                     unsigned int Cn, Bool interp, Vect val)
  26. {
  27.     register int i, j, k, h, lo, hi;
  28.  
  29.     if (Cn>n-1) Cn = n-1;       /* Anything greater gives one polynomial */
  30.     for (k=0; t> a[k]; k++);    /* Find enclosing knot interval */
  31.     for (h=k; t==a[k]; k++);    /* May want to use fewer legs */
  32.     if (k>n) {k = n; if (h>k) h = k;}
  33.     h = 1+Cn - (k-h); k--;
  34.     lo = k-Cn; hi = k+1+Cn;
  35.  
  36.     if (interp) {               /* Lagrange interpolation steps */
  37.         int drop=0;
  38.         if (lo<0) {lo = 0; drop += Cn-k;
  39.                    if (hi-lo<Cn) {drop += Cn-hi; hi = Cn;}}
  40.         if (hi>n) {hi = n; drop += k+1+Cn-n;
  41.                    if (hi-lo<Cn) {drop += lo-(n-Cn); lo = n-Cn;}}
  42.         for (i=lo; i<=hi; i++) V_Op(work[i],=,,p[i],m);
  43.         for (j=1; j<=Cn; j++) {
  44.             for (i=lo; i<=hi-j; i++) {
  45.                 lerp(t,a[i],a[i+j],work[i],work[i+1],m,work[i]);
  46.             }
  47.         }
  48.         h = 1+Cn-drop;
  49.     } else {                    /* Prepare for B-spline steps */
  50.         if (lo<0) {h += lo; lo = 0;}
  51.         for (i=lo; i<=lo+h; i++) V_Op(work[i],=,,p[i],m);
  52.         if (h<0) h = 0;
  53.     }
  54.     for (j=0; j<h; j++) {
  55.         int tmp = 1+Cn-j;
  56.         for (i=h-1; i>=j; i--) {
  57.             lerp(t,a[lo+i],a[lo+i+tmp],work[lo+i],work[lo+i+1],m,work[lo+i+1]);
  58.         }
  59.     }
  60.     V_Op(val,=,,work[lo+h],m);
  61.     return (k);
  62. }
  63.