home *** CD-ROM | disk | FTP | other *** search
/ The Fred Fish Collection 1.5 / ffcollection-1-5-1992-11.iso / ff_progs / educ / splines.lzh / SPLINES / INTERPL.C < prev    next >
C/C++ Source or Header  |  1991-08-16  |  8KB  |  284 lines
  1. /* The routines in this file are copyright (c) 1987 by Helene (Lee) Taran.
  2.  * Permission is granted for use and free distribution as long as the
  3.  * original author's name is included with the code.
  4.  */
  5.  
  6. #include "spline.h"
  7.  
  8. #define ALTSIGN(n,x) ((n) % 2 ? (x) : -(x))
  9. #define VECTOR(n) ((REAL_POINT *)calloc((n),sizeof(REAL_POINT)))
  10. #define MATRIX(n) ((double *)calloc((n)*(n),sizeof(double)))
  11.  
  12. double G[MAXG];       /* vector of g values */
  13.  
  14. void *calloc();
  15.  
  16. /* Init_GValues should be called once before drawing any curves.  Computes
  17.  * the first MAXG number of Gn values using iteration to save time and 
  18.  * stack space.
  19.  */
  20. Init_GValues()
  21. {
  22.   int n;
  23.   G[0] = 0.0;
  24.   G[1] = 1.0;
  25.   for (n = 2; n < MAXG; n++) 
  26.     G[n] = (4 * G[n-1]  - G[n-2]);
  27. }
  28.  
  29. Print_GValues() /* only for debugging */
  30. {
  31.   int i;
  32.   for (i = 0; i < MAXG; i++)
  33.     printf("G[%d] = %lf\n",i,G[i]);
  34. }
  35.  
  36.  
  37. /* Init_OpenSpline_Matrices : allocates an initializes the L and D matrices
  38.  * that are used to create open splines.
  39.  */
  40. Init_OpenSpline_Matrices(n,L,D)
  41. int n;       /* size of each of the square matrices */
  42. double *L;    /* lower triangular matrix used in LDL^t decomposition */
  43. double *D;    /* diagonal matrix used in LDL^t decomposition */
  44. {
  45.    int row,col;
  46.    for (row = 0; row < n; row++)
  47.        for (col = 0; col < n; col++) {
  48.  
  49.       if (row == col) {  /* assign values to the diagonal elements */
  50.          AssignValue(n,L,row,col,1.0);
  51.          AssignValue(n,D,row,col, G[row+2]/G[row+1]);
  52.       } 
  53.      else if (row < col) {  /* assign values to the upper triangular region */
  54.          AssignValue(n,L,row,col,0.0);
  55.          AssignValue(n,D,row,col,0.0);
  56.      }
  57.      else { /* assign values to the lower triangular region */
  58.          AssignValue(n,D,row,col,0.0);
  59.          if (col < row-1) AssignValue(n,L,row,col,0.0);
  60.          else AssignValue(n,L,row,col,G[row]/G[row+1]);
  61.      }
  62.   }
  63. }
  64.   
  65.  
  66. AssignValue(n,matrix,row,col, value)
  67. int n;                 /* dimension of NxN matrix */                  
  68. double *matrix;         /* <n> by <n> matrix that will obtain the value */
  69. int row,col;           /* where */
  70. double value;           /* value being assigned */
  71. {
  72.    /* each row is sequence of consequtive columns; first loc = 0,0 */
  73.    *(matrix + n * row + col) = value;
  74. }
  75.  
  76. double Value(n,matrix,row,col)
  77. int n;                  /* dimension of NxN matrix */                  
  78. double *matrix;         /* <n> by <n> matrix that will obtain the value */
  79. int row,col;            /* where */
  80. {
  81.    /* each row is a sequence of consequtive columns;  first loc = 0,0 */
  82.   return(*(matrix + n * row + col));
  83. }
  84.  
  85.  
  86. double TValue(n,matrix,row,col)  /* return the transpose value */
  87. int n;                   /* dimension of NxN matrix */                  
  88. double *matrix;          /* <n> by <n> matrix that will obtain the value */
  89. int row,col;             /* where */
  90. {
  91.   return(*(matrix + n * col + row));
  92. }
  93.           
  94.  
  95. Init_ClosedSpline_Matrices(n,L,D)
  96. int n;       /* size of each of the square matrices */
  97. double *L;    /* lower triangular matrix used in LDL^t decomposition */
  98. double *D;    /* diagonal matrix used in LDL^t decomposition */
  99. {
  100.    int row,col;
  101.  
  102.    for (row = 0; row < n; row++)
  103.      for (col =0; col < n; col++)  {
  104.  
  105.      if (row == col) { /* assign diagonal values */
  106.         AssignValue(n,L,row,col,1.0);
  107.         if (row == n-1) 
  108.           AssignValue(n,D,row,col, (G[row+2] - G[row] - ALTSIGN(row,1) * 2)/G[row+1]);
  109.         else AssignValue(n,D,row,col,G[row+2]/G[row+1]);
  110.      }
  111.      else if (row < col) { /* assign to upper triangular region */
  112.          AssignValue(n,L,row,col,0.0);
  113.          AssignValue(n,D,row,col,0.0);
  114.      }
  115.      else {  /* assign to the lower triangular region */
  116.          AssignValue(n,D,row,col,0.0);
  117.          if (row < n-1) 
  118.         if (col == row-1) AssignValue(n,L,row,col,G[row]/G[row+1]);
  119.             else AssignValue(n,L,row,col,0.0);
  120.          else if (col < n-2) AssignValue(n,L,row,col,ALTSIGN(col,-1)/G[col+2]);    
  121.          else AssignValue(n,L,row,col,(G[row]+ALTSIGN(col,-1))/G[row+1]);
  122.       }
  123.    }
  124. }
  125.  
  126.  
  127. PrintMatrix(n,matrix) /* for debugging only */
  128. int n;
  129. double *matrix;
  130. {
  131.   int row, col;
  132.   for (row = 0; row < n; row++) {
  133.     for (col = 0; col < n; col++)
  134.       printf("%lf ",Value(n,matrix,row,col));
  135.     printf("\n");
  136.   }
  137. }
  138.  
  139.  
  140. /* solves for x in the equation LDL^t * x = b and puts the result in <b> */
  141. Solve(n,L,D,b)  
  142. int n;            /* number of elements in <x> and <b> */
  143. double *L;        /* lower triangular matrix : assumed to be initialized */
  144. double *D;        /* diagonal matrix assumed to be initialized */
  145. REAL_POINT b[];   /* vector of known point values & result destination */
  146. {
  147.    int row, col;
  148.  
  149.    for (row = 1; row < n; row++) /* forward substitution */
  150.        for (col = 0; col < row; col++) {
  151.            b[row].x -= (Value(n,L,row,col) * b[col].x);
  152.            b[row].y -= (Value(n,L,row,col) * b[col].y);
  153.        }
  154.  
  155.    b[n-1].x /= Value(n,D,n-1,n-1);
  156.    b[n-1].y /= Value(n,D,n-1,n-1);
  157.  
  158.    for (row = n-2; row >= 0; row--) { /* back substitution */
  159.        b[row].x /= Value(n,D,row,row);
  160.        b[row].y  /= Value(n,D,row,row);
  161.        for (col = row +1; col < n; col++ ) {
  162.            b[row].x -= (TValue(n,L,row,col) * b[col].x);
  163.            b[row].y -= (TValue(n,L,row,col) * b[col].y);
  164.        }
  165.    }
  166. }
  167.  
  168.  
  169. PrintPoint(p)
  170. REAL_POINT *p;
  171. {
  172.   printf("%lf, %lf\n", p->x, p->y);
  173. }
  174.  
  175. PrintVector(n,v)
  176. int n;
  177. REAL_POINT v[];
  178. {
  179.    int i;
  180.    for (i = 0; i < n; i++) PrintPoint(&v[i]);
  181.    printf("\n");
  182. }
  183.  
  184.  
  185. FillVector(n,v,list,xadj,yadj)
  186. int n;                /* size of vectors <= length of list */
  187. REAL_POINT v[];       /* vector to receive the x,y values in list */
  188. DLISTPTR list;        /* first element to get the values from  */ 
  189. float xadj, yadj;     /* how to adjust the x,y values  (1 if no adjustment) */
  190. {
  191.    DLISTPTR tmp = list; 
  192.    int i;
  193.     
  194.    if (n > 0) {
  195.       for (i = 0; i < n; i++) {
  196.           v[i].x = xadj * (POINT(tmp)->x) ;
  197.           v[i].y = yadj * (POINT(tmp)->y) ;
  198.           tmp = NEXT(tmp);
  199.      }
  200.   }
  201. }
  202.  
  203.  
  204. /* ListVector : takes a vector of points and makes it into a list 
  205.  * can be passed to the any of the bspline drawing routines.
  206.  */
  207. DLISTPTR ListVector(n,v)
  208. int n;                     /* n <= length of <v>  */
  209. REAL_POINT v[];
  210. {
  211.    DLISTPTR list = (DLISTPTR)calloc(1,sizeof(DLIST_ELEMENT));
  212.    DLISTPTR element;
  213.    int i;
  214.  
  215.    Init_List(list);
  216.    for (i = 0; i < n; i++) {
  217.        element = (DLISTPTR)calloc(1,sizeof(DLIST_ELEMENT));
  218.        element->contents = &(v[i]);
  219.        INSERT_LAST(element,list);  /* insert in same order */
  220.    }
  221.    return(list);
  222. }
  223.  
  224.  
  225. Draw_Open_Ispline(w,control_points)
  226. struct Window *w;
  227. DLISTPTR control_points;
  228. {
  229.    int n = LENGTH(control_points) -2;
  230.    DLISTPTR newpoints;
  231.  
  232.    REAL_POINT *b = VECTOR(n);
  233.    double *L  = MATRIX(n);
  234.    double *D  = MATRIX(n);
  235.  
  236.    FillVector(n, b, NEXT(FIRST(control_points)), 6.0, 6.0);
  237.    b[0].x -= POINT(FIRST(control_points))->x ;
  238.    b[0].y -= POINT(FIRST(control_points))->y ;
  239.    b[n-1].x -= POINT(LAST(control_points))->x ;
  240.    b[n-1].y -= POINT(LAST(control_points))->y ;
  241.  
  242.    Init_OpenSpline_Matrices(n,L,D);
  243.    Solve(n,L,D,b);
  244.  
  245.    newpoints =  ListVector(n,b);
  246.    MOVE_FIRST(control_points,newpoints);
  247.    MOVE_LAST(control_points, newpoints);
  248.  
  249.    Draw_Natural_Bspline(w,newpoints);
  250.  
  251.    MOVE_FIRST(newpoints,control_points);
  252.    MOVE_LAST(newpoints,control_points);
  253.    ClearList(newpoints,FALSE);      
  254.    free(newpoints);  free(L); free(D); free(b);
  255.    
  256. }
  257.  
  258.  
  259. Draw_Closed_Ispline(w,control_points)
  260. struct Window *w;
  261. DLISTPTR control_points;
  262. {
  263.    int n = LENGTH(control_points);
  264.    DLISTPTR newpoints;
  265.  
  266.    REAL_POINT *b = VECTOR(n);
  267.    double *L  = MATRIX(n);
  268.    double *D  = MATRIX(n);
  269.  
  270.    FillVector(n, b, FIRST(control_points), 6.0, 6.0);
  271.  
  272.    Init_ClosedSpline_Matrices(n,L,D);
  273.    Solve(n,L,D,b);
  274.  
  275.    newpoints =  ListVector(n,b);
  276.    Draw_Closed_Bspline(w,newpoints);
  277.  
  278.    ClearList(newpoints,FALSE);      
  279.    free(newpoints);  free(L); free(D); free(b);
  280.    
  281. }
  282.  
  283.  
  284.