contrib/mul/vil3d/vil3d_tricub_interp.txx Source File

Go to the documentation of this file.
00001 #ifndef vil3d_tricub_interp_txx_
00002 #define vil3d_tricub_interp_txx_
00003 //:
00004 // \file
00005 // \brief Tricubic interpolation functions for 3D images
00006 
00007 #include "vil3d_tricub_interp.h"
00008 #include <vcl_limits.h>
00009 #include <vil/vil_round.h>
00010 
00011 
00012 namespace
00013 {
00014   void get_tricubic_coeff( double t,
00015                            double &c0, double &c1, double &c2, double &c3 )
00016   {
00017     c0 = ((2-t)*t-1)*t;    // -1
00018     c1 = (3*t-5)*t*t+2;    //  0
00019     c2 = ((4-3*t)*t+1)*t;  // +1
00020     c3 = (t-1)*t*t;        // +2
00021   }
00022 }
00023 
00024 
00025 template<class T>
00026 // double vil3d_tricub_interp_unsafe(double x, double y, double z,
00027 double vil3d_tricub_interp_raw(double x, double y, double z,
00028                                const T* data,
00029                                vcl_ptrdiff_t xstep, vcl_ptrdiff_t ystep, vcl_ptrdiff_t zstep)
00030 {
00031     int p1x=int(x);
00032     double normx = x-p1x;
00033     int p1y=int(y);
00034     double normy = y-p1y;
00035     int p1z=int(z);
00036     double normz = z-p1z;
00037 
00038     const T* pix1 = data + p1y*ystep + p1x*xstep + p1z*zstep;
00039 
00040     // use separability.
00041     // the s's are for the x-direction, the t's for the y-direction,
00042     // the u's are for the z-direction.
00043     double s0, s1, s2, s3;
00044     get_tricubic_coeff( normx, s0, s1, s2, s3 );
00045     double t0, t1, t2, t3;
00046     get_tricubic_coeff( normy, t0, t1, t2, t3 );
00047     double u0, u1, u2, u3;
00048     get_tricubic_coeff( normz, u0, u1, u2, u3 );
00049 
00050 #define vil3d_I(dx,dy,dz) (pix1[(dx)*xstep+(dy)*ystep+(dz)*zstep])
00051 
00052     double xi00 = s0*vil3d_I(-1,-1,-1) + s1*vil3d_I(+0,-1,-1) + s2*vil3d_I(+1,-1,-1) + s3*vil3d_I(+2,-1,-1);
00053     double xi10 = s0*vil3d_I(-1,+0,-1) + s1*vil3d_I(+0,+0,-1) + s2*vil3d_I(+1,+0,-1) + s3*vil3d_I(+2,+0,-1);
00054     double xi20 = s0*vil3d_I(-1,+1,-1) + s1*vil3d_I(+0,+1,-1) + s2*vil3d_I(+1,+1,-1) + s3*vil3d_I(+2,+1,-1);
00055     double xi30 = s0*vil3d_I(-1,+2,-1) + s1*vil3d_I(+0,+2,-1) + s2*vil3d_I(+1,+2,-1) + s3*vil3d_I(+2,+2,-1);
00056 
00057     double xi01 = s0*vil3d_I(-1,-1,+0) + s1*vil3d_I(+0,-1,+0) + s2*vil3d_I(+1,-1,+0) + s3*vil3d_I(+2,-1,+0);
00058     double xi11 = s0*vil3d_I(-1,+0,+0) + s1*vil3d_I(+0,+0,+0) + s2*vil3d_I(+1,+0,+0) + s3*vil3d_I(+2,+0,+0);
00059     double xi21 = s0*vil3d_I(-1,+1,+0) + s1*vil3d_I(+0,+1,+0) + s2*vil3d_I(+1,+1,+0) + s3*vil3d_I(+2,+1,+0);
00060     double xi31 = s0*vil3d_I(-1,+2,+0) + s1*vil3d_I(+0,+2,+0) + s2*vil3d_I(+1,+2,+0) + s3*vil3d_I(+2,+2,+0);
00061 
00062     double xi02 = s0*vil3d_I(-1,-1,+1) + s1*vil3d_I(+0,-1,+1) + s2*vil3d_I(+1,-1,+1) + s3*vil3d_I(+2,-1,+1);
00063     double xi12 = s0*vil3d_I(-1,+0,+1) + s1*vil3d_I(+0,+0,+1) + s2*vil3d_I(+1,+0,+1) + s3*vil3d_I(+2,+0,+1);
00064     double xi22 = s0*vil3d_I(-1,+1,+1) + s1*vil3d_I(+0,+1,+1) + s2*vil3d_I(+1,+1,+1) + s3*vil3d_I(+2,+1,+1);
00065     double xi32 = s0*vil3d_I(-1,+2,+1) + s1*vil3d_I(+0,+2,+1) + s2*vil3d_I(+1,+2,+1) + s3*vil3d_I(+2,+2,+1);
00066 
00067     double xi03 = s0*vil3d_I(-1,-1,+2) + s1*vil3d_I(+0,-1,+2) + s2*vil3d_I(+1,-1,+2) + s3*vil3d_I(+2,-1,+2);
00068     double xi13 = s0*vil3d_I(-1,+0,+2) + s1*vil3d_I(+0,+0,+2) + s2*vil3d_I(+1,+0,+2) + s3*vil3d_I(+2,+0,+2);
00069     double xi23 = s0*vil3d_I(-1,+1,+2) + s1*vil3d_I(+0,+1,+2) + s2*vil3d_I(+1,+1,+2) + s3*vil3d_I(+2,+1,+2);
00070     double xi33 = s0*vil3d_I(-1,+2,+2) + s1*vil3d_I(+0,+2,+2) + s2*vil3d_I(+1,+2,+2) + s3*vil3d_I(+2,+2,+2);
00071 #undef vil3d_I
00072 
00073     double val0 = xi00*t0 + xi10*t1 + xi20*t2 + xi30*t3;
00074     double val1 = xi01*t0 + xi11*t1 + xi21*t2 + xi31*t3;
00075     double val2 = xi02*t0 + xi12*t1 + xi22*t2 + xi32*t3;
00076     double val3 = xi03*t0 + xi13*t1 + xi23*t2 + xi33*t3;
00077 
00078     double val = 0.125 * ( val0*u0 + val1*u1 + val2*u2 + val3*u3 );
00079 
00080     return val;
00081 }
00082 
00083 #pragma optimize( "tpgsy", off )
00084 
00085 template<class T>
00086 double vil3d_tricub_interp_safe_trilinear_extend(double x, double y, double z,
00087                                                  const T* data,
00088                                                  int nx, int ny, int nz,
00089                                                  vcl_ptrdiff_t xstep, vcl_ptrdiff_t ystep, vcl_ptrdiff_t zstep)
00090 {
00091   if( x>=1 && y>=1 && z>=1 &&
00092       x<=nx-3 && y<=ny-3 && z<=nz-3 )
00093     return vil3d_tricub_interp_raw(x,y,z,data,xstep,ystep,zstep);
00094 
00095   int p1x=vil_round_floor(x);
00096   double normx = x - p1x;
00097   int p1y=vil_round_floor(y);
00098   double normy = y - p1y;
00099   int p1z=vil_round_floor(z);
00100   double normz = z - p1z;
00101 
00102   const T* pix1 = data + p1y*ystep + p1x*xstep + p1z*zstep;
00103   double xi00 = vcl_numeric_limits<double>::infinity();
00104   double xi10 = vcl_numeric_limits<double>::infinity();
00105   double xi20 = vcl_numeric_limits<double>::infinity();
00106   double xi30 = vcl_numeric_limits<double>::infinity();
00107   double xi01 = vcl_numeric_limits<double>::infinity();
00108   double xi11 = vcl_numeric_limits<double>::infinity();
00109   double xi21 = vcl_numeric_limits<double>::infinity();
00110   double xi31 = vcl_numeric_limits<double>::infinity();
00111   double xi02 = vcl_numeric_limits<double>::infinity();
00112   double xi12 = vcl_numeric_limits<double>::infinity();
00113   double xi22 = vcl_numeric_limits<double>::infinity();
00114   double xi32 = vcl_numeric_limits<double>::infinity();
00115   double xi03 = vcl_numeric_limits<double>::infinity();
00116   double xi13 = vcl_numeric_limits<double>::infinity();
00117   double xi23 = vcl_numeric_limits<double>::infinity();
00118   double xi33 = vcl_numeric_limits<double>::infinity();
00119 
00120   double val0 = vcl_numeric_limits<double>::infinity();
00121   double val1 = vcl_numeric_limits<double>::infinity();
00122   double val2 = vcl_numeric_limits<double>::infinity();
00123   double val3 = vcl_numeric_limits<double>::infinity();
00124 
00125   double val  = vcl_numeric_limits<double>::infinity();
00126 
00127 #define vil3d_I(dx,dy,dz) (pix1[(dx)*xstep+(dy)*ystep+(dz)*zstep])
00128   if( x < 0 ) x = 0.0;
00129   if( y < 0 ) y = 0.0;
00130   if( z < 0 ) z = 0.0;
00131   if( x > nx-1 ) x = nx-1;
00132   if( y > ny-1 ) y = ny-1;
00133   if( z > nz-1 ) z = nz-1;
00134 
00135   if( z<1 || z> nz-3)
00136   {
00137     if( y<1  || y> ny-3)
00138     {
00139       if( x<1  || x>nx-3)
00140       {
00141         xi11 = vil3d_I(+0,+0,+0) + (vil3d_I(+1,+0,+0)-vil3d_I(+0,+0,+0))*normx;
00142         xi21 = vil3d_I(+0,+1,+0) + (vil3d_I(+1,+1,+0)-vil3d_I(+0,+1,+0))*normx;
00143 
00144         xi12 = vil3d_I(+0,+0,+1) + (vil3d_I(+1,+0,+1)-vil3d_I(+0,+0,+1))*normx;
00145         xi22 = vil3d_I(+0,+1,+1) + (vil3d_I(+1,+1,+1)-vil3d_I(+0,+1,+1))*normx;
00146       }
00147       else
00148       {
00149         double s0, s1, s2, s3;
00150         get_tricubic_coeff( normx, s0, s1, s2, s3 );
00151 
00152         xi11 = 0.5 * ( s0*vil3d_I(-1,+0,+0) +
00153                        s1*vil3d_I(+0,+0,+0) +
00154                        s2*vil3d_I(+1,+0,+0) +
00155                        s3*vil3d_I(+2,+0,+0) );
00156         xi21 = 0.5 * ( s0*vil3d_I(-1,+1,+0) +
00157                        s1*vil3d_I(+0,+1,+0) +
00158                        s2*vil3d_I(+1,+1,+0) +
00159                        s3*vil3d_I(+2,+1,+0) );
00160 
00161         xi12 = 0.5 * ( s0*vil3d_I(-1,+0,+1) +
00162                        s1*vil3d_I(+0,+0,+1) +
00163                        s2*vil3d_I(+1,+0,+1) +
00164                        s3*vil3d_I(+2,+0,+1) );
00165         xi22 = 0.5 * ( s0*vil3d_I(-1,+1,+1) +
00166                        s1*vil3d_I(+0,+1,+1) +
00167                        s2*vil3d_I(+1,+1,+1) +
00168                        s3*vil3d_I(+2,+1,+1) );
00169       } // end if/else x
00170 
00171       val1 = xi11 + (xi21 - xi11)*normy;
00172       val2 = xi12 + (xi22 - xi12)*normy;
00173 
00174     } // end of if y
00175     else
00176     {
00177       if( x<1  || x>nx-3)
00178       {
00179         xi01 = vil3d_I(+0,-1,+0) + (vil3d_I(+1,-1,+0)-vil3d_I(+0,-1,+0))*normx;
00180         xi11 = vil3d_I(+0,+0,+0) + (vil3d_I(+1,+0,+0)-vil3d_I(+0,+0,+0))*normx;
00181         xi21 = vil3d_I(+0,+1,+0) + (vil3d_I(+1,+1,+0)-vil3d_I(+0,+1,+0))*normx;
00182         xi31 = vil3d_I(+0,+2,+0) + (vil3d_I(+1,+2,+0)-vil3d_I(+0,+2,+0))*normx;
00183 
00184         xi02 = vil3d_I(+0,-1,+1) + (vil3d_I(+1,-1,+1)-vil3d_I(+0,-1,+1))*normx;
00185         xi12 = vil3d_I(+0,+0,+1) + (vil3d_I(+1,+0,+1)-vil3d_I(+0,+0,+1))*normx;
00186         xi22 = vil3d_I(+0,+1,+1) + (vil3d_I(+1,+1,+1)-vil3d_I(+0,+1,+1))*normx;
00187         xi32 = vil3d_I(+0,+2,+1) + (vil3d_I(+1,+2,+1)-vil3d_I(+0,+2,+1))*normx;
00188       }
00189       else
00190       {
00191         double s0, s1, s2, s3;
00192         get_tricubic_coeff( normx, s0, s1, s2, s3 );
00193 
00194         xi01 = 0.5 * ( s0*vil3d_I(-1,-1,+0) +
00195                        s1*vil3d_I(+0,-1,+0) +
00196                        s2*vil3d_I(+1,-1,+0) +
00197                        s3*vil3d_I(+2,-1,+0) );
00198         xi11 = 0.5 * ( s0*vil3d_I(-1,+0,+0) +
00199                        s1*vil3d_I(+0,+0,+0) +
00200                        s2*vil3d_I(+1,+0,+0) +
00201                        s3*vil3d_I(+2,+0,+0) );
00202         xi21 = 0.5 * ( s0*vil3d_I(-1,+1,+0) +
00203                        s1*vil3d_I(+0,+1,+0) +
00204                        s2*vil3d_I(+1,+1,+0) +
00205                        s3*vil3d_I(+2,+1,+0) );
00206         xi31 = 0.5 * ( s0*vil3d_I(-1,+2,+0) +
00207                        s1*vil3d_I(+0,+2,+0) +
00208                        s2*vil3d_I(+1,+2,+0) +
00209                        s3*vil3d_I(+2,+2,+0) );
00210 
00211         xi02 = 0.5 * ( s0*vil3d_I(-1,-1,+1) +
00212                        s1*vil3d_I(+0,-1,+1) +
00213                        s2*vil3d_I(+1,-1,+1) +
00214                        s3*vil3d_I(+2,-1,+1) );
00215         xi12 = 0.5 * ( s0*vil3d_I(-1,+0,+1) +
00216                        s1*vil3d_I(+0,+0,+1) +
00217                        s2*vil3d_I(+1,+0,+1) +
00218                        s3*vil3d_I(+2,+0,+1) );
00219         xi22 = 0.5 * ( s0*vil3d_I(-1,+1,+1) +
00220                        s1*vil3d_I(+0,+1,+1) +
00221                        s2*vil3d_I(+1,+1,+1) +
00222                        s3*vil3d_I(+2,+1,+1) );
00223         xi32 = 0.5 * ( s0*vil3d_I(-1,+2,+1) +
00224                        s1*vil3d_I(+0,+2,+1) +
00225                        s2*vil3d_I(+1,+2,+1) +
00226                        s3*vil3d_I(+2,+2,+1) );
00227       } // end if/else x
00228 
00229       double t0, t1, t2, t3;
00230       get_tricubic_coeff( normy, t0, t1, t2, t3 );
00231       val1 = 0.5 * ( xi01*t0 + xi11*t1 + xi21*t2 + xi31*t3 );
00232       val2 = 0.5 * ( xi02*t0 + xi12*t1 + xi22*t2 + xi32*t3 );
00233     } // end of else y
00234 
00235     val = val1 + (val2 - val1) * normz;
00236   } //end of if z
00237   else
00238   {
00239     if( y<1  || y> ny-3)
00240     {
00241       if( x<1  || x>nx-3)
00242       {
00243         xi10 = vil3d_I(+0,+0,-1) + (vil3d_I(+1,+0,-1)-vil3d_I(+0,+0,-1))*normx;
00244         xi20 = vil3d_I(+0,+1,-1) + (vil3d_I(+1,+1,-1)-vil3d_I(+0,+1,-1))*normx;
00245 
00246         xi11 = vil3d_I(+0,+0,+0) + (vil3d_I(+1,+0,+0)-vil3d_I(+0,+0,+0))*normx;
00247         xi21 = vil3d_I(+0,+1,+0) + (vil3d_I(+1,+1,+0)-vil3d_I(+0,+1,+0))*normx;
00248 
00249         xi12 = vil3d_I(+0,+0,+1) + (vil3d_I(+1,+0,+1)-vil3d_I(+0,+0,+1))*normx;
00250         xi22 = vil3d_I(+0,+1,+1) + (vil3d_I(+1,+1,+1)-vil3d_I(+0,+1,+1))*normx;
00251 
00252         xi13 = vil3d_I(+0,+0,+2) + (vil3d_I(+1,+0,+2)-vil3d_I(+0,+0,+2))*normx;
00253         xi23 = vil3d_I(+0,+1,+2) + (vil3d_I(+1,+1,+2)-vil3d_I(+0,+1,+2))*normx;
00254       }
00255       else
00256       {
00257         double s0, s1, s2, s3;
00258         get_tricubic_coeff( normx, s0, s1, s2, s3 );
00259 
00260         xi10 = 0.5 * ( s0*vil3d_I(-1,+0,-1) +
00261                        s1*vil3d_I(+0,+0,-1) +
00262                        s2*vil3d_I(+1,+0,-1) +
00263                        s3*vil3d_I(+2,+0,-1) );
00264         xi20 = 0.5 * ( s0*vil3d_I(-1,+1,-1) +
00265                        s1*vil3d_I(+0,+1,-1) +
00266                        s2*vil3d_I(+1,+1,-1) +
00267                        s3*vil3d_I(+2,+1,-1) );
00268 
00269         xi11 = 0.5 * ( s0*vil3d_I(-1,+0,+0) +
00270                        s1*vil3d_I(+0,+0,+0) +
00271                        s2*vil3d_I(+1,+0,+0) +
00272                        s3*vil3d_I(+2,+0,+0) );
00273         xi21 = 0.5 * ( s0*vil3d_I(-1,+1,+0) +
00274                        s1*vil3d_I(+0,+1,+0) +
00275                        s2*vil3d_I(+1,+1,+0) +
00276                        s3*vil3d_I(+2,+1,+0) );
00277 
00278         xi12 = 0.5 * ( s0*vil3d_I(-1,+0,+1) +
00279                        s1*vil3d_I(+0,+0,+1) +
00280                        s2*vil3d_I(+1,+0,+1) +
00281                        s3*vil3d_I(+2,+0,+1) );
00282         xi22 = 0.5 * ( s0*vil3d_I(-1,+1,+1) +
00283                        s1*vil3d_I(+0,+1,+1) +
00284                        s2*vil3d_I(+1,+1,+1) +
00285                        s3*vil3d_I(+2,+1,+1) );
00286 
00287         xi13 = 0.5 * ( s0*vil3d_I(-1,+0,+2) +
00288                        s1*vil3d_I(+0,+0,+2) +
00289                        s2*vil3d_I(+1,+0,+2) +
00290                        s3*vil3d_I(+2,+0,+2) );
00291         xi23 = 0.5 * ( s0*vil3d_I(-1,+1,+2) +
00292                        s1*vil3d_I(+0,+1,+2) +
00293                        s2*vil3d_I(+1,+1,+2) +
00294                        s3*vil3d_I(+2,+1,+2) );
00295       } // end if/else x
00296 
00297       val0 = xi10 + (xi20 - xi10)*normy;
00298       val1 = xi11 + (xi21 - xi11)*normy;
00299       val2 = xi12 + (xi22 - xi12)*normy;
00300       val3 = xi13 + (xi23 - xi13)*normy;
00301     } // end of if y
00302     else
00303     {
00304       if( x<1  || x>nx-3)
00305       {
00306         xi00 = vil3d_I(+0,-1,-1) + (vil3d_I(+1,-1,-1)-vil3d_I(+0,-1,-1))*normx;
00307         xi10 = vil3d_I(+0,+0,-1) + (vil3d_I(+1,+0,-1)-vil3d_I(+0,+0,-1))*normx;
00308         xi20 = vil3d_I(+0,+1,-1) + (vil3d_I(+1,+1,-1)-vil3d_I(+0,+1,-1))*normx;
00309         xi30 = vil3d_I(+0,+2,-1) + (vil3d_I(+1,+2,-1)-vil3d_I(+0,+2,-1))*normx;
00310 
00311         xi01 = vil3d_I(+0,-1,+0) + (vil3d_I(+1,-1,+0)-vil3d_I(+0,-1,+0))*normx;
00312         xi11 = vil3d_I(+0,+0,+0) + (vil3d_I(+1,+0,+0)-vil3d_I(+0,+0,+0))*normx;
00313         xi21 = vil3d_I(+0,+1,+0) + (vil3d_I(+1,+1,+0)-vil3d_I(+0,+1,+0))*normx;
00314         xi31 = vil3d_I(+0,+2,+0) + (vil3d_I(+1,+2,+0)-vil3d_I(+0,+2,+0))*normx;
00315 
00316         xi02 = vil3d_I(+0,-1,+1) + (vil3d_I(+1,-1,+1)-vil3d_I(+0,-1,+1))*normx;
00317         xi12 = vil3d_I(+0,+0,+1) + (vil3d_I(+1,+0,+1)-vil3d_I(+0,+0,+1))*normx;
00318         xi22 = vil3d_I(+0,+1,+1) + (vil3d_I(+1,+1,+1)-vil3d_I(+0,+1,+1))*normx;
00319         xi32 = vil3d_I(+0,+2,+1) + (vil3d_I(+1,+2,+1)-vil3d_I(+0,+2,+1))*normx;
00320 
00321         xi03 = vil3d_I(+0,-1,+2) + (vil3d_I(+1,-1,+2)-vil3d_I(+0,-1,+2))*normx;
00322         xi13 = vil3d_I(+0,+0,+2) + (vil3d_I(+1,+0,+2)-vil3d_I(+0,+0,+2))*normx;
00323         xi23 = vil3d_I(+0,+1,+2) + (vil3d_I(+1,+1,+2)-vil3d_I(+0,+1,+2))*normx;
00324         xi33 = vil3d_I(+0,+2,+2) + (vil3d_I(+1,+2,+2)-vil3d_I(+0,+2,+2))*normx;
00325       }
00326       else
00327       {
00328         double s0, s1, s2, s3;
00329         get_tricubic_coeff( normx, s0, s1, s2, s3 );
00330 
00331         xi00 = 0.5 * ( s0*vil3d_I(-1,-1,-1) +
00332                        s1*vil3d_I(+0,-1,-1) +
00333                        s2*vil3d_I(+1,-1,-1) +
00334                        s3*vil3d_I(+2,-1,-1) );
00335         xi10 = 0.5 * ( s0*vil3d_I(-1,+0,-1) +
00336                        s1*vil3d_I(+0,+0,-1) +
00337                        s2*vil3d_I(+1,+0,-1) +
00338                        s3*vil3d_I(+2,+0,-1) );
00339         xi20 = 0.5 * ( s0*vil3d_I(-1,+1,-1) +
00340                        s1*vil3d_I(+0,+1,-1) +
00341                        s2*vil3d_I(+1,+1,-1) +
00342                        s3*vil3d_I(+2,+1,-1) );
00343         xi30 = 0.5 * ( s0*vil3d_I(-1,+2,-1) +
00344                        s1*vil3d_I(+0,+2,-1) +
00345                        s2*vil3d_I(+1,+2,-1) +
00346                        s3*vil3d_I(+2,+2,-1) );
00347 
00348         xi01 = 0.5 * ( s0*vil3d_I(-1,-1,+0) +
00349                        s1*vil3d_I(+0,-1,+0) +
00350                        s2*vil3d_I(+1,-1,+0) +
00351                        s3*vil3d_I(+2,-1,+0) );
00352         xi11 = 0.5 * ( s0*vil3d_I(-1,+0,+0) +
00353                        s1*vil3d_I(+0,+0,+0) +
00354                        s2*vil3d_I(+1,+0,+0) +
00355                        s3*vil3d_I(+2,+0,+0) );
00356         xi21 = 0.5 * ( s0*vil3d_I(-1,+1,+0) +
00357                        s1*vil3d_I(+0,+1,+0) +
00358                        s2*vil3d_I(+1,+1,+0) +
00359                        s3*vil3d_I(+2,+1,+0) );
00360         xi31 = 0.5 * ( s0*vil3d_I(-1,+2,+0) +
00361                        s1*vil3d_I(+0,+2,+0) +
00362                        s2*vil3d_I(+1,+2,+0) +
00363                        s3*vil3d_I(+2,+2,+0) );
00364 
00365         xi02 = 0.5 * ( s0*vil3d_I(-1,-1,+1) +
00366                        s1*vil3d_I(+0,-1,+1) +
00367                        s2*vil3d_I(+1,-1,+1) +
00368                        s3*vil3d_I(+2,-1,+1) );
00369         xi12 = 0.5 * ( s0*vil3d_I(-1,+0,+1) +
00370                        s1*vil3d_I(+0,+0,+1) +
00371                        s2*vil3d_I(+1,+0,+1) +
00372                        s3*vil3d_I(+2,+0,+1) );
00373         xi22 = 0.5 * ( s0*vil3d_I(-1,+1,+1) +
00374                        s1*vil3d_I(+0,+1,+1) +
00375                        s2*vil3d_I(+1,+1,+1) +
00376                        s3*vil3d_I(+2,+1,+1) );
00377         xi32 = 0.5 * ( s0*vil3d_I(-1,+2,+1) +
00378                        s1*vil3d_I(+0,+2,+1) +
00379                        s2*vil3d_I(+1,+2,+1) +
00380                        s3*vil3d_I(+2,+2,+1) );
00381 
00382         xi03 = 0.5 * ( s0*vil3d_I(-1,-1,+2) +
00383                        s1*vil3d_I(+0,-1,+2) +
00384                        s2*vil3d_I(+1,-1,+2) +
00385                        s3*vil3d_I(+2,-1,+2) );
00386         xi13 = 0.5 * ( s0*vil3d_I(-1,+0,+2) +
00387                        s1*vil3d_I(+0,+0,+2) +
00388                        s2*vil3d_I(+1,+0,+2) +
00389                        s3*vil3d_I(+2,+0,+2) );
00390         xi23 = 0.5 * ( s0*vil3d_I(-1,+1,+2) +
00391                        s1*vil3d_I(+0,+1,+2) +
00392                        s2*vil3d_I(+1,+1,+2) +
00393                        s3*vil3d_I(+2,+1,+2) );
00394         xi33 = 0.5 * ( s0*vil3d_I(-1,+2,+2) +
00395                        s1*vil3d_I(+0,+2,+2) +
00396                        s2*vil3d_I(+1,+2,+2) +
00397                        s3*vil3d_I(+2,+2,+2) );
00398       } // end if/else x
00399 
00400       double t0, t1, t2, t3;
00401       get_tricubic_coeff( normy, t0, t1, t2, t3 );
00402       val0 = 0.5 * ( xi00*t0 + xi10*t1 + xi20*t2 + xi30*t3 );
00403       val1 = 0.5 * ( xi01*t0 + xi11*t1 + xi21*t2 + xi31*t3 );
00404       val2 = 0.5 * ( xi02*t0 + xi12*t1 + xi22*t2 + xi32*t3 );
00405       val3 = 0.5 * ( xi03*t0 + xi13*t1 + xi23*t2 + xi33*t3 );
00406     } // end of else y
00407 
00408     double u0, u1, u2, u3;
00409     get_tricubic_coeff( normz, u0, u1, u2, u3 );
00410 
00411     val = 0.5 * ( val0*u0 + val1*u1 + val2*u2 + val3*u3 );
00412   } //end of else z
00413 
00414   return val;
00415 }
00416 
00417 #define VIL3D_TRICUB_INTERP_INSTANTIATE(T)                              \
00418 template double vil3d_tricub_interp_raw (double x, double y, double z, const T* data,  \
00419                                          vcl_ptrdiff_t xstep, vcl_ptrdiff_t ystep, vcl_ptrdiff_t zstep); \
00420 template double vil3d_tricub_interp_safe_trilinear_extend (double x, double y, double z, const T* data,  \
00421                                                            int nx, int ny, int nz, \
00422                                                            vcl_ptrdiff_t xstep, vcl_ptrdiff_t ystep, vcl_ptrdiff_t zstep)
00423 
00424 #endif // vil3d_tricub_interp_txx_