home *** CD-ROM | disk | FTP | other *** search
/ Fresh Fish 7 / FreshFishVol7.bin / new / misc / sci / splines / interpl.c < prev    next >
C/C++ Source or Header  |  1994-09-16  |  8KB  |  335 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. /*
  7.  * Reformatted, and modified to compile without warnings and errors
  8.  * under SAS/C -6.5x by gduncan@philips.oz.au (GMD). This included
  9.  * proto generation, renaming of some literals to avoid name collisions,
  10.  * and Amiga Version string (also displayed in Title bar).
  11.  * No original version number, this version arbitrarily named 1.1. 
  12.  * - otherwise no functional changes.
  13.  * GMD - Sep 94
  14.  */
  15.  
  16. #include "all.h"
  17.  
  18. #define ALTSIGN(n,x) ((n) % 2 ? (x) : -(x))
  19. #define VECTOR(n) ((REAL_POINT *)calloc((n),sizeof(REAL_POINT)))
  20. #define MATRIX(n) ((double *)calloc((n)*(n),sizeof(double)))
  21.  
  22. double G[MAXG];            /* vector of g values */
  23.  
  24.  
  25. /* Init_GValues should be called once before drawing any curves.  Computes
  26.  * the first MAXG number of Gn values using iteration to save time and 
  27.  * stack space.
  28.  */
  29.  
  30. void
  31. Init_GValues ()
  32. {
  33.   int n;
  34.   G[0] = 0.0;
  35.   G[1] = 1.0;
  36.   for (n = 2; n < MAXG; n++)
  37.     G[n] = (4 * G[n - 1] - G[n - 2]);
  38. }
  39.  
  40. void
  41. Print_GValues ()        /* only for debugging */
  42. {
  43.   int i;
  44.   for (i = 0; i < MAXG; i++)
  45.     printf ("G[%d] = %lf\n", i, G[i]);
  46. }
  47.  
  48.  
  49. /* Init_OpenSpline_Matrices : allocates an initializes the L and D matrices
  50.  * that are used to create open splines.
  51.  */
  52.  
  53. void
  54. Init_OpenSpline_Matrices (n, L, D)
  55.      int n;            /* size of each of the square matrices */
  56.      double *L;            /* lower triangular matrix used in LDL^t decomposition */
  57.      double *D;            /* diagonal matrix used in LDL^t decomposition */
  58. {
  59.   int row, col;
  60.   for (row = 0; row < n; row++)
  61.     for (col = 0; col < n; col++)
  62.       {
  63.  
  64.     if (row == col)
  65.       {            /* assign values to the diagonal elements */
  66.         AssignValue (n, L, row, col, 1.0);
  67.         AssignValue (n, D, row, col, G[row + 2] / G[row + 1]);
  68.       }
  69.     else if (row < col)
  70.       {            /* assign values to the upper triangular region */
  71.         AssignValue (n, L, row, col, 0.0);
  72.         AssignValue (n, D, row, col, 0.0);
  73.       }
  74.     else
  75.       {            /* assign values to the lower triangular region */
  76.         AssignValue (n, D, row, col, 0.0);
  77.         if (col < row - 1)
  78.           AssignValue (n, L, row, col, 0.0);
  79.         else
  80.           AssignValue (n, L, row, col, G[row] / G[row + 1]);
  81.       }
  82.       }
  83. }
  84.  
  85. void
  86. AssignValue (n, matrix, row, col, 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.      double value;        /* value being assigned */
  91. {
  92.   /* each row is sequence of consequtive columns; first loc = 0,0 */
  93.   *(matrix + n * row + col) = value;
  94. }
  95.  
  96. double
  97. Value (n, matrix, row, col)
  98.      int n;            /* dimension of NxN matrix */
  99.      double *matrix;        /* <n> by <n> matrix that will obtain the value */
  100.      int row, col;        /* where */
  101. {
  102.   /* each row is a sequence of consequtive columns;  first loc = 0,0 */
  103.   return (*(matrix + n * row + col));
  104. }
  105.  
  106.  
  107. double
  108. TValue (n, matrix, row, col)    /* return the transpose value */
  109.      int n;            /* dimension of NxN matrix */
  110.      double *matrix;        /* <n> by <n> matrix that will obtain the value */
  111.      int row, col;        /* where */
  112. {
  113.   return (*(matrix + n * col + row));
  114. }
  115.  
  116. void
  117. Init_ClosedSpline_Matrices (n, L, D)
  118.      int n;            /* size of each of the square matrices */
  119.      double *L;            /* lower triangular matrix used in LDL^t decomposition */
  120.      double *D;            /* diagonal matrix used in LDL^t decomposition */
  121. {
  122.   int row, col;
  123.  
  124.   for (row = 0; row < n; row++)
  125.     for (col = 0; col < n; col++)
  126.       {
  127.  
  128.     if (row == col)
  129.       {            /* assign diagonal values */
  130.         AssignValue (n, L, row, col, 1.0);
  131.         if (row == n - 1)
  132.           AssignValue (n, D, row, col, (G[row + 2] - G[row] - ALTSIGN (row, 1) * 2) / G[row + 1]);
  133.         else
  134.           AssignValue (n, D, row, col, G[row + 2] / G[row + 1]);
  135.       }
  136.     else if (row < col)
  137.       {            /* assign to upper triangular region */
  138.         AssignValue (n, L, row, col, 0.0);
  139.         AssignValue (n, D, row, col, 0.0);
  140.       }
  141.     else
  142.       {            /* assign to the lower triangular region */
  143.         AssignValue (n, D, row, col, 0.0);
  144.         if (row < n - 1)
  145.           if (col == row - 1)
  146.         AssignValue (n, L, row, col, G[row] / G[row + 1]);
  147.           else
  148.         AssignValue (n, L, row, col, 0.0);
  149.         else if (col < n - 2)
  150.           AssignValue (n, L, row, col, ALTSIGN (col, -1) / G[col + 2]);
  151.         else
  152.           AssignValue (n, L, row, col, (G[row] + ALTSIGN (col, -1)) / G[row + 1]);
  153.       }
  154.       }
  155. }
  156.  
  157. void
  158. PrintMatrix (n, matrix)        /* for debugging only */
  159.      int n;
  160.      double *matrix;
  161. {
  162.   int row, col;
  163.   for (row = 0; row < n; row++)
  164.     {
  165.       for (col = 0; col < n; col++)
  166.     printf ("%lf ", Value (n, matrix, row, col));
  167.       printf ("\n");
  168.     }
  169. }
  170.  
  171.  
  172. /* solves for x in the equation LDL^t * x = b and puts the result in <b> */
  173.  
  174. void
  175. Solve (n, L, D, b)
  176.      int n;            /* number of elements in <x> and <b> */
  177.      double *L;            /* lower triangular matrix : assumed to be initialized */
  178.      double *D;            /* diagonal matrix assumed to be initialized */
  179.      REAL_POINT b[];        /* vector of known point values & result destination */
  180. {
  181.   int row, col;
  182.  
  183.   for (row = 1; row < n; row++)    /* forward substitution */
  184.     for (col = 0; col < row; col++)
  185.       {
  186.     b[row].x -= (Value (n, L, row, col) * b[col].x);
  187.     b[row].y -= (Value (n, L, row, col) * b[col].y);
  188.       }
  189.  
  190.   b[n - 1].x /= Value (n, D, n - 1, n - 1);
  191.   b[n - 1].y /= Value (n, D, n - 1, n - 1);
  192.  
  193.   for (row = n - 2; row >= 0; row--)
  194.     {                /* back substitution */
  195.       b[row].x /= Value (n, D, row, row);
  196.       b[row].y /= Value (n, D, row, row);
  197.       for (col = row + 1; col < n; col++)
  198.     {
  199.       b[row].x -= (TValue (n, L, row, col) * b[col].x);
  200.       b[row].y -= (TValue (n, L, row, col) * b[col].y);
  201.     }
  202.     }
  203. }
  204.  
  205. void
  206. PrintPoint (p)
  207.      REAL_POINT *p;
  208. {
  209.   printf ("%lf, %lf\n", p->x, p->y);
  210. }
  211.  
  212. void
  213. PrintVector (n, v)
  214.      int n;
  215.      REAL_POINT v[];
  216. {
  217.   int i;
  218.   for (i = 0; i < n; i++)
  219.     PrintPoint (&v[i]);
  220.   printf ("\n");
  221. }
  222.  
  223. void
  224. FillVector (n, v, list, xadj, yadj)
  225.      int n;            /* size of vectors <= length of list */
  226.      REAL_POINT v[];        /* vector to receive the x,y values in list */
  227.      DLISTPTR list;        /* first element to get the values from  */
  228.      float xadj, yadj;        /* how to adjust the x,y values  (1 if no adjustment) */
  229. {
  230.   DLISTPTR tmp = list;
  231.   int i;
  232.  
  233.   if (n > 0)
  234.     {
  235.       for (i = 0; i < n; i++)
  236.     {
  237.       v[i].x = xadj * (POINT (tmp)->x);
  238.       v[i].y = yadj * (POINT (tmp)->y);
  239.       tmp = NEXT (tmp);
  240.     }
  241.     }
  242. }
  243.  
  244.  
  245. /* ListVector : takes a vector of points and makes it into a list 
  246.  * can be passed to the any of the bspline drawing routines.
  247.  */
  248.  
  249.  
  250. DLISTPTR
  251. ListVector (n, v)
  252.      int n;            /* n <= length of <v>  */
  253.      REAL_POINT v[];
  254. {
  255.   DLISTPTR list = (DLISTPTR) calloc (1, sizeof (DLIST_ELEMENT));
  256.   DLISTPTR element;
  257.   int i;
  258.  
  259.   Init_List (list);
  260.   for (i = 0; i < n; i++)
  261.     {
  262.       element = (DLISTPTR) calloc (1, sizeof (DLIST_ELEMENT));
  263.       element->contents = &(v[i]);
  264.       INSERT_LAST (element, list);    /* insert in same order */
  265.     }
  266.   return (list);
  267. }
  268.  
  269.  
  270. void
  271. Draw_Open_Ispline (w, control_points)
  272.      WINDOW *w;
  273.      DLISTPTR control_points;
  274. {
  275.   int n = LENGTH (control_points) - 2;
  276.   DLISTPTR newpoints;
  277.  
  278.   REAL_POINT *b = VECTOR (n);
  279.   double *L = MATRIX (n);
  280.   double *D = MATRIX (n);
  281.  
  282.   FillVector (n, b, NEXT (FIRST (control_points)), 6.0, 6.0);
  283.   b[0].x -= POINT (FIRST (control_points))->x;
  284.   b[0].y -= POINT (FIRST (control_points))->y;
  285.   b[n - 1].x -= POINT (LAST (control_points))->x;
  286.   b[n - 1].y -= POINT (LAST (control_points))->y;
  287.  
  288.   Init_OpenSpline_Matrices (n, L, D);
  289.   Solve (n, L, D, b);
  290.  
  291.   newpoints = ListVector (n, b);
  292.   MOVE_FIRST (control_points, newpoints);
  293.   MOVE_LAST (control_points, newpoints);
  294.  
  295.   Draw_Natural_Bspline (w, newpoints);
  296.  
  297.   MOVE_FIRST (newpoints, control_points);
  298.   MOVE_LAST (newpoints, control_points);
  299.   ClearList (newpoints, FALSE);
  300.   free (newpoints);
  301.   free (L);
  302.   free (D);
  303.   free (b);
  304.  
  305. }
  306.  
  307. void
  308. Draw_Closed_Ispline (w, control_points)
  309.      WINDOW *w;
  310.      DLISTPTR control_points;
  311. {
  312.   int n = LENGTH (control_points);
  313.   DLISTPTR newpoints;
  314.  
  315.   REAL_POINT *b = VECTOR (n);
  316.   double *L = MATRIX (n);
  317.   double *D = MATRIX (n);
  318.  
  319.   FillVector (n, b, FIRST (control_points), 6.0, 6.0);
  320.  
  321.   Init_Clos