home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / exampleCode / inventor / inventorTemplates1.1.2 / SplineClass.c++ < prev    next >
Encoding:
C/C++ Source or Header  |  1994-08-02  |  3.8 KB  |  184 lines
  1. /*
  2.  * Copyright (C) 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. #include <stdio.h>
  18. #include <stdlib.h>
  19.  
  20. #include "SplineClass.h"
  21.  
  22. //
  23. // SPLINE STUFF!!!!!
  24. // I stole these mostly from Drew Olbrich
  25. //
  26.  
  27. SplineClass::SplineClass(int num, double *xvals, double *fvals)
  28. {
  29.   int i;
  30.   double t;
  31.   double *h;
  32.   double *u;
  33.   double *ta;
  34.   double *td;
  35.   double *tb;
  36.   double *tf;
  37.   screwed = 0;
  38.  
  39.   if (num < 3) {
  40.       fprintf(stderr, "too few points\n");
  41.       screwed = 1;
  42.       return;
  43.   }
  44.  
  45.   x = (double *) malloc((num*sizeof(double)));
  46.   if (x == NULL) 
  47.   {
  48.     fprintf(stderr, "Could not allocate spline data.\n");
  49.     screwed = 1;
  50.     return;
  51.   }
  52.   for (i = 0; i < num; i++)
  53.     x[i] = xvals[i];
  54.  
  55.   n = num;
  56.   min_x = xvals[0];
  57.   max_x = xvals[n - 1];
  58.  
  59.   a = (double *) malloc(n*sizeof(double));
  60.   b = (double *) malloc(n*sizeof(double));
  61.   c = (double *) malloc(n*sizeof(double));
  62.   d = (double *) malloc(n*sizeof(double));
  63.  
  64.   if (d == NULL)
  65.   {
  66.     fprintf(stderr, "Could not allocate spline coefficient data.\n");
  67.     screwed = 1;
  68.     return;
  69.   }
  70.  
  71.   h = (double *) malloc((n+1)*sizeof(double));
  72.   u = (double *) malloc((n+1)*sizeof(double));
  73.  
  74.   ta = (double *) malloc((n+1)*sizeof(double));
  75.   td = (double *) malloc((n+1)*sizeof(double));
  76.   tb = (double *) malloc((n+1)*sizeof(double));
  77.   tf = (double *) malloc((n+1)*sizeof(double));
  78.  
  79.   if (tf == NULL)
  80.   {
  81.     fprintf(stderr, "Could not allocate temporary spline data.\n");
  82.     screwed = 1;
  83.     return;
  84.   }
  85.  
  86.   for (i = 0; i < n - 1; i++) {
  87.     h[i] = x[i + 1] - x[i];
  88.     u[i] = (fvals[i + 1] - fvals[i])/h[i];
  89.     a[i] = fvals[i];
  90.   }
  91.  
  92.   for (i = 0; i < (n - 2); i++)
  93.     td[i] = 2.0*(h[i] + h[i + 1]);
  94.  
  95.   for (i = 0; i < (n - 2); i++)
  96.     ta[i] = tb[i] = h[i + 1];
  97.  
  98.   for (i = 0; i < (n - 2); i++)
  99.     tf[i] = 3.0*(u[i + 1] - u[i]);
  100.  
  101.   for (i = 1; i < (n - 2); i++)
  102.   {
  103.     t = ta[i - 1]/td[i - 1];
  104.     td[i] = td[i] - t*tb[i - 1];
  105.     tf[i] = tf[i] - t*tf[i - 1];
  106.   }
  107.  
  108.   c[0] = 0.0;
  109.   c[n - 1] = 0.0;
  110.   c[n - 2] = tf[n - 3]/td[n - 3];
  111.   for (i = (n - 4); i >= 0; i--)
  112.     c[i + 1] = (tf[i] - tb[i]*c[i + 2])/td[i];
  113.  
  114.   for (i = 0; i < (n - 1); i++) {
  115.     b[i] = u[i] - h[i]*(2.0*c[i] + c[i + 1])/3.0;
  116.     d[i] = (c[i + 1] - c[i])/(h[i]*3.0);
  117.   }
  118.  
  119.   free(h);
  120.   free(u);
  121.  
  122.   free(ta);
  123.   free(td);
  124.   free(tb);
  125.   free(tf);
  126.  
  127.   screwed = 0;
  128. }
  129.  
  130. short SplineClass::isOk()
  131. {
  132.     return screwed;
  133. }
  134.  
  135. double SplineClass::evaluate(double t)
  136. {
  137.   int i;
  138.   int high, low, mid;
  139.   double result, qq;
  140.  
  141.   if (screwed) {
  142.       fprintf(stderr, "screwed!\n");
  143.       return 0.0;
  144.   }
  145.  
  146.   // I am changing this behavior so before min_x you get min_x and
  147.   // after max_x you get max_x
  148.   //  if ((t < min_x) || (t > max_x))
  149.   //    return 0.0;
  150.  
  151.   low = 0;
  152.   high = (int) n - 1;
  153.  
  154.   if (t < min_x) return a[low];
  155.   if (t > max_x) return a[high];
  156.  
  157.   i = -1;
  158.   while (i == -1 && low <= high) {
  159.     mid = (low + high) >> 1;
  160.     if (t < x[mid])
  161.       high = mid - 1;
  162.     else if (mid + 1 < n && t > x[mid + 1])
  163.       low = mid + 1;
  164.     else
  165.       i = mid;
  166.   }
  167.   qq = t - x[i];
  168.  
  169.   result = a[i]    
  170.     + qq*(b[i] + qq*(c[i] + qq*d[i]));
  171.   
  172.   return result;
  173. }
  174.  
  175.  
  176. SplineClass::~SplineClass()
  177. {
  178.     free(x);
  179.     free(a);
  180.     free(b);
  181.     free(c);
  182.     free(d);
  183. }
  184.