BICPL: splines.h Source File

00001 #ifndef DEF_SPLINES 00002 #define DEF_SPLINES 00003 00004 /* ---------------------------------------------------------------------------- 00005 @COPYRIGHT : 00006 Copyright 1993,1994,1995 David MacDonald, 00007 McConnell Brain Imaging Centre, 00008 Montreal Neurological Institute, McGill University. 00009 Permission to use, copy, modify, and distribute this 00010 software and its documentation for any purpose and without 00011 fee is hereby granted, provided that the above copyright 00012 notice appear in all copies. The author and McGill University 00013 make no representations about the suitability of this 00014 software for any purpose. It is provided "as is" without 00015 express or implied warranty. 00016 ---------------------------------------------------------------------------- */ 00017 00018 #ifndef lint 00019 static char splines_rcsid[] = "$Header: /software/source//libraries/bicpl/Include/bicpl/splines.h,v 1.1 2000/10/07 10:40:37 stever Exp $"; 00020 #endif 00021 00022 #include <volume_io.h> 00023 00024 #define LINEAR_COEF_00 1.0 00025 #define LINEAR_COEF_01 0.0 00026 #define LINEAR_COEF_10 -1.0 00027 #define LINEAR_COEF_11 1.0 00028 00029 #define LINEAR_COEF_0( v0, v1 ) \ 00030 (v0) 00031 #define LINEAR_COEF_1( v0, v1 ) \ 00032 ((v1) - (v0)) 00033 00034 #define LINEAR_UNIVAR( v0, v1, u ) \ 00035 ( (v0) + (u) * ((v1) - (v0)) ) 00036 00037 /* ------------------ QUADRATIC INTERPOLATING SPLINES ----------------------- */ 00038 00039 #define QUADRATIC_COEF_00 0.5 00040 #define QUADRATIC_COEF_01 0.5 00041 #define QUADRATIC_COEF_02 0.0 00042 #define QUADRATIC_COEF_10 (-1.0) 00043 #define QUADRATIC_COEF_11 1.0 00044 #define QUADRATIC_COEF_12 0.0 00045 #define QUADRATIC_COEF_20 0.5 00046 #define QUADRATIC_COEF_21 (-1.0) 00047 #define QUADRATIC_COEF_22 0.5 00048 00049 #define QUADRATIC_COEF_OF_V0( u ) \ 00050 ( QUADRATIC_COEF_00 + (u) * ( \ 00051 QUADRATIC_COEF_10 + (u) * \ 00052 QUADRATIC_COEF_20) ) 00053 #define QUADRATIC_COEF_OF_V1( u ) \ 00054 ( QUADRATIC_COEF_01 + (u) * ( \ 00055 QUADRATIC_COEF_11 + (u) * \ 00056 QUADRATIC_COEF_21) ) 00057 #define QUADRATIC_COEF_OF_V2( u ) \ 00058 ( QUADRATIC_COEF_02 + (u) * ( \ 00059 QUADRATIC_COEF_12 + (u) * \ 00060 QUADRATIC_COEF_22) ) 00061 00062 #define COMPUTE_QUADRATIC_COEFFS( u, c0, c1, c2 ) \ 00063 { \ 00064 (c0) = QUADRATIC_COEF_OF_V0( u ); \ 00065 (c1) = QUADRATIC_COEF_OF_V1( u ); \ 00066 (c2) = QUADRATIC_COEF_OF_V2( u ); \ 00067 } 00068 00069 #define QUADRATIC_UNIVAR( v0, v1, v2, u ) \ 00070 ( (v0) * QUADRATIC_COEF_OF_V0(u) + \ 00071 (v1) * QUADRATIC_COEF_OF_V1(u) + \ 00072 (v2) * QUADRATIC_COEF_OF_V2(u) ) 00073 00074 #define QUADRATIC_DERIV_COEF_OF_V0( u ) \ 00075 ( QUADRATIC_COEF_10 + 2.0 * (u) * \ 00076 QUADRATIC_COEF_20 ) 00077 #define QUADRATIC_DERIV_COEF_OF_V1( u ) \ 00078 ( QUADRATIC_COEF_11 + 2.0 * (u) * \ 00079 QUADRATIC_COEF_21 ) 00080 #define QUADRATIC_DERIV_COEF_OF_V2( u ) \ 00081 ( QUADRATIC_COEF_12 + 2.0 * (u) * \ 00082 QUADRATIC_COEF_22 ) 00083 00084 #define COMPUTE_QUADRATIC_DERIV_COEFFS( u, c0, c1, c2 ) \ 00085 { \ 00086 (c0) = QUADRATIC_DERIV_COEF_OF_V0( u ); \ 00087 (c1) = QUADRATIC_DERIV_COEF_OF_V1( u ); \ 00088 (c2) = QUADRATIC_DERIV_COEF_OF_V2( u ); \ 00089 } 00090 00091 #define QUADRATIC_UNIVAR_DERIV( v0, v1, v2, u ) \ 00092 ( (v0) * QUADRATIC_DERIV_COEF_OF_V0(u) + \ 00093 (v1) * QUADRATIC_DERIV_COEF_OF_V1(u) + \ 00094 (v2) * QUADRATIC_DERIV_COEF_OF_V2(u) ) 00095 00096 #define QUADRATIC_DERIV2_COEF_OF_V0( u ) \ 00097 ( 2.0 * QUADRATIC_COEF_20 ) 00098 #define QUADRATIC_DERIV2_COEF_OF_V1( u ) \ 00099 ( 2.0 * QUADRATIC_COEF_21 ) 00100 #define QUADRATIC_DERIV2_COEF_OF_V2( u ) \ 00101 ( 2.0 * QUADRATIC_COEF_22 ) 00102 00103 #define COMPUTE_QUADRATIC_DERIV2_COEFFS( u, c0, c1, c2 ) \ 00104 { \ 00105 (c0) = QUADRATIC_DERIV2_COEF_OF_V0( u ); \ 00106 (c1) = QUADRATIC_DERIV2_COEF_OF_V1( u ); \ 00107 (c2) = QUADRATIC_DERIV2_COEF_OF_V2( u ); \ 00108 } 00109 00110 #define QUADRATIC_UNIVAR_DERIV2( v0, v1, v2, u ) \ 00111 ( (v0) * QUADRATIC_DERIV2_COEF_OF_V0(u) + \ 00112 (v1) * QUADRATIC_DERIV2_COEF_OF_V1(u) + \ 00113 (v2) * QUADRATIC_DERIV2_COEF_OF_V2(u) ) 00114 00115 #define DOT3( a1, b1, c1, a2, b2, c2 ) \ 00116 ( (a1) * (a2) + (b1) * (b2) + (c1) * (c2) ) 00117 00118 #define MULT3( r1, r2, r3, v1, v2, v3, v00, v01, v02, v10, v11, v12, v20, v21, v22 ) \ 00119 { (r1) = DOT3( v1, v2, v3, v00, v01, v02 ); \ 00120 (r2) = DOT3( v1, v2, v3, v10, v11, v12 ); \ 00121 (r3) = DOT3( v1, v2, v3, v20, v21, v22 ); \ 00122 } 00123 00124 #define QUADRATIC_BIVAR( v00, v01, v02, v10, v11, v12, v20, v21, v22, u_parm, v_parm, val ) \ 00125 { \ 00126 Real wv0, wv1, wv2; \ 00127 Real v0, v1, v2; \ 00128 \ 00129 COMPUTE_QUADRATIC_COEFFS( v_parm, wv0, wv1, wv2 ); \ 00130 \ 00131 MULT3( v0, v1, v2, wv0, wv1, wv2, \ 00132 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00133 \ 00134 (val) = QUADRATIC_UNIVAR( v0, v1, v2, u_parm ); \ 00135 } 00136 00137 #define QUADRATIC_BIVAR_DERIV( v00, v01, v02, v10, v11, v12, v20, v21, v22, u, v, val, du, dv ) \ 00138 { \ 00139 Real wu0, wu1, wu2; \ 00140 Real wv0, wv1, wv2; \ 00141 Real wdv0, wdv1, wdv2; \ 00142 Real v0, v1, v2; \ 00143 Real dv0, dv1, dv2; \ 00144 \ 00145 COMPUTE_QUADRATIC_COEFFS( u_parm, wu0, wu1, wu2 ); \ 00146 COMPUTE_QUADRATIC_COEFFS( v_parm, wv0, wv1, wv2 ); \ 00147 COMPUTE_QUADRATIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2 ); \ 00148 \ 00149 MULT3( v0, v1, v2, wv0, wv1, wv2, \ 00150 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00151 \ 00152 MULT3( dv0, dv1, dv2, wdv0, wdv1, wdv2, \ 00153 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00154 \ 00155 (dv) = DOT3( wu0, wu1, wu2, dv0, dv1, dv2 ); \ 00156 (du) = QUADRATIC_UNIVAR_DERIV( v0, v1, v2, u_parm ); \ 00157 (val) = DOT3( wu0, wu1, wu2, v0, v1, v2 ); \ 00158 } 00159 00160 #define QUADRATIC_BIVAR_DERIV2( v00, v01, v02, v10, v11, v12, v20, v21, v22, u_parm, v_parm, val, du, dv, duu, duv, dvv ) \ 00161 { \ 00162 Real wu0, wu1, wu2; \ 00163 Real wdu0, wdu1, wdu2; \ 00164 Real wv0, wv1, wv2; \ 00165 Real wdv0, wdv1, wdv2; \ 00166 Real wdvv0, wdvv1, wdvv2; \ 00167 Real v0, v1, v2; \ 00168 Real dv0, dv1, dv2; \ 00169 Real dvv0, dvv1, dvv2; \ 00170 \ 00171 COMPUTE_QUADRATIC_COEFFS( u_parm, wu0, wu1, wu2 ); \ 00172 COMPUTE_QUADRATIC_DERIV_COEFFS( u_parm, wdu0, wdu1, wdu2 ); \ 00173 COMPUTE_QUADRATIC_COEFFS( v_parm, wv0, wv1, wv2 ); \ 00174 COMPUTE_QUADRATIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2 ); \ 00175 COMPUTE_QUADRATIC_DERIV2_COEFFS( v_parm, wdvv0, wdvv1, wdvv2 ); \ 00176 \ 00177 MULT3( v0, v1, v2, wv0, wv1, wv2, \ 00178 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00179 MULT3( dv0, dv1, dv2, wdv0, wdv1, wdv2, \ 00180 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00181 MULT3( dvv0, dvv1, dvv2, wdvv0, wdvv1, wdvv2, \ 00182 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00183 \ 00184 (val) = DOT3( wu0, wu1, wu2, v0, v1, v2 ); \ 00185 (dv) = DOT3( wu0, wu1, wu2, dv0, dv1, dv2 ); \ 00186 (dvv) = DOT3( wu0, wu1, wu2, dvv0, dvv1, dvv2 ); \ 00187 (du) = DOT3( wdu0, wdu1, wdu2, v0, v1, v2 ); \ 00188 (duv) = DOT3( wdu0, wdu1, wdu2, dv0, dv1, dv2 ); \ 00189 (duu) = QUADRATIC_UNIVAR_DERIV2( v0, v1, v2, u_parm ); \ 00190 } 00191 00192 #define QUADRATIC_TRIVAR( c, u_parm, v_parm, w_parm, val ) \ 00193 { \ 00194 Real wv0, wv1, wv2; \ 00195 Real ww0, ww1, ww2; \ 00196 Real v00, v01, v02, v10, v11, v12, v20, v21, v22; \ 00197 Real v0, v1, v2; \ 00198 \ 00199 COMPUTE_QUADRATIC_COEFFS( v_parm, wv0, wv1, wv2 ); \ 00200 COMPUTE_QUADRATIC_COEFFS( w_parm, ww0, ww1, ww2 ); \ 00201 \ 00202 MULT3( v00, v01, v02, ww0, ww1, ww2, \ 00203 GLUE(c,000), GLUE(c,001), GLUE(c,002), \ 00204 GLUE(c,010), GLUE(c,011), GLUE(c,012), \ 00205 GLUE(c,020), GLUE(c,021), GLUE(c,022) ); \ 00206 MULT3( v10, v11, v12, ww0, ww1, ww2, \ 00207 GLUE(c,100), GLUE(c,101), GLUE(c,102), \ 00208 GLUE(c,110), GLUE(c,111), GLUE(c,112), \ 00209 GLUE(c,120), GLUE(c,121), GLUE(c,122) ); \ 00210 MULT3( v20, v21, v22, ww0, ww1, ww2, \ 00211 GLUE(c,200), GLUE(c,201), GLUE(c,202), \ 00212 GLUE(c,210), GLUE(c,211), GLUE(c,212), \ 00213 GLUE(c,220), GLUE(c,221), GLUE(c,222) ); \ 00214 \ 00215 MULT3( v0, v1, v2, wv0, wv1, wv2, \ 00216 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00217 \ 00218 (val) = QUADRATIC_UNIVAR( v0, v1, v2, u_parm ); \ 00219 } 00220 00221 #define QUADRATIC_TRIVAR_DERIV( c, u_parm, v_parm, w_parm, val, du, dv, dw ) \ 00222 { \ 00223 Real wu0, wu1, wu2; \ 00224 Real wv0, wv1, wv2; \ 00225 Real wdv0, wdv1, wdv2; \ 00226 Real ww0, ww1, ww2; \ 00227 Real wdw0, wdw1, wdw2; \ 00228 Real v00, v01, v02, v10, v11, v12, v20, v21, v22; \ 00229 Real dw00, dw01, dw02, dw10, dw11, dw12, dw20, dw21, dw22; \ 00230 Real v0, v1, v2; \ 00231 Real dv0, dv1, dv2; \ 00232 Real dw0, dw1, dw2; \ 00233 \ 00234 COMPUTE_QUADRATIC_COEFFS( u_parm, wu0, wu1, wu2 ); \ 00235 COMPUTE_QUADRATIC_COEFFS( v_parm, wv0, wv1, wv2 ); \ 00236 COMPUTE_QUADRATIC_COEFFS( w_parm, ww0, ww1, ww2 ); \ 00237 COMPUTE_QUADRATIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2 ); \ 00238 COMPUTE_QUADRATIC_DERIV_COEFFS( w_parm, wdw0, wdw1, wdw2 ); \ 00239 \ 00240 MULT3( v00, v01, v02, ww0, ww1, ww2, \ 00241 GLUE(c,000), GLUE(c,001), GLUE(c,002), \ 00242 GLUE(c,010), GLUE(c,011), GLUE(c,012), \ 00243 GLUE(c,020), GLUE(c,021), GLUE(c,022) ); \ 00244 MULT3( v10, v11, v12, ww0, ww1, ww2, \ 00245 GLUE(c,100), GLUE(c,101), GLUE(c,102), \ 00246 GLUE(c,110), GLUE(c,111), GLUE(c,112), \ 00247 GLUE(c,120), GLUE(c,121), GLUE(c,122) ); \ 00248 MULT3( v20, v21, v22, ww0, ww1, ww2, \ 00249 GLUE(c,200), GLUE(c,201), GLUE(c,202), \ 00250 GLUE(c,210), GLUE(c,211), GLUE(c,212), \ 00251 GLUE(c,220), GLUE(c,221), GLUE(c,222) ); \ 00252 \ 00253 MULT3( dw00, dw01, dw02, wdw0, wdw1, wdw2, \ 00254 GLUE(c,000), GLUE(c,001), GLUE(c,002), \ 00255 GLUE(c,010), GLUE(c,011), GLUE(c,012), \ 00256 GLUE(c,020), GLUE(c,021), GLUE(c,022) ); \ 00257 MULT3( dw10, dw11, dw12, wdw0, wdw1, wdw2, \ 00258 GLUE(c,100), GLUE(c,101), GLUE(c,102), \ 00259 GLUE(c,110), GLUE(c,111), GLUE(c,112), \ 00260 GLUE(c,120), GLUE(c,121), GLUE(c,122) ); \ 00261 MULT3( dw20, dw21, dw22, wdw0, wdw1, wdw2, \ 00262 GLUE(c,200), GLUE(c,201), GLUE(c,202), \ 00263 GLUE(c,210), GLUE(c,211), GLUE(c,212), \ 00264 GLUE(c,220), GLUE(c,221), GLUE(c,222) ); \ 00265 \ 00266 MULT3( v0, v1, v2, wv0, wv1, wv2, \ 00267 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00268 MULT3( dv0, dv1, dv2, wdv0, wdv1, wdv2, \ 00269 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00270 MULT3( dw0, dw1, dw2, wv0, wv1, wv2, \ 00271 dw00, dw01, dw02, dw10, dw11, dw12, dw20, dw21, dw22 ); \ 00272 \ 00273 (val) = DOT3( v0, v1, v2, wu0, wu1, wu2 ); \ 00274 (du) = QUADRATIC_UNIVAR_DERIV( v0, v1, v2, u_parm ); \ 00275 (dv) = DOT3( dv0, dv1, dv2, wu0, wu1, wu2 ); \ 00276 (dw) = DOT3( dw0, dw1, dw2, wu0, wu1, wu2 ); \ 00277 } 00278 00279 #define QUADRATIC_TRIVAR_DERIV2( c, u_parm, v_parm, w_parm, val, du, dv, dw, duu, duv, duw, dvv, dvw, dww ) \ 00280 { \ 00281 Real wu0, wu1, wu2; \ 00282 Real wdu0, wdu1, wdu2; \ 00283 Real wv0, wv1, wv2; \ 00284 Real wdv0, wdv1, wdv2; \ 00285 Real wdvv0, wdvv1, wdvv2; \ 00286 Real ww0, ww1, ww2; \ 00287 Real wdw0, wdw1, wdw2; \ 00288 Real wdww0, wdww1, wdww2; \ 00289 Real v00, v01, v02, v10, v11, v12, v20, v21, v22; \ 00290 Real dw00, dw01, dw02, dw10, dw11, dw12, dw20, dw21, dw22; \ 00291 Real dww00, dww01, dww02, dww10, dww11, dww12, dww20, dww21,dww22;\ 00292 Real v0, v1, v2; \ 00293 Real dv0, dv1, dv2; \ 00294 Real dvv0, dvv1, dvv2; \ 00295 Real dw0, dw1, dw2; \ 00296 Real dww0, dww1, dww2; \ 00297 Real dvw0, dvw1, dvw2; \ 00298 \ 00299 COMPUTE_QUADRATIC_COEFFS( u_parm, wu0, wu1, wu2 ); \ 00300 COMPUTE_QUADRATIC_COEFFS( v_parm, wv0, wv1, wv2 ); \ 00301 COMPUTE_QUADRATIC_COEFFS( w_parm, ww0, ww1, ww2 ); \ 00302 COMPUTE_QUADRATIC_DERIV_COEFFS( u_parm, wdu0, wdu1, wdu2 ); \ 00303 COMPUTE_QUADRATIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2 ); \ 00304 COMPUTE_QUADRATIC_DERIV_COEFFS( w_parm, wdw0, wdw1, wdw2 ); \ 00305 COMPUTE_QUADRATIC_DERIV2_COEFFS( v_parm, wdvv0, wdvv1, wdvv2 ); \ 00306 COMPUTE_QUADRATIC_DERIV2_COEFFS( w_parm, wdww0, wdww1, wdww2 ); \ 00307 \ 00308 MULT3( v00, v01, v02, ww0, ww1, ww2, \ 00309 GLUE(c,000), GLUE(c,001), GLUE(c,002), \ 00310 GLUE(c,010), GLUE(c,011), GLUE(c,012), \ 00311 GLUE(c,020), GLUE(c,021), GLUE(c,022) ); \ 00312 MULT3( v10, v11, v12, ww0, ww1, ww2, \ 00313 GLUE(c,100), GLUE(c,101), GLUE(c,102), \ 00314 GLUE(c,110), GLUE(c,111), GLUE(c,112), \ 00315 GLUE(c,120), GLUE(c,121), GLUE(c,122) ); \ 00316 MULT3( v20, v21, v22, ww0, ww1, ww2, \ 00317 GLUE(c,200), GLUE(c,201), GLUE(c,202), \ 00318 GLUE(c,210), GLUE(c,211), GLUE(c,212), \ 00319 GLUE(c,220), GLUE(c,221), GLUE(c,222) ); \ 00320 \ 00321 MULT3( dw00, dw01, dw02, wdw0, wdw1, wdw2, \ 00322 GLUE(c,000), GLUE(c,001), GLUE(c,002), \ 00323 GLUE(c,010), GLUE(c,011), GLUE(c,012), \ 00324 GLUE(c,020), GLUE(c,021), GLUE(c,022) ); \ 00325 MULT3( dw10, dw11, dw12, wdw0, wdw1, wdw2, \ 00326 GLUE(c,100), GLUE(c,101), GLUE(c,102), \ 00327 GLUE(c,110), GLUE(c,111), GLUE(c,112), \ 00328 GLUE(c,120), GLUE(c,121), GLUE(c,122) ); \ 00329 MULT3( dw20, dw21, dw22, wdw0, wdw1, wdw2, \ 00330 GLUE(c,200), GLUE(c,201), GLUE(c,202), \ 00331 GLUE(c,210), GLUE(c,211), GLUE(c,212), \ 00332 GLUE(c,220), GLUE(c,221), GLUE(c,222) ); \ 00333 \ 00334 MULT3( dww00, dww01, dww02, wdww0, wdww1, wdww2, \ 00335 GLUE(c,000), GLUE(c,001), GLUE(c,002), \ 00336 GLUE(c,010), GLUE(c,011), GLUE(c,012), \ 00337 GLUE(c,020), GLUE(c,021), GLUE(c,022) ); \ 00338 MULT3( dww10, dww11, dww12, wdww0, wdww1, wdww2, \ 00339 GLUE(c,100), GLUE(c,101), GLUE(c,102), \ 00340 GLUE(c,110), GLUE(c,111), GLUE(c,112), \ 00341 GLUE(c,120), GLUE(c,121), GLUE(c,122) ); \ 00342 MULT3( dww20, dww21, dww22, wdww0, wdww1, wdww2, \ 00343 GLUE(c,200), GLUE(c,201), GLUE(c,202), \ 00344 GLUE(c,210), GLUE(c,211), GLUE(c,212), \ 00345 GLUE(c,220), GLUE(c,221), GLUE(c,222) ); \ 00346 \ 00347 MULT3( v0, v1, v2, wv0, wv1, wv2, \ 00348 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00349 MULT3( dv0, dv1, dv2, wdv0, wdv1, wdv2, \ 00350 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00351 MULT3( dvv0, dvv1, dvv2, wdvv0, wdvv1, wdvv2, \ 00352 v00, v01, v02, v10, v11, v12, v20, v21, v22 ); \ 00353 MULT3( dw0, dw1, dw2, wv0, wv1, wv2, \ 00354 dw00, dw01, dw02, dw10, dw11, dw12, dw20, dw21, dw22 ); \ 00355 MULT3( dww0, dww1, dww2, wv0, wv1, wv2, \ 00356 dww00, dww01, dww02, dww10, dww11, dww12, dww20, dww21, dww22);\ 00357 MULT3( dvw0, dvw1, dvw2, wdv0, wdv1, wdv2, \ 00358 dw00, dw01, dw02, dw10, dw11, dw12, dw20, dw21, dw22 ); \ 00359 \ 00360 (val) = DOT3( v0, v1, v2, wu0, wu1, wu2 ); \ 00361 (dv) = DOT3( dv0, dv1, dv2, wu0, wu1, wu2 ); \ 00362 (dw) = DOT3( dw0, dw1, dw2, wu0, wu1, wu2 ); \ 00363 (dvv) = DOT3( dvv0, dvv1, dvv2, wu0, wu1, wu2 ); \ 00364 (dvw) = DOT3( dvw0, dvw1, dvw2, wu0, wu1, wu2 ); \ 00365 (dww) = DOT3( dww0, dww1, dww2, wu0, wu1, wu2 ); \ 00366 (du) = DOT3( v0, v1, v2, wdu0, wdu1, wdu2 ); \ 00367 (duv) = DOT3( dv0, dv1, dv2, wdu0, wdu1, wdu2 ); \ 00368 (duw) = DOT3( dw0, dw1, dw2, wdu0, wdu1, wdu2 ); \ 00369 (duu) = QUADRATIC_UNIVAR_DERIV2( v0, v1, v2, u_parm ); \ 00370 } 00371 00372 /* -------------------- CUBIC INTERPOLATING SPLINES ----------------------- */ 00373 00374 #define CUBIC_COEF_00 0.0 00375 #define CUBIC_COEF_01 1.0 00376 #define CUBIC_COEF_02 0.0 00377 #define CUBIC_COEF_03 0.0 00378 00379 #define CUBIC_COEF_10 (-0.5) 00380 #define CUBIC_COEF_11 0.0 00381 #define CUBIC_COEF_12 0.5 00382 #define CUBIC_COEF_13 0.0 00383 00384 #define CUBIC_COEF_20 1.0 00385 #define CUBIC_COEF_21 (-2.5) 00386 #define CUBIC_COEF_22 2.0 00387 #define CUBIC_COEF_23 (-0.5) 00388 00389 #define CUBIC_COEF_30 (-0.5) 00390 #define CUBIC_COEF_31 1.5 00391 #define CUBIC_COEF_32 (-1.5) 00392 #define CUBIC_COEF_33 0.5 00393 00394 #define CUBIC_COEF_OF_V0( u ) \ 00395 ( CUBIC_COEF_00 + (u) * ( \ 00396 CUBIC_COEF_10 + (u) * (\ 00397 CUBIC_COEF_20 + (u) * \ 00398 CUBIC_COEF_30)) ) 00399 #define CUBIC_COEF_OF_V1( u ) \ 00400 ( CUBIC_COEF_01 + (u) * ( \ 00401 CUBIC_COEF_11 + (u) * (\ 00402 CUBIC_COEF_21 + (u) * \ 00403 CUBIC_COEF_31)) ) 00404 #define CUBIC_COEF_OF_V2( u ) \ 00405 ( CUBIC_COEF_02 + (u) * ( \ 00406 CUBIC_COEF_12 + (u) * (\ 00407 CUBIC_COEF_22 + (u) * \ 00408 CUBIC_COEF_32)) ) 00409 #define CUBIC_COEF_OF_V3( u ) \ 00410 ( CUBIC_COEF_03 + (u) * ( \ 00411 CUBIC_COEF_13 + (u) * (\ 00412 CUBIC_COEF_23 + (u) * \ 00413 CUBIC_COEF_33)) ) 00414 00415 #define COMPUTE_CUBIC_COEFFS( u, c0, c1, c2, c3 ) \ 00416 { \ 00417 (c0) = CUBIC_COEF_OF_V0( u ); \ 00418 (c1) = CUBIC_COEF_OF_V1( u ); \ 00419 (c2) = CUBIC_COEF_OF_V2( u ); \ 00420 (c3) = CUBIC_COEF_OF_V3( u ); \ 00421 } 00422 00423 #define CUBIC_UNIVAR( v0, v1, v2, v3, u ) \ 00424 ( (v0) * CUBIC_COEF_OF_V0(u) + \ 00425 (v1) * CUBIC_COEF_OF_V1(u) + \ 00426 (v2) * CUBIC_COEF_OF_V2(u) + \ 00427 (v3) * CUBIC_COEF_OF_V3(u) ) 00428 00429 #define CUBIC_DERIV_COEF_OF_V0( u ) \ 00430 ( CUBIC_COEF_10 + (u) * \ 00431 (2.0 * CUBIC_COEF_20 + (u) * \ 00432 3.0 * CUBIC_COEF_30 ) ) 00433 #define CUBIC_DERIV_COEF_OF_V1( u ) \ 00434 ( CUBIC_COEF_11 + (u) * \ 00435 (2.0 * CUBIC_COEF_21 + (u) * \ 00436 3.0 * CUBIC_COEF_31 ) ) 00437 #define CUBIC_DERIV_COEF_OF_V2( u ) \ 00438 ( CUBIC_COEF_12 + (u) * \ 00439 (2.0 * CUBIC_COEF_22 + (u) * \ 00440 3.0 * CUBIC_COEF_32 ) ) 00441 #define CUBIC_DERIV_COEF_OF_V3( u ) \ 00442 ( CUBIC_COEF_13 + (u) * \ 00443 (2.0 * CUBIC_COEF_23 + (u) * \ 00444 3.0 * CUBIC_COEF_33 ) ) 00445 00446 #define COMPUTE_CUBIC_DERIV_COEFFS( u, c0, c1, c2, c3 ) \ 00447 { \ 00448 (c0) = CUBIC_DERIV_COEF_OF_V0( u ); \ 00449 (c1) = CUBIC_DERIV_COEF_OF_V1( u ); \ 00450 (c2) = CUBIC_DERIV_COEF_OF_V2( u ); \ 00451 (c3) = CUBIC_DERIV_COEF_OF_V3( u ); \ 00452 } 00453 00454 #define CUBIC_UNIVAR_DERIV( v0, v1, v2, v3, u ) \ 00455 ( (v0) * CUBIC_DERIV_COEF_OF_V0(u) + \ 00456 (v1) * CUBIC_DERIV_COEF_OF_V1(u) + \ 00457 (v2) * CUBIC_DERIV_COEF_OF_V2(u) + \ 00458 (v3) * CUBIC_DERIV_COEF_OF_V3(u) ) 00459 00460 #define CUBIC_DERIV2_COEF_OF_V0( u ) \ 00461 ( 2.0 * CUBIC_COEF_20 + (u) * \ 00462 6.0 * CUBIC_COEF_30 ) 00463 #define CUBIC_DERIV2_COEF_OF_V1( u ) \ 00464 ( 2.0 * CUBIC_COEF_21 + (u) * \ 00465 6.0 * CUBIC_COEF_31 ) 00466 #define CUBIC_DERIV2_COEF_OF_V2( u ) \ 00467 ( 2.0 * CUBIC_COEF_22 + (u) * \ 00468 6.0 * CUBIC_COEF_32 ) 00469 #define CUBIC_DERIV2_COEF_OF_V3( u ) \ 00470 ( 2.0 * CUBIC_COEF_23 + (u) * \ 00471 6.0 * CUBIC_COEF_33 ) 00472 00473 #define COMPUTE_CUBIC_DERIV2_COEFFS( u, c0, c1, c2, c3 ) \ 00474 { \ 00475 (c0) = CUBIC_DERIV2_COEF_OF_V0( u ); \ 00476 (c1) = CUBIC_DERIV2_COEF_OF_V1( u ); \ 00477 (c2) = CUBIC_DERIV2_COEF_OF_V2( u ); \ 00478 (c3) = CUBIC_DERIV2_COEF_OF_V3( u ); \ 00479 } 00480 00481 #define CUBIC_UNIVAR_DERIV2( v0, v1, v2, v3, u ) \ 00482 ( (v0) * CUBIC_DERIV2_COEF_OF_V0(u) + \ 00483 (v1) * CUBIC_DERIV2_COEF_OF_V1(u) + \ 00484 (v2) * CUBIC_DERIV2_COEF_OF_V2(u) + \ 00485 (v3) * CUBIC_DERIV2_COEF_OF_V3(u) ) 00486 00487 #define DOT4( a1, b1, c1, d1, a2, b2, c2, d2 ) \ 00488 ( (a1) * (a2) + (b1) * (b2) + (c1) * (c2) + (d1) * (d2) ) 00489 00490 #define MULT4( r1, r2, r3, r4, w1, w2, w3, w4, v ) \ 00491 { (r1) = DOT4(w1,w2,w3,w4,GLUE(v,00),GLUE(v,01),GLUE(v,02),GLUE(v,03));\ 00492 (r2) = DOT4(w1,w2,w3,w4,GLUE(v,10),GLUE(v,11),GLUE(v,12),GLUE(v,13));\ 00493 (r3) = DOT4(w1,w2,w3,w4,GLUE(v,20),GLUE(v,21),GLUE(v,22),GLUE(v,23));\ 00494 (r4) = DOT4(w1,w2,w3,w4,GLUE(v,30),GLUE(v,31),GLUE(v,32),GLUE(v,33));\ 00495 } 00496 00497 #define CUBIC_BIVAR(cv,u_parm,v_parm,val) \ 00498 { \ 00499 Real wv0, wv1, wv2, wv3; \ 00500 Real v0, v1, v2, v3; \ 00501 \ 00502 COMPUTE_CUBIC_COEFFS( v_parm, wv0, wv1, wv2, wv3 ); \ 00503 \ 00504 MULT4( v0, v1, v2, v3, wv0, wv1, wv2, wv3,cv); \ 00505 \ 00506 (val) = CUBIC_UNIVAR( v0, v1, v2, v3, u_parm ); \ 00507 } 00508 00509 #define CUBIC_BIVAR_DERIV( cv, u_parm, v_parm, val, du, dv ) \ 00510 { \ 00511 Real wu0, wu1, wu2, wu3; \ 00512 Real wv0, wv1, wv2, wv3; \ 00513 Real wdv0, wdv1, wdv2, wdv3; \ 00514 Real v0, v1, v2, v3; \ 00515 Real dv0, dv1, dv2, dv3; \ 00516 \ 00517 COMPUTE_CUBIC_COEFFS( u_parm, wu0, wu1, wu2, wu3 ); \ 00518 COMPUTE_CUBIC_COEFFS( v_parm, wv0, wv1, wv2, wv3 ); \ 00519 COMPUTE_CUBIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2, wdv3 ); \ 00520 \ 00521 MULT4( v0, v1, v2, v3, wv0, wv1, wv2, wv3,cv); \ 00522 MULT4( dv0, dv1, dv2, dv3, wdv0, wdv1, wdv2, wdv3,cv); \ 00523 \ 00524 (dv) = DOT4( wu0, wu1, wu2, wu3, dv0, dv1, dv2, dv3 ); \ 00525 (du) = CUBIC_UNIVAR_DERIV( v0, v1, v2, v3, u_parm ); \ 00526 (val) = DOT4( wu0, wu1, wu2, wu3, v0, v1, v2, v3 ); \ 00527 } 00528 00529 #define CUBIC_BIVAR_DERIV2( cv, u_parm, v_parm, val, du, dv, duu, duv, dvv ) \ 00530 { \ 00531 Real wu0, wu1, wu2, wu3; \ 00532 Real wdu0, wdu1, wdu2, wdu3; \ 00533 Real wv0, wv1, wv2, wv3; \ 00534 Real wdv0, wdv1, wdv2, wdv3; \ 00535 Real wdvv0, wdvv1, wdvv2, wdvv3; \ 00536 Real v0, v1, v2, v3; \ 00537 Real dv0, dv1, dv2, dv3; \ 00538 Real dvv0, dvv1, dvv2, dvv3; \ 00539 \ 00540 COMPUTE_CUBIC_COEFFS( u_parm, wu0, wu1, wu2, wu3 ); \ 00541 COMPUTE_CUBIC_DERIV_COEFFS( u_parm, wdu0, wdu1, wdu2, wdu3 ); \ 00542 COMPUTE_CUBIC_COEFFS( v_parm, wv0, wv1, wv2, wv3 ); \ 00543 COMPUTE_CUBIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2, wdv3 ); \ 00544 COMPUTE_CUBIC_DERIV2_COEFFS( v_parm, wdvv0, wdvv1, wdvv2, wdvv3 ); \ 00545 \ 00546 MULT4( v0, v1, v2, v3, wv0, wv1, wv2, wv3,cv); \ 00547 MULT4( dv0, dv1, dv2, dv3, wdv0, wdv1, wdv2, wdv3,cv); \ 00548 MULT4( dvv0, dvv1, dvv2, dvv3, wdvv0, wdvv1, wdvv2, wdvv3,cv); \ 00549 \ 00550 (val) = DOT4( wu0, wu1, wu2, wu3, v0, v1, v2, v3 ); \ 00551 (dv) = DOT4( wu0, wu1, wu2, wu3, dv0, dv1, dv2, dv3 ); \ 00552 (dvv) = DOT4( wu0, wu1, wu2, wu3, dvv0, dvv1, dvv2, dvv3 ); \ 00553 (du) = DOT4( wdu0, wdu1, wdu2, wdu3, v0, v1, v2, v3 ); \ 00554 (duv) = DOT4( wdu0, wdu1, wdu2, wdu3, dv0, dv1, dv2, dv3 ); \ 00555 (duu) = CUBIC_UNIVAR_DERIV2( v0, v1, v2, v3, u_parm ); \ 00556 } 00557 00558 #define CUBIC_TRIVAR( c, u_parm, v_parm, w_parm, val ) \ 00559 { \ 00560 Real wv0, wv1, wv2, wv3; \ 00561 Real ww0, ww1, ww2, ww3; \ 00562 Real v00, v01, v02, v03, v10, v11, v12, v13; \ 00563 Real v20, v21, v22, v23, v30, v31, v32, v33; \ 00564 Real v0, v1, v2, v3; \ 00565 \ 00566 COMPUTE_CUBIC_COEFFS( v_parm, wv0, wv1, wv2, wv3 ); \ 00567 COMPUTE_CUBIC_COEFFS( w_parm, ww0, ww1, ww2, ww3 ); \ 00568 \ 00569 MULT4( v00, v01, v02, v03, ww0, ww1, ww2, ww3,GLUE(c,0));\ 00570 MULT4( v10, v11, v12, v13, ww0, ww1, ww2, ww3,GLUE(c,1));\ 00571 MULT4( v20, v21, v22, v23, ww0, ww1, ww2, ww3,GLUE(c,2));\ 00572 MULT4( v30, v31, v32, v33, ww0, ww1, ww2, ww3,GLUE(c,3));\ 00573 \ 00574 MULT4( v0, v1, v2, v3, wv0, wv1, wv2, wv3,v); \ 00575 \ 00576 (val) = CUBIC_UNIVAR( v0, v1, v2, v3, u_parm ); \ 00577 } 00578 00579 #define CUBIC_TRIVAR_DERIV( c, u_parm, v_parm, w_parm, val, deriv_u, deriv_v, deriv_w ) \ 00580 { \ 00581 Real wu0, wu1, wu2, wu3; \ 00582 Real wv0, wv1, wv2, wv3; \ 00583 Real wdv0, wdv1, wdv2, wdv3; \ 00584 Real ww0, ww1, ww2, ww3; \ 00585 Real wdw0, wdw1, wdw2, wdw3; \ 00586 Real v00, v01, v02, v03, v10, v11, v12, v13; \ 00587 Real v20, v21, v22, v23, v30, v31, v32, v33; \ 00588 Real dw00, dw01, dw02, dw03, dw10, dw11, dw12, dw13; \ 00589 Real dw20, dw21, dw22, dw23, dw30, dw31, dw32, dw33; \ 00590 Real v0, v1, v2, v3; \ 00591 Real dv0, dv1, dv2, dv3; \ 00592 Real dw0, dw1, dw2, dw3; \ 00593 \ 00594 COMPUTE_CUBIC_COEFFS( u_parm, wu0, wu1, wu2, wu3 ); \ 00595 COMPUTE_CUBIC_COEFFS( v_parm, wv0, wv1, wv2, wv3 ); \ 00596 COMPUTE_CUBIC_COEFFS( w_parm, ww0, ww1, ww2, ww3 ); \ 00597 COMPUTE_CUBIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2, wdv3 ); \ 00598 COMPUTE_CUBIC_DERIV_COEFFS( w_parm, wdw0, wdw1, wdw2, wdw3 ); \ 00599 \ 00600 MULT4( v00, v01, v02, v03, ww0, ww1, ww2, ww3,GLUE(c,0)); \ 00601 MULT4( v10, v11, v12, v13, ww0, ww1, ww2, ww3,GLUE(c,1)); \ 00602 MULT4( v20, v21, v22, v23, ww0, ww1, ww2, ww3,GLUE(c,2)); \ 00603 MULT4( v30, v31, v32, v33, ww0, ww1, ww2, ww3,GLUE(c,3)); \ 00604 \ 00605 MULT4( dw00, dw01, dw02, dw03, wdw0, wdw1, wdw2, wdw3,GLUE(c,0));\ 00606 MULT4( dw10, dw11, dw12, dw13, wdw0, wdw1, wdw2, wdw3,GLUE(c,1));\ 00607 MULT4( dw20, dw21, dw22, dw23, wdw0, wdw1, wdw2, wdw3,GLUE(c,2));\ 00608 MULT4( dw30, dw31, dw32, dw33, wdw0, wdw1, wdw2, wdw3,GLUE(c,3));\ 00609 \ 00610 MULT4( v0, v1, v2, v3, wv0, wv1, wv2, wv3,v); \ 00611 MULT4( dv0, dv1, dv2, dv3, wdv0, wdv1, wdv2, wdv3,v);\ 00612 MULT4( dw0, dw1, dw2, dw3, wv0, wv1, wv2, wv3,dw);\ 00613 \ 00614 (val) = DOT4( v0, v1, v2, v3, wu0, wu1, wu2, wu3 ); \ 00615 (deriv_u) = CUBIC_UNIVAR_DERIV( v0, v1, v2, v3, u_parm ); \ 00616 (deriv_v) = DOT4( dv0, dv1, dv2, dv3, wu0, wu1, wu2, wu3 ); \ 00617 (deriv_w) = DOT4( dw0, dw1, dw2, dw3, wu0, wu1, wu2, wu3 ); \ 00618 } 00619 00620 #define CUBIC_TRIVAR_DERIV2( c, u_parm, v_parm, w_parm, val, deriv_u, deriv_v, deriv_w, deriv_uu, deriv_uv, deriv_uw, deriv_vv, deriv_vw, deriv_ww ) \ 00621 { \ 00622 Real wu0, wu1, wu2, wu3; \ 00623 Real wdu0, wdu1, wdu2, wdu3; \ 00624 Real wv0, wv1, wv2, wv3; \ 00625 Real wdv0, wdv1, wdv2, wdv3; \ 00626 Real wdvv0, wdvv1, wdvv2, wdvv3; \ 00627 Real ww0, ww1, ww2, ww3; \ 00628 Real wdw0, wdw1, wdw2, wdw3; \ 00629 Real wdww0, wdww1, wdww2, wdww3; \ 00630 Real v00, v01, v02, v03, v10, v11, v12, v13; \ 00631 Real v20, v21, v22, v23, v30, v31, v32, v33; \ 00632 Real dw00, dw01, dw02, dw03, dw10, dw11, dw12, dw13; \ 00633 Real dw20, dw21, dw22, dw23, dw30, dw31, dw32, dw33; \ 00634 Real dww00, dww01, dww02, dww03, dww10, dww11, dww12, dww13; \ 00635 Real dww20, dww21, dww22, dww23, dww30, dww31, dww32, dww33; \ 00636 Real v0, v1, v2, v3; \ 00637 Real dv0, dv1, dv2, dv3; \ 00638 Real dvv0, dvv1, dvv2, dvv3; \ 00639 Real dw0, dw1, dw2, dw3; \ 00640 Real dww0, dww1, dww2, dww3; \ 00641 Real dvw0, dvw1, dvw2, dvw3; \ 00642 \ 00643 COMPUTE_CUBIC_COEFFS( u_parm, wu0, wu1, wu2, wu3 ); \ 00644 COMPUTE_CUBIC_COEFFS( v_parm, wv0, wv1, wv2, wv3 ); \ 00645 COMPUTE_CUBIC_COEFFS( w_parm, ww0, ww1, ww2, ww3 ); \ 00646 COMPUTE_CUBIC_DERIV_COEFFS( u_parm, wdu0, wdu1, wdu2, wdu3 ); \ 00647 COMPUTE_CUBIC_DERIV_COEFFS( v_parm, wdv0, wdv1, wdv2, wdv3 ); \ 00648 COMPUTE_CUBIC_DERIV_COEFFS( w_parm, wdw0, wdw1, wdw2, wdw3 ); \ 00649 COMPUTE_CUBIC_DERIV2_COEFFS( v_parm, wdvv0, wdvv1, wdvv2, wdvv3 ); \ 00650 COMPUTE_CUBIC_DERIV2_COEFFS( w_parm, wdww0, wdww1, wdww2, wdww3 ); \ 00651 \ 00652 MULT4( v00, v01, v02, v03, ww0, ww1, ww2, ww3,GLUE(c,0)); \ 00653 MULT4( v10, v11, v12, v13, ww0, ww1, ww2, ww3,GLUE(c,1)); \ 00654 MULT4( v20, v21, v22, v23, ww0, ww1, ww2, ww3,GLUE(c,2)); \ 00655 MULT4( v30, v31, v32, v33, ww0, ww1, ww2, ww3,GLUE(c,3)); \ 00656 \ 00657 MULT4( dw00, dw01, dw02, dw03, wdw0, wdw1, wdw2, wdw3,GLUE(c,0)); \ 00658 MULT4( dw10, dw11, dw12, dw13, wdw0, wdw1, wdw2, wdw3,GLUE(c,1)); \ 00659 MULT4( dw20, dw21, dw22, dw23, wdw0, wdw1, wdw2, wdw3,GLUE(c,2)); \ 00660 MULT4( dw30, dw31, dw32, dw33, wdw0, wdw1, wdw2, wdw3,GLUE(c,3)); \ 00661 \ 00662 MULT4( dww00, dww01, dww02, dww03, wdww0, wdww1, wdww2, wdww3,GLUE(c,0)); \ 00663 MULT4( dww10, dww11, dww12, dww13, wdww0, wdww1, wdww2, wdww3,GLUE(c,1)); \ 00664 MULT4( dww20, dww21, dww22, dww23, wdww0, wdww1, wdww2, wdww3,GLUE(c,2)); \ 00665 MULT4( dww30, dww31, dww32, dww33, wdww0, wdww1, wdww2, wdww3,GLUE(c,3)); \ 00666 \ 00667 MULT4( v0, v1, v2, v3, wv0, wv1, wv2, wv3,v); \ 00668 MULT4( dv0, dv1, dv2, dv3, wdv0, wdv1, wdv2, wdv3,v); \ 00669 MULT4( dvv0, dvv1, dvv2, dvv3, wdvv0, wdvv1, wdvv2, wdvv3,v); \ 00670 MULT4( dw0, dw1, dw2, dw3, wv0, wv1, wv2, wv3,dw); \ 00671 MULT4( dww0, dww1, dww2, dww3, wv0, wv1, wv2, wv3,dww); \ 00672 MULT4( dvw0, dvw1, dvw2, dvw3, wdv0, wdv1, wdv2, wdv3,dw); \ 00673 \ 00674 (val) = DOT4( v0, v1, v2, v3, wu0, wu1, wu2, wu3 ); \ 00675 (deriv_v) = DOT4( dv0, dv1, dv2, dv3, wu0, wu1, wu2, wu3 ); \ 00676 (deriv_w) = DOT4( dw0, dw1, dw2, dw3, wu0, wu1, wu2, wu3 ); \ 00677 (deriv_vv) = DOT4( dvv0, dvv1, dvv2, dvv3, wu0, wu1, wu2, wu3 ); \ 00678 (deriv_vw) = DOT4( dvw0, dvw1, dvw2, dvw3, wu0, wu1, wu2, wu3 ); \ 00679 (deriv_ww) = DOT4( dww0, dww1, dww2, dww3, wu0, wu1, wu2, wu3 ); \ 00680 (deriv_u) = DOT4( v0, v1, v2, v3, wdu0, wdu1, wdu2, wdu3 ); \ 00681 (deriv_uv) = DOT4( dv0, dv1, dv2, dv3, wdu0, wdu1, wdu2, wdu3 ); \ 00682 (deriv_uw) = DOT4( dw0, dw1, dw2, dw3, wdu0, wdu1, wdu2, wdu3 ); \ 00683 (deriv_uu) = CUBIC_UNIVAR_DERIV2( v0, v1, v2, v3, u_parm ); \ 00684 } 00685 00686 #endif