00001 
00002 #include <stdlib.h>
00003 
00004 /*---------------------------------------------------------------------------
00005   Compute spherical harmonic expansion up to order 2 (quadratic terms)
00006   npt     : number of points to evaluate at
00007   wt      : coefficient vector [9]  (r**2 = x**2+y**2)
00008             [0] = 1
00009             [1] = x  [2] = y   [3] = z
00010             [4] = 2*z**2-r**2  [5] = xz [6] = yz  [7] = x**2-y**2 [8] = xy
00011   x, y, z : vectors of points at which to evaluate [npt]
00012   v       : output vector [npt]
00013 -----------------------------------------------------------------------------*/
00014 
00015 void warp3D_sharm2( int npt , float *wt ,
00016                     float *x , float *y , float *z , float *v )
00017 {
00018    int i ;
00019    float xq,yq,zq , a , b1,b2,b3 , c4,c5,c6,c7,c8 ;
00020 
00021    if( npt < 1 || wt == NULL || x == NULL ||
00022                   y  == NULL || z == NULL || v == NULL ) return ;
00023 
00024    a  = wt[0] ; b1 = wt[1] ; b2 = wt[2] ; b3 = wt[3] ;
00025    c4 = wt[4] ; c5 = wt[5] ; c6 = wt[6] ; c7 = wt[7] ; c8 = wt[8] ;
00026 
00027    for( i=0 ; i < npt ; i++ ){
00028       xq = x[i]*x[i] ; yq = y[i]*y[i] ; zq = z[i]*z[i] ;
00029       v[i] =  a + b2*y[i]
00030             + c4 * (2.0*zq-xq-yq)
00031             + c7 * (xq-yq)
00032             + (c5*x[i] + c6*y[i] + b3) * z[i]
00033             + (c8*y[i] + b1) * x[i] ;
00034    }
00035 }
00036 
00037 /*---------------------------------------------------------------------------*/
00038 
00039 void warp3D_sharm2_grad( int npt , float *wt ,
00040                          float *x , float *y , float *z ,
00041                          float *gx, float *gy, float *gz )
00042 {
00043    int i ;
00044    float b1,b2,b3 , c4,c5,c6,c7,c8 ;
00045    float gg,hh ;
00046 
00047    if( npt < 1   || wt == NULL || x == NULL ||
00048                     y  == NULL || z == NULL ||
00049       gx == NULL || gy == NULL || gz == NULL  ) return ;
00050 
00051                 b1 = wt[1] ; b2 = wt[2] ; b3 = wt[3] ;
00052    c4 = wt[4] ; c5 = wt[5] ; c6 = wt[6] ; c7 = wt[7] ; c8 = wt[8] ;
00053 
00054    gg = 2.0*(c7-c4) ; hh = -2.0*(c4+c7) ;
00055 
00056    for( i=0 ; i < npt ; i++ ){
00057       gx[i] = b1 + gg*x[i] + c8*y[i] + c5*z[i] ;
00058       gy[i] = b2 + c8*x[i] + hh*y[i] + c6*z[i] ;
00059       gz[i] = b3 + c5*x[i] + c6*y[i] + c4*z[i] ;
00060    }
00061 }
00062 
00063 
00064 /*---------------------------------------------------------------------------
00065   Compute spherical harmonic expansion up to order 3 (cubic terms)
00066   npt     : number of points to evaluate at
00067   wt      : coefficient vector [16]  (r**2 = x**2+y**2)
00068             [0] = 1
00069             [1] = x  [2] = y   [3] = z
00070             [4] = 2*z**2-r**2  [5] = xz [6] = yz  [7] = x**2-y**2 [8] = xy
00071             [9]  = z*(2z**2-3*r**2)
00072             [10] = x*(4*z**2-r**2)
00073             [11] = y*(4*z**2-r**2)
00074             [12] = z*(x**2-y**2)
00075             [13] = xyz
00076             [14] = x*(x**2-3*y**2)
00077             [15] = y*(3*x**2-y**2)
00078   x, y, z : vectors of points at which to evaluate [npt]
00079   v       : output vector [npt]
00080 -----------------------------------------------------------------------------*/
00081 
00082 void warp3D_sharm3( int npt , float *wt ,
00083                     float *x , float *y , float *z , float *v )
00084 {
00085    int i ;
00086    float xq,yq,zq,rr , a , b1,b2,b3 , c4,c5,c6,c7,c8 , d9,d10,d11,d12,d13,d14,d15;
00087 
00088    if( npt < 1 || wt == NULL || x == NULL ||
00089                   y  == NULL || z == NULL || v == NULL ) return ;
00090 
00091    a  = wt[0] ; b1 = wt[1] ; b2 = wt[2] ; b3 = wt[3] ;
00092    c4 = wt[4] ; c5 = wt[5] ; c6 = wt[6] ; c7 = wt[7] ; c8 = wt[8] ;
00093 
00094    d9  = wt[9]  ; d10 = wt[10] ; d11 = wt[11] ; d12 = wt[12] ;
00095    d13 = wt[13] ; d14 = wt[14] ; d15 = wt[15] ;
00096 
00097    for( i=0 ; i < npt ; i++ ){
00098       xq = x[i]*x[i] ; yq = y[i]*y[i] ; zq = z[i]*z[i] ; rr = xq+yq ;
00099       v[i] =  a
00100             + ( b2 + d15 * (3.0*xq-yq) ) * y[i]
00101             + c4 * (2.0*zq-xq-yq)
00102             + (c7 + d12*z[i] ) * (xq-yq)
00103             + (c5*x[i] + c6*y[i] + b3 + d9*(2.0*zq-3.0*rr) ) * z[i]
00104             + (c8*y[i] + b1 + d13*y[i]*z[i] + d14*(xq-3.0*yq) ) * x[i]
00105             + (d10*x[i] + d11*y[i])*(4.0*zq-rr) ;
00106    }
00107 }
00108 
00109 /*---------------------------------------------------------------------------
00110   Code for this was generated by Maple's codegen facility
00111 -----------------------------------------------------------------------------*/
00112 
00113 void warp3D_sharm3_grad( int npt , float *wt ,
00114                          float *x , float *y , float *z ,
00115                          float *gx, float *gy, float *gz )
00116 {
00117    int i ;
00118    float xq,yq,zq,rr ,  b1,b2,b3 , c4,c5,c6,c7,c8 , d9,d10,d11,d12,d13,d14,d15;
00119    float df0,df1,df2,df3 , t9,t21,t28,t29,t32,t34,t40,t41,t43 ;
00120 
00121    if( npt < 1 || wt == NULL || x  == NULL ||
00122                   y  == NULL || z  == NULL ||
00123                   gx == NULL || gy == NULL || gz == NULL ) return ;
00124 
00125                 b1 = wt[1] ; b2 = wt[2] ; b3 = wt[3] ;
00126    c4 = wt[4] ; c5 = wt[5] ; c6 = wt[6] ; c7 = wt[7] ; c8 = wt[8] ;
00127 
00128    d9  = wt[9]  ; d10 = wt[10] ; d11 = wt[11] ; d12 = wt[12] ;
00129    d13 = wt[13] ; d14 = wt[14] ; d15 = wt[15] ;
00130 
00131    for( i=0 ; i < npt ; i++ ){
00132      xq = x[i]*x[i];
00133      yq = y[i]*y[i];
00134      zq = z[i]*z[i];
00135      rr = xq+yq;
00136      t9 = d12*z[i];
00137      t21 = d13*y[i];
00138      t28 = d10*x[i];
00139      t29 = d11*y[i];
00140      t32 = 4.0*zq-rr;
00141      t34 = d9*z[i];
00142      df3 = 3.0*t34-t28-t29;
00143      df2 = 2.0*(c4+t34)+4.0*(t28+t29);
00144      t40 = d15*y[i];
00145      t41 = d14*x[i];
00146      t43 = df3;
00147      df1 = -t40-c4-c7-t9-3.0*t41+t43;
00148      df0 = 3.0*t40-c4+c7+t9+t41+t43;
00149      gx[i] = c5*z[i]+c8*y[i]+b1+t21*z[i]+d14*(xq-3.0*yq)+d10*t32+2.0*df0*x[i];
00150      gy[i] = b2+d15*(3.0*xq-yq)+c6*z[i]+(c8+d13*z[i])*x[i]+d11*t32+2.0*df1*y[i];
00151      gz[i] = d12*(xq-yq)+c5*x[i]+c6*y[i]+b3+d9*(2.0*zq-3.0*rr)+t21*x[i]+2.0*df2*z[i];
00152   }
00153 }