home *** CD-ROM | disk | FTP | other *** search
/ RISC DISC 3 / RISC_DISC_3.iso / resources / etexts / gems / gemsv / ch3_3 / tricubic.c
Encoding:
Text File  |  1995-03-04  |  2.2 KB  |  87 lines
  1. typedef struct
  2. {
  3.   float           x, y, z;
  4. }               Point;
  5.  
  6. /*
  7.  * TriCubic - tri-cubic interpolation at point, p.
  8.  *   inputs:
  9.  *     p - the interpolation point.
  10.  *     volume - a pointer to the float volume data, stored in x,
  11.  *              y, then z order (x index increasing fastest).
  12.  *     xDim, yDim, zDim - dimensions of the array of volume data.
  13.  *   returns:
  14.  *     the interpolated value at p.
  15.  *   note:
  16.  *     NO range checking is done in this function.
  17.  */
  18.  
  19.  
  20. float           TriCubic (Point p, float *volume, int xDim, int yDim, int zDim)
  21. {
  22.   int             x, y, z;
  23.   register int    i, j, k;
  24.   float           dx, dy, dz;
  25.   register float *pv;
  26.   float           u[4], v[4], w[4];
  27.   float           r[4], q[4];
  28.   float           vox = 0;
  29.   int             xyDim;
  30.  
  31.   xyDim = xDim * yDim;
  32.  
  33.   x = (int) p.x, y = (int) p.y, z = (int) p.z;
  34.   if (x < 0 || x >= xDim || y < 0 || y >= yDim || z < 0 || z >= zDim)
  35.     return (0);
  36.  
  37.   dx = p.x - (float) x, dy = p.y - (float) y, dz = p.z - (float) z;
  38.   pv = volume + (x - 1) + (y - 1) * xDim + (z - 1) * xyDim;
  39.  
  40. # define CUBE(x)   ((x) * (x) * (x))
  41. # define SQR(x)    ((x) * (x))
  42. /*
  43.  #define DOUBLE(x) ((x) + (x))
  44.  #define HALF(x)   ...
  45.  *
  46.  * may also be used to reduce the number of floating point
  47.  * multiplications. The IEEE standard allows for DOUBLE/HALF
  48.  * operations.
  49.  */
  50.  
  51.   /* factors for Catmull-Rom interpolation */
  52.  
  53.   u[0] = -0.5 * CUBE (dx) + SQR (dx) - 0.5 * dx;
  54.   u[1] = 1.5 * CUBE (dx) - 2.5 * SQR (dx) + 1;
  55.   u[2] = -1.5 * CUBE (dx) + 2 * SQR (dx) + 0.5 * dx;
  56.   u[3] = 0.5 * CUBE (dx) - 0.5 * SQR (dx);
  57.  
  58.   v[0] = -0.5 * CUBE (dy) + SQR (dy) - 0.5 * dy;
  59.   v[1] = 1.5 * CUBE (dy) - 2.5 * SQR (dy) + 1;
  60.   v[2] = -1.5 * CUBE (dy) + 2 * SQR (dy) + 0.5 * dy;
  61.   v[3] = 0.5 * CUBE (dy) - 0.5 * SQR (dy);
  62.  
  63.   w[0] = -0.5 * CUBE (dz) + SQR (dz) - 0.5 * dz;
  64.   w[1] = 1.5 * CUBE (dz) - 2.5 * SQR (dz) + 1;
  65.   w[2] = -1.5 * CUBE (dz) + 2 * SQR (dz) + 0.5 * dz;
  66.   w[3] = 0.5 * CUBE (dz) - 0.5 * SQR (dz);
  67.  
  68.   for (k = 0; k < 4; k++)
  69.   {
  70.     q[k] = 0;
  71.     for (j = 0; j < 4; j++)
  72.     {
  73.       r[j] = 0;
  74.       for (i = 0; i < 4; i++)
  75.       {
  76.         r[j] += u[i] * *pv;
  77.         pv++;
  78.       }
  79.       q[k] += v[j] * r[j];
  80.       pv += xDim - 4;
  81.     }
  82.     vox += w[k] * q[k];
  83.     pv += xyDim - 4 * xDim;
  84.   }
  85.   return (vox < 0 ? 0.0 : vox);
  86. }
  87.