00001 #include "mrilib.h"
00002 
00003 /*** NOT 7D SAFE ***/
00004 
00005 /****************************************************************************/
00006 /** functions to "warp" 3D images -- not very efficient, but quite general **/
00007 /****************************************************************************/
00008 
00009 /*--------------------------------------------------------------------------*/
00010 
00011 static int wtype = MRI_LINEAR ;
00012 void mri_warp3D_method( int mode ){ wtype = mode ; } /** set interpolation **/
00013 
00014 /*--------------------------------------------------------------------------*/
00015 
00016 static byte *womask = NULL ;  /* 19 Nov 2004 */
00017 
00018 void mri_warp3D_set_womask( MRI_IMAGE *wim )
00019 {
00020    womask = (wim == NULL || wim->kind != MRI_byte) ? NULL
00021                                                    : MRI_BYTE_PTR(wim) ;
00022 }
00023 
00024 #define SKIP(i,j,k) (womask!=NULL && womask[(i)+(j)*nxnew+(k)*nxynew]==0)
00025 
00026 /*--------------------------------------------------------------------------*/
00027 
00028 static int zout = 1 ;
00029 void mri_warp3D_zerout( int zzzz ){ zout  = zzzz ; }      /** outside = 0? **/
00030 
00031 /*--------------------------------------------------------------------------*/
00032 /* Macros for interpolation and access to data arrays. */
00033 
00034 #define FAR(i,j,k) far[(i)+(j)*nx+(k)*nxy]         /* input point */
00035 #define NAR(i,j,k) nar[(i)+(j)*nxnew+(k)*nxynew]   /* output point */
00036 
00037 /** clip input mm to range [0..nn] (for indexing) **/
00038 
00039 #define CLIP(mm,nn) if(mm < 0)mm=0; else if(mm > nn)mm=nn
00040 
00041 /* define Lagrange cubic interpolation polynomials */
00042 
00043 #define P_M1(x)  (-(x)*((x)-1)*((x)-2))
00044 #define P_00(x)  (3*((x)+1)*((x)-1)*((x)-2))
00045 #define P_P1(x)  (-3*(x)*((x)+1)*((x)-2))
00046 #define P_P2(x)  ((x)*((x)+1)*((x)-1))
00047 #define P_FACTOR 4.62962963e-3            /* 1/216 = final scaling factor */
00048 
00049 /*---------------------------------------------------------------------------*/
00050 /*! Transform a 3D image geometrically, using cubic interpolation.
00051     See comments for mri_warp3D() for details about the parameters.
00052 -----------------------------------------------------------------------------*/
00053 
00054 INLINE MRI_IMAGE *mri_warp3D_cubic(
00055                     MRI_IMAGE *im, int nxnew, int nynew, int nznew,
00056                     void wf(float,float,float,float *,float *,float *) )
00057 {
00058    MRI_IMAGE *imfl , *newImg ;
00059    float *far , *nar ;
00060    float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
00061    int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ;
00062    float bot,top,val ;
00063    float nxh,nyh,nzh ;
00064 
00065    int ix_m1,ix_00,ix_p1,ix_p2 ;         /* interpolation indices */
00066    int jy_m1,jy_00,jy_p1,jy_p2 ;         /* (input image) */
00067    int kz_m1,kz_00,kz_p1,kz_p2 ;
00068 
00069    float wt_m1,wt_00,wt_p1,wt_p2 ;   /* interpolation weights */
00070 
00071    float f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, /* interpolants */
00072          f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00,
00073          f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1,
00074          f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2,
00075          f_km1    , f_k00    , f_kp1    , f_kp2     ;
00076 
00077 ENTRY("mri_warp3D_cubic") ;
00078 
00079    /*--- sanity check inputs ---*/
00080 
00081    if( im == NULL || wf == NULL ) RETURN(NULL) ;
00082 
00083    /*-- dimensional analysis --*/
00084 
00085    nx = im->nx ; nx1 = nx-1 ;          /* input image dimensions */
00086    ny = im->ny ; ny1 = ny-1 ;
00087    nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ;
00088 
00089    nxnew = (nxnew > 0) ? nxnew : nx ;  /* output image dimensions */
00090    nynew = (nynew > 0) ? nynew : ny ;
00091    nznew = (nznew > 0) ? nznew : nz ;
00092    nxynew = nxnew*nynew ;
00093 
00094    /*----- Allow for different input image types, by breaking them into
00095            components, doing each one separately, and then reassembling.
00096            This is needed since we only warp float-valued images below.  -----*/
00097 
00098    switch( im->kind ){
00099 
00100      case MRI_float:                              /* use input directly */
00101        imfl = im ; break ;
00102 
00103      default:{                                    /* floatize input */
00104        imfl = mri_to_float(im) ;
00105        newImg  = mri_warp3D_cubic( imfl , nxnew,nynew,nznew , wf ) ;
00106        mri_free(imfl) ;
00107        imfl = mri_to_mri(im->kind,newImg) ;
00108        if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
00109        RETURN(newImg) ;
00110      }
00111 
00112      case MRI_rgb:{                                /* break into 3 pieces */
00113        MRI_IMARR *imar = mri_rgb_to_3float(im) ;
00114        MRI_IMAGE *rim,*gim,*bim ;
00115        rim = mri_warp3D_cubic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00116                      mri_free( IMARR_SUBIM(imar,0) ) ;
00117        gim = mri_warp3D_cubic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00118                      mri_free( IMARR_SUBIM(imar,1) ) ;
00119        bim = mri_warp3D_cubic( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
00120                      mri_free( IMARR_SUBIM(imar,2) ) ;
00121        FREE_IMARR(imar) ;
00122        newImg = mri_3to_rgb( rim,gim,bim ) ;
00123        mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
00124      }
00125 
00126      case MRI_complex:{                             /* break into 2 pieces */
00127        MRI_IMARR *imar = mri_complex_to_pair(im) ;
00128        MRI_IMAGE *rim, *iim ;
00129        rim = mri_warp3D_cubic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00130                      mri_free( IMARR_SUBIM(imar,0) ) ;
00131        iim = mri_warp3D_cubic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00132                      mri_free( IMARR_SUBIM(imar,1) ) ;
00133        FREE_IMARR(imar) ;
00134        newImg = mri_pair_to_complex( rim , iim ) ;
00135        mri_free(rim); mri_free(iim); RETURN(newImg);
00136      }
00137 
00138    } /* end of special cases of input datum */
00139 
00140    /*----- at this point, imfl is in float format -----*/
00141 
00142    far = MRI_FLOAT_PTR( imfl ) ;                         /* input image data */
00143 
00144    newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ;  /* make output image */
00145    nar = MRI_FLOAT_PTR( newImg ) ;                         /* output image data */
00146 
00147    bot = top = far[0] ;                             /* find input data range */
00148    for( ii=1 ; ii < imfl->nvox ; ii++ ){
00149           if( far[ii] > top ) top = far[ii] ;
00150      else if( far[ii] < bot ) bot = far[ii] ;
00151    }
00152 
00153    /*** loop over output points and warp to them ***/
00154 
00155    nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
00156    for( kk=0 ; kk < nznew ; kk++ ){
00157     zpr = kk ;
00158     for( jj=0 ; jj < nynew ; jj++ ){
00159      ypr = jj ;
00160      for( ii=0 ; ii < nxnew ; ii++ ){
00161        xpr = ii ;
00162        if( SKIP(ii,jj,kk) ) continue ;    /* 19 Nov 2004 */
00163        wf( xpr,ypr,zpr , &xx,&yy,&zz ) ;  /* get xx,yy,zz in original image */
00164 
00165        if( zout &&
00166            (xx < -0.5 || xx > nxh ||
00167             yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){
00168 
00169           NAR(ii,jj,kk) = 0.0 ; continue ;
00170        }
00171 
00172        ix = floor(xx) ;  fx = xx - ix ;   /* integer and       */
00173        jy = floor(yy) ;  fy = yy - jy ;   /* fractional coords */
00174        kz = floor(zz) ;  fz = zz - kz ;
00175 
00176        ix_m1 = ix-1    ; ix_00 = ix      ; ix_p1 = ix+1    ; ix_p2 = ix+2    ;
00177        CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ;
00178 
00179        jy_m1 = jy-1    ; jy_00 = jy      ; jy_p1 = jy+1    ; jy_p2 = jy+2    ;
00180        CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ;
00181 
00182        kz_m1 = kz-1    ; kz_00 = kz      ; kz_p1 = kz+1    ; kz_p2 = kz+2    ;
00183        CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ;
00184 
00185        wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ;  /* interpolation weights */
00186        wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
00187 
00188 #undef  XINT
00189 #define XINT(j,k) wt_m1*FAR(ix_m1,j,k)+wt_00*FAR(ix_00,j,k)  \
00190                  +wt_p1*FAR(ix_p1,j,k)+wt_p2*FAR(ix_p2,j,k)
00191 
00192        /* interpolate to location ix+fx at each jy,kz level */
00193 
00194        f_jm1_km1 = XINT(jy_m1,kz_m1) ; f_j00_km1 = XINT(jy_00,kz_m1) ;
00195        f_jp1_km1 = XINT(jy_p1,kz_m1) ; f_jp2_km1 = XINT(jy_p2,kz_m1) ;
00196        f_jm1_k00 = XINT(jy_m1,kz_00) ; f_j00_k00 = XINT(jy_00,kz_00) ;
00197        f_jp1_k00 = XINT(jy_p1,kz_00) ; f_jp2_k00 = XINT(jy_p2,kz_00) ;
00198        f_jm1_kp1 = XINT(jy_m1,kz_p1) ; f_j00_kp1 = XINT(jy_00,kz_p1) ;
00199        f_jp1_kp1 = XINT(jy_p1,kz_p1) ; f_jp2_kp1 = XINT(jy_p2,kz_p1) ;
00200        f_jm1_kp2 = XINT(jy_m1,kz_p2) ; f_j00_kp2 = XINT(jy_00,kz_p2) ;
00201        f_jp1_kp2 = XINT(jy_p1,kz_p2) ; f_jp2_kp2 = XINT(jy_p2,kz_p2) ;
00202 
00203        /* interpolate to jy+fy at each kz level */
00204 
00205        wt_m1 = P_M1(fy) ; wt_00 = P_00(fy) ;
00206        wt_p1 = P_P1(fy) ; wt_p2 = P_P2(fy) ;
00207 
00208        f_km1 =  wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1
00209               + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 ;
00210 
00211        f_k00 =  wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00
00212               + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 ;
00213 
00214        f_kp1 =  wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1
00215               + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 ;
00216 
00217        f_kp2 =  wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2
00218               + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 ;
00219 
00220        /* interpolate to kz+fz to get output */
00221 
00222        wt_m1 = P_M1(fz) ; wt_00 = P_00(fz) ;
00223        wt_p1 = P_P1(fz) ; wt_p2 = P_P2(fz) ;
00224 
00225        val = P_FACTOR * (  wt_m1 * f_km1 + wt_00 * f_k00
00226                          + wt_p1 * f_kp1 + wt_p2 * f_kp2 ) ;
00227 
00228             if( val > top ) val = top ;   /* clip to input data range */
00229        else if( val < bot ) val = bot ;
00230 
00231        NAR(ii,jj,kk) = val ;
00232      }
00233     }
00234    }
00235 
00236    /*** cleanup and return ***/
00237 
00238    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00239    RETURN(newImg);
00240 }
00241 
00242 /*---------------------------------------------------------------------------*/
00243 /*! Transform a 3D image geometrically, using linear interpolation.
00244     See comments for mri_warp3D() for details about the parameters.
00245 -----------------------------------------------------------------------------*/
00246 
00247 INLINE MRI_IMAGE *mri_warp3D_linear(
00248                     MRI_IMAGE *im, int nxnew, int nynew, int nznew,
00249                     void wf(float,float,float,float *,float *,float *) )
00250 {
00251    MRI_IMAGE *imfl , *newImg ;
00252    float *far , *nar ;
00253    float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
00254    int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ;
00255    float val ;
00256    float nxh,nyh,nzh ;
00257 
00258    int ix_00,ix_p1 ;         /* interpolation indices */
00259    int jy_00,jy_p1 ;         /* (input image) */
00260    int kz_00,kz_p1 ;
00261 
00262    float wt_00,wt_p1 ;   /* interpolation weights */
00263 
00264    float f_j00_k00, f_jp1_k00, f_j00_kp1, f_jp1_kp1, f_k00, f_kp1 ;
00265 
00266 ENTRY("mri_warp3D_linear") ;
00267 
00268    /*--- sanity check inputs ---*/
00269 
00270    if( im == NULL || wf == NULL ) RETURN(NULL) ;
00271 
00272    /*-- dimensional analysis --*/
00273 
00274    nx = im->nx ; nx1 = nx-1 ;          /* input image dimensions */
00275    ny = im->ny ; ny1 = ny-1 ;
00276    nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ;
00277 
00278    nxnew = (nxnew > 0) ? nxnew : nx ;  /* output image dimensions */
00279    nynew = (nynew > 0) ? nynew : ny ;
00280    nznew = (nznew > 0) ? nznew : nz ;
00281    nxynew = nxnew*nynew ;
00282 
00283    /*----- allow for different input image types, by breaking them into
00284            components, doing each one separately, and then reassembling -----*/
00285 
00286    switch( im->kind ){
00287 
00288      case MRI_float:                        /* use input directly */
00289        imfl = im ; break ;
00290 
00291      default:{                              /* floatize-input */
00292        imfl = mri_to_float(im) ;
00293        newImg  = mri_warp3D_linear( imfl , nxnew,nynew,nznew , wf ) ;
00294        mri_free(imfl) ;
00295        imfl = mri_to_mri(im->kind,newImg) ;
00296        if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
00297        RETURN(newImg) ;
00298      }
00299 
00300      case MRI_rgb:{
00301        MRI_IMARR *imar = mri_rgb_to_3float(im) ;
00302        MRI_IMAGE *rim,*gim,*bim ;
00303        rim = mri_warp3D_linear( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00304                       mri_free( IMARR_SUBIM(imar,0) ) ;
00305        gim = mri_warp3D_linear( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00306                       mri_free( IMARR_SUBIM(imar,1) ) ;
00307        bim = mri_warp3D_linear( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
00308                       mri_free( IMARR_SUBIM(imar,2) ) ;
00309        FREE_IMARR(imar) ;
00310        newImg = mri_3to_rgb( rim,gim,bim ) ;
00311        mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
00312      }
00313 
00314      case MRI_complex:{
00315        MRI_IMARR *imar = mri_complex_to_pair(im) ;
00316        MRI_IMAGE *rim, *iim ;
00317        rim = mri_warp3D_linear( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00318                       mri_free( IMARR_SUBIM(imar,0) ) ;
00319        iim = mri_warp3D_linear( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00320                       mri_free( IMARR_SUBIM(imar,1) ) ;
00321        FREE_IMARR(imar) ;
00322        newImg = mri_pair_to_complex( rim , iim ) ;
00323        mri_free(rim); mri_free(iim); RETURN(newImg);
00324      }
00325 
00326    } /* end of special cases of input datum */
00327 
00328    /*----- at this point, imfl is in float format -----*/
00329 
00330    far = MRI_FLOAT_PTR( imfl ) ;                         /* input image data */
00331 
00332    newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ;  /* make output image */
00333    nar = MRI_FLOAT_PTR( newImg ) ;                         /* output image data */
00334 
00335    /*** loop over output points and warp to them ***/
00336 
00337    nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
00338    for( kk=0 ; kk < nznew ; kk++ ){
00339     zpr = kk ;
00340     for( jj=0 ; jj < nynew ; jj++ ){
00341      ypr = jj ;
00342      for( ii=0 ; ii < nxnew ; ii++ ){
00343        xpr = ii ;
00344        if( SKIP(ii,jj,kk) ) continue ;    /* 19 Nov 2004 */
00345        wf( xpr,ypr,zpr , &xx,&yy,&zz ) ;  /* get xx,yy,zz in original image */
00346 
00347        if( zout &&
00348            (xx < -0.5 || xx > nxh ||
00349             yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){
00350 
00351           NAR(ii,jj,kk) = 0.0 ; continue ;
00352        }
00353 
00354        ix = floor(xx) ;  fx = xx - ix ;   /* integer and       */
00355        jy = floor(yy) ;  fy = yy - jy ;   /* fractional coords */
00356        kz = floor(zz) ;  fz = zz - kz ;
00357 
00358        ix_00 = ix      ; ix_p1 = ix+1    ;
00359        CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ;
00360 
00361        jy_00 = jy      ; jy_p1 = jy+1    ;
00362        CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ;
00363 
00364        kz_00 = kz      ; kz_p1 = kz+1    ;
00365        CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ;
00366 
00367        wt_00 = 1.0-fx ; wt_p1 = fx ;
00368 
00369 #undef  XINT
00370 #define XINT(j,k) wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k)
00371 
00372        /* interpolate to location ix+fx at each jy,kz level */
00373 
00374        f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ;
00375        f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ;
00376 
00377        /* interpolate to jy+fy at each kz level */
00378 
00379        wt_00 = 1.0-fy ; wt_p1 = fy ;
00380 
00381        f_k00 =  wt_00 * f_j00_k00 + wt_p1 * f_jp1_k00 ;
00382 
00383        f_kp1 =  wt_00 * f_j00_kp1 + wt_p1 * f_jp1_kp1 ;
00384 
00385        /* interpolate to kz+fz to get output */
00386 
00387        val = (1.0-fz) * f_k00 + fz * f_kp1 ;
00388 
00389        NAR(ii,jj,kk) = val ;
00390      }
00391     }
00392    }
00393 
00394    /*** cleanup and return ***/
00395 
00396    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00397    RETURN(newImg);
00398 }
00399 
00400 /*---------------------------------------------------------------------------*/
00401 /*! Transform a 3D image geometrically, using NN 'interpolation'.
00402     See comments for mri_warp3D() for details about the parameters.
00403 -----------------------------------------------------------------------------*/
00404 
00405 INLINE MRI_IMAGE *mri_warp3D_NN(
00406              MRI_IMAGE *im, int nxnew, int nynew, int nznew,
00407              void wf(float,float,float,float *,float *,float *) )
00408 {
00409    MRI_IMAGE *imfl , *newImg ;
00410    float *far , *nar ;
00411    float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
00412    int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ;
00413    float nxh,nyh,nzh ;
00414 
00415 ENTRY("mri_warp3D_NN") ;
00416 
00417    /*--- sanity check inputs ---*/
00418 
00419    if( im == NULL || wf == NULL ) RETURN(NULL) ;
00420 
00421    /*-- dimensional analysis --*/
00422 
00423    nx = im->nx ; nx1 = nx-1 ;          /* input image dimensions */
00424    ny = im->ny ; ny1 = ny-1 ;
00425    nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ;
00426 
00427    nxnew = (nxnew > 0) ? nxnew : nx ;  /* output image dimensions */
00428    nynew = (nynew > 0) ? nynew : ny ;
00429    nznew = (nznew > 0) ? nznew : nz ;
00430    nxynew = nxnew*nynew ;
00431 
00432    /*----- allow for different input image types, by breaking them into
00433            components, doing each one separately, and then reassembling -----*/
00434 
00435    switch( im->kind ){
00436 
00437      case MRI_float:                        /* use input directly */
00438        imfl = im ; break ;
00439 
00440      default:{                              /* floatize-input */
00441        imfl = mri_to_float(im) ;
00442        newImg  = mri_warp3D_NN( imfl , nxnew,nynew,nznew , wf ) ;
00443        mri_free(imfl) ;
00444        imfl = mri_to_mri(im->kind,newImg) ;
00445        if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
00446        RETURN(newImg) ;
00447      }
00448 
00449      case MRI_rgb:{
00450        MRI_IMARR *imar = mri_rgb_to_3float(im) ;
00451        MRI_IMAGE *rim,*gim,*bim ;
00452        rim = mri_warp3D_NN( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00453                   mri_free( IMARR_SUBIM(imar,0) ) ;
00454        gim = mri_warp3D_NN( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00455                   mri_free( IMARR_SUBIM(imar,1) ) ;
00456        bim = mri_warp3D_NN( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
00457                   mri_free( IMARR_SUBIM(imar,2) ) ;
00458        FREE_IMARR(imar) ;
00459        newImg = mri_3to_rgb( rim,gim,bim ) ;
00460        mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
00461      }
00462 
00463      case MRI_complex:{
00464        MRI_IMARR *imar = mri_complex_to_pair(im) ;
00465        MRI_IMAGE *rim, *iim ;
00466        rim = mri_warp3D_NN( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00467                   mri_free( IMARR_SUBIM(imar,0) ) ;
00468        iim = mri_warp3D_NN( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00469                   mri_free( IMARR_SUBIM(imar,1) ) ;
00470        FREE_IMARR(imar) ;
00471        newImg = mri_pair_to_complex( rim , iim ) ;
00472        mri_free(rim); mri_free(iim); RETURN(newImg);
00473      }
00474 
00475    } /* end of special cases of input datum */
00476 
00477    /*----- at this point, imfl is in float format -----*/
00478 
00479    far = MRI_FLOAT_PTR( imfl ) ;                         /* input image data */
00480 
00481    newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ;  /* make output image */
00482    nar = MRI_FLOAT_PTR( newImg ) ;                         /* output image data */
00483 
00484    /*** loop over output points and warp to them ***/
00485 
00486    nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
00487    for( kk=0 ; kk < nznew ; kk++ ){
00488     zpr = kk ;
00489     for( jj=0 ; jj < nynew ; jj++ ){
00490      ypr = jj ;
00491      for( ii=0 ; ii < nxnew ; ii++ ){
00492        xpr = ii ;
00493        if( SKIP(ii,jj,kk) ) continue ;    /* 19 Nov 2004 */
00494        wf( xpr,ypr,zpr , &xx,&yy,&zz ) ;  /* get xx,yy,zz in original image */
00495 
00496        if( zout &&
00497            (xx < -0.5 || xx > nxh ||
00498             yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){
00499 
00500           NAR(ii,jj,kk) = 0.0 ; continue ;
00501        }
00502 
00503        ix = rint(xx) ; jy = rint(yy) ; kz = rint(zz) ;
00504        CLIP(ix,nx1)  ; CLIP(jy,ny1)  ; CLIP(kz,nz1)  ;
00505 
00506        NAR(ii,jj,kk) = FAR(ix,jy,kz) ;
00507      }
00508     }
00509    }
00510 
00511    /*** cleanup and return ***/
00512 
00513    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00514    RETURN(newImg);
00515 }
00516 
00517 /* define quintic interpolation polynomials (Lagrange) */
00518 
00519 #define Q_M2(x)  (x*(x*x-1.0)*(2.0-x)*(x-3.0)*0.008333333)
00520 #define Q_M1(x)  (x*(x*x-4.0)*(x-1.0)*(x-3.0)*0.041666667)
00521 #define Q_00(x)  ((x*x-4.0)*(x*x-1.0)*(3.0-x)*0.083333333)
00522 #define Q_P1(x)  (x*(x*x-4.0)*(x+1.0)*(x-3.0)*0.083333333)
00523 #define Q_P2(x)  (x*(x*x-1.0)*(x+2.0)*(3.0-x)*0.041666667)
00524 #define Q_P3(x)  (x*(x*x-1.0)*(x*x-4.0)*0.008333333)
00525 
00526 /*---------------------------------------------------------------------------*/
00527 /*! Transform a 3D image geometrically, using quintic interpolation.
00528     See comments for mri_warp3D() for details about the parameters.
00529      - RWCox - 06 Aug 2003
00530 -----------------------------------------------------------------------------*/
00531 
00532 INLINE MRI_IMAGE *mri_warp3D_quintic(
00533                     MRI_IMAGE *im, int nxnew, int nynew, int nznew,
00534                     void wf(float,float,float,float *,float *,float *) )
00535 {
00536    MRI_IMAGE *imfl , *newImg ;
00537    float *far , *nar ;
00538    float xpr,ypr,zpr, xx,yy,zz, fx,fy,fz ;
00539    int ii,jj,kk, nx,ny,nz,nxy, nx1,ny1,nz1, ix,jy,kz, nxynew ;
00540    float bot,top,val ;
00541    float nxh,nyh,nzh ;
00542 
00543    int ix_m2,ix_m1,ix_00,ix_p1,ix_p2,ix_p3 ; /* interpolation indices */
00544    int jy_m2,jy_m1,jy_00,jy_p1,jy_p2,jy_p3 ; /* (input image) */
00545    int kz_m2,kz_m1,kz_00,kz_p1,kz_p2,kz_p3 ;
00546 
00547    float wt_m2,wt_m1,wt_00,wt_p1,wt_p2,wt_p3 ; /* interpolation weights */
00548 
00549    float f_jm2_km2, f_jm1_km2, f_j00_km2, f_jp1_km2, f_jp2_km2, f_jp3_km2,
00550          f_jm2_km1, f_jm1_km1, f_j00_km1, f_jp1_km1, f_jp2_km1, f_jp3_km1,
00551          f_jm2_k00, f_jm1_k00, f_j00_k00, f_jp1_k00, f_jp2_k00, f_jp3_k00,
00552          f_jm2_kp1, f_jm1_kp1, f_j00_kp1, f_jp1_kp1, f_jp2_kp1, f_jp3_kp1,
00553          f_jm2_kp2, f_jm1_kp2, f_j00_kp2, f_jp1_kp2, f_jp2_kp2, f_jp3_kp2,
00554          f_jm2_kp3, f_jm1_kp3, f_j00_kp3, f_jp1_kp3, f_jp2_kp3, f_jp3_kp3,
00555          f_km2    , f_km1    , f_k00    , f_kp1    , f_kp2    , f_kp3     ;
00556 
00557 ENTRY("mri_warp3D_quinitc") ;
00558 
00559    /*--- sanity check inputs ---*/
00560 
00561    if( im == NULL || wf == NULL ) RETURN(NULL) ;
00562 
00563    /*-- dimensional analysis --*/
00564 
00565    nx = im->nx ; nx1 = nx-1 ;          /* input image dimensions */
00566    ny = im->ny ; ny1 = ny-1 ;
00567    nz = im->nz ; nz1 = nz-1 ; nxy = nx*ny ;
00568 
00569    nxnew = (nxnew > 0) ? nxnew : nx ;  /* output image dimensions */
00570    nynew = (nynew > 0) ? nynew : ny ;
00571    nznew = (nznew > 0) ? nznew : nz ;
00572    nxynew = nxnew*nynew ;
00573 
00574    /*----- Allow for different input image types, by breaking them into
00575            components, doing each one separately, and then reassembling.
00576            This is needed since we only warp float-valued images below.  -----*/
00577 
00578    switch( im->kind ){
00579 
00580      case MRI_float:                              /* use input directly */
00581        imfl = im ; break ;
00582 
00583      default:{                                    /* floatize input */
00584        imfl = mri_to_float(im) ;
00585        newImg  = mri_warp3D_quintic( imfl , nxnew,nynew,nznew , wf ) ;
00586        mri_free(imfl) ;
00587        imfl = mri_to_mri(im->kind,newImg) ;
00588        if( imfl != NULL ){ mri_free(newImg); newImg = imfl; }
00589        RETURN(newImg) ;
00590      }
00591 
00592      case MRI_rgb:{                                /* break into 3 pieces */
00593        MRI_IMARR *imar = mri_rgb_to_3float(im) ;
00594        MRI_IMAGE *rim,*gim,*bim ;
00595        rim = mri_warp3D_quintic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00596                      mri_free( IMARR_SUBIM(imar,0) ) ;
00597        gim = mri_warp3D_quintic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00598                      mri_free( IMARR_SUBIM(imar,1) ) ;
00599        bim = mri_warp3D_quintic( IMARR_SUBIM(imar,2), nxnew,nynew,nznew, wf ) ;
00600                      mri_free( IMARR_SUBIM(imar,2) ) ;
00601        FREE_IMARR(imar) ;
00602        newImg = mri_3to_rgb( rim,gim,bim ) ;
00603        mri_free(rim); mri_free(gim); mri_free(bim); RETURN(newImg);
00604      }
00605 
00606      case MRI_complex:{                             /* break into 2 pieces */
00607        MRI_IMARR *imar = mri_complex_to_pair(im) ;
00608        MRI_IMAGE *rim, *iim ;
00609        rim = mri_warp3D_quintic( IMARR_SUBIM(imar,0), nxnew,nynew,nznew, wf ) ;
00610                      mri_free( IMARR_SUBIM(imar,0) ) ;
00611        iim = mri_warp3D_quintic( IMARR_SUBIM(imar,1), nxnew,nynew,nznew, wf ) ;
00612                      mri_free( IMARR_SUBIM(imar,1) ) ;
00613        FREE_IMARR(imar) ;
00614        newImg = mri_pair_to_complex( rim , iim ) ;
00615        mri_free(rim); mri_free(iim); RETURN(newImg);
00616      }
00617 
00618    } /* end of special cases of input datum */
00619 
00620    /*----- at this point, imfl is in float format -----*/
00621 
00622    far = MRI_FLOAT_PTR( imfl ) ;                         /* input image data */
00623 
00624    newImg = mri_new_vol( nxnew,nynew,nznew, MRI_float ) ;  /* make output image */
00625    nar = MRI_FLOAT_PTR( newImg ) ;                         /* output image data */
00626 
00627    bot = top = far[0] ;                             /* find input data range */
00628    for( ii=1 ; ii < imfl->nvox ; ii++ ){
00629           if( far[ii] > top ) top = far[ii] ;
00630      else if( far[ii] < bot ) bot = far[ii] ;
00631    }
00632 
00633    /*** loop over output points and warp to them ***/
00634 
00635    nxh = nx-0.5 ; nyh = ny-0.5 ; nzh = nz-0.5 ;
00636    for( kk=0 ; kk < nznew ; kk++ ){
00637     zpr = kk ;
00638     for( jj=0 ; jj < nynew ; jj++ ){
00639      ypr = jj ;
00640      for( ii=0 ; ii < nxnew ; ii++ ){
00641        xpr = ii ;
00642        if( SKIP(ii,jj,kk) ) continue ;    /* 19 Nov 2004 */
00643        wf( xpr,ypr,zpr , &xx,&yy,&zz ) ;  /* get xx,yy,zz in original image */
00644 
00645        if( zout &&
00646            (xx < -0.5 || xx > nxh ||
00647             yy < -0.5 || yy > nyh || zz < -0.5 || zz > nzh ) ){
00648 
00649           NAR(ii,jj,kk) = 0.0 ; continue ;
00650        }
00651 
00652        ix = floor(xx) ;  fx = xx - ix ;   /* integer and       */
00653        jy = floor(yy) ;  fy = yy - jy ;   /* fractional coords */
00654        kz = floor(zz) ;  fz = zz - kz ;   /* in input image    */
00655 
00656        /* compute indexes from which to interpolate (-2,-1,0,+1,+2,+3),
00657           but clipped to lie inside input image volume                 */
00658 
00659        ix_m1 = ix-1    ; ix_00 = ix      ; ix_p1 = ix+1    ; ix_p2 = ix+2    ;
00660        CLIP(ix_m1,nx1) ; CLIP(ix_00,nx1) ; CLIP(ix_p1,nx1) ; CLIP(ix_p2,nx1) ;
00661        ix_m2 = ix-2    ; ix_p3 = ix+3 ;
00662        CLIP(ix_m2,nx1) ; CLIP(ix_p3,nx1) ;
00663 
00664        jy_m1 = jy-1    ; jy_00 = jy      ; jy_p1 = jy+1    ; jy_p2 = jy+2    ;
00665        CLIP(jy_m1,ny1) ; CLIP(jy_00,ny1) ; CLIP(jy_p1,ny1) ; CLIP(jy_p2,ny1) ;
00666        jy_m2 = jy-2    ; jy_p3 = jy+3 ;
00667        CLIP(jy_m2,ny1) ; CLIP(jy_p3,ny1) ;
00668 
00669        kz_m1 = kz-1    ; kz_00 = kz      ; kz_p1 = kz+1    ; kz_p2 = kz+2    ;
00670        CLIP(kz_m1,nz1) ; CLIP(kz_00,nz1) ; CLIP(kz_p1,nz1) ; CLIP(kz_p2,nz1) ;
00671        kz_m2 = kz-2    ; kz_p3 = kz+3 ;
00672        CLIP(kz_m2,nz1) ; CLIP(kz_p3,nz1) ;
00673 
00674        wt_m1 = Q_M1(fx) ; wt_00 = Q_00(fx) ;  /* interpolation weights */
00675        wt_p1 = Q_P1(fx) ; wt_p2 = Q_P2(fx) ;  /* in x-direction        */
00676        wt_m2 = Q_M2(fx) ; wt_p3 = Q_P3(fx) ;
00677 
00678 #undef  XINT
00679 #define XINT(j,k) wt_m2*FAR(ix_m2,j,k)+wt_m1*FAR(ix_m1,j,k) \
00680                  +wt_00*FAR(ix_00,j,k)+wt_p1*FAR(ix_p1,j,k) \
00681                  +wt_p2*FAR(ix_p2,j,k)+wt_p3*FAR(ix_p3,j,k)
00682 
00683        /* interpolate to location ix+fx at each jy,kz level */
00684 
00685        f_jm2_km2 = XINT(jy_m2,kz_m2) ; f_jm1_km2 = XINT(jy_m1,kz_m2) ;
00686        f_j00_km2 = XINT(jy_00,kz_m2) ; f_jp1_km2 = XINT(jy_p1,kz_m2) ;
00687        f_jp2_km2 = XINT(jy_p2,kz_m2) ; f_jp3_km2 = XINT(jy_p3,kz_m2) ;
00688 
00689        f_jm2_km1 = XINT(jy_m2,kz_m1) ; f_jm1_km1 = XINT(jy_m1,kz_m1) ;
00690        f_j00_km1 = XINT(jy_00,kz_m1) ; f_jp1_km1 = XINT(jy_p1,kz_m1) ;
00691        f_jp2_km1 = XINT(jy_p2,kz_m1) ; f_jp3_km1 = XINT(jy_p3,kz_m1) ;
00692 
00693        f_jm2_k00 = XINT(jy_m2,kz_00) ; f_jm1_k00 = XINT(jy_m1,kz_00) ;
00694        f_j00_k00 = XINT(jy_00,kz_00) ; f_jp1_k00 = XINT(jy_p1,kz_00) ;
00695        f_jp2_k00 = XINT(jy_p2,kz_00) ; f_jp3_k00 = XINT(jy_p3,kz_00) ;
00696 
00697        f_jm2_kp1 = XINT(jy_m2,kz_p1) ; f_jm1_kp1 = XINT(jy_m1,kz_p1) ;
00698        f_j00_kp1 = XINT(jy_00,kz_p1) ; f_jp1_kp1 = XINT(jy_p1,kz_p1) ;
00699        f_jp2_kp1 = XINT(jy_p2,kz_p1) ; f_jp3_kp1 = XINT(jy_p3,kz_p1) ;
00700 
00701        f_jm2_kp2 = XINT(jy_m2,kz_p2) ; f_jm1_kp2 = XINT(jy_m1,kz_p2) ;
00702        f_j00_kp2 = XINT(jy_00,kz_p2) ; f_jp1_kp2 = XINT(jy_p1,kz_p2) ;
00703        f_jp2_kp2 = XINT(jy_p2,kz_p2) ; f_jp3_kp2 = XINT(jy_p3,kz_p2) ;
00704 
00705        f_jm2_kp3 = XINT(jy_m2,kz_p3) ; f_jm1_kp3 = XINT(jy_m1,kz_p3) ;
00706        f_j00_kp3 = XINT(jy_00,kz_p3) ; f_jp1_kp3 = XINT(jy_p1,kz_p3) ;
00707        f_jp2_kp3 = XINT(jy_p2,kz_p3) ; f_jp3_kp3 = XINT(jy_p3,kz_p3) ;
00708 
00709        /* interpolate to jy+fy at each kz level */
00710 
00711        wt_m1 = Q_M1(fy) ; wt_00 = Q_00(fy) ; wt_p1 = Q_P1(fy) ;
00712        wt_p2 = Q_P2(fy) ; wt_m2 = Q_M2(fy) ; wt_p3 = Q_P3(fy) ;
00713 
00714        f_km2 =  wt_m2 * f_jm2_km2 + wt_m1 * f_jm1_km2 + wt_00 * f_j00_km2
00715               + wt_p1 * f_jp1_km2 + wt_p2 * f_jp2_km2 + wt_p3 * f_jp3_km2 ;
00716 
00717        f_km1 =  wt_m2 * f_jm2_km1 + wt_m1 * f_jm1_km1 + wt_00 * f_j00_km1
00718               + wt_p1 * f_jp1_km1 + wt_p2 * f_jp2_km1 + wt_p3 * f_jp3_km1 ;
00719 
00720        f_k00 =  wt_m2 * f_jm2_k00 + wt_m1 * f_jm1_k00 + wt_00 * f_j00_k00
00721               + wt_p1 * f_jp1_k00 + wt_p2 * f_jp2_k00 + wt_p3 * f_jp3_k00 ;
00722 
00723        f_kp1 =  wt_m2 * f_jm2_kp1 + wt_m1 * f_jm1_kp1 + wt_00 * f_j00_kp1
00724               + wt_p1 * f_jp1_kp1 + wt_p2 * f_jp2_kp1 + wt_p3 * f_jp3_kp1 ;
00725 
00726        f_kp2 =  wt_m2 * f_jm2_kp2 + wt_m1 * f_jm1_kp2 + wt_00 * f_j00_kp2
00727               + wt_p1 * f_jp1_kp2 + wt_p2 * f_jp2_kp2 + wt_p3 * f_jp3_kp2 ;
00728 
00729        f_kp3 =  wt_m2 * f_jm2_kp3 + wt_m1 * f_jm1_kp3 + wt_00 * f_j00_kp3
00730               + wt_p1 * f_jp1_kp3 + wt_p2 * f_jp2_kp3 + wt_p3 * f_jp3_kp3 ;
00731 
00732        /* interpolate to kz+fz to get output */
00733 
00734        wt_m1 = Q_M1(fz) ; wt_00 = Q_00(fz) ; wt_p1 = Q_P1(fz) ;
00735        wt_p2 = Q_P2(fz) ; wt_m2 = Q_M2(fz) ; wt_p3 = Q_P3(fz) ;
00736 
00737        val =  wt_m2 * f_km2 + wt_m1 * f_km1 + wt_00 * f_k00
00738             + wt_p1 * f_kp1 + wt_p2 * f_kp2 + wt_p3 * f_kp3 ;
00739 
00740             if( val > top ) val = top ;   /* clip to input data range */
00741        else if( val < bot ) val = bot ;
00742 
00743        NAR(ii,jj,kk) = val ;
00744      }
00745     }
00746    }
00747 
00748    /*** cleanup and return ***/
00749 
00750    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00751    RETURN(newImg);
00752 }
00753 
00754 /*--------------------------------------------------------------------------*/
00755 /*! Generic warping function.  Calls the specific one for the desired
00756     interpolation method, which was set by mri_warp3D_method().
00757       - im = input image (3D).
00758       - nxnew,nynew,nznew = size of output image grid.
00759       - wf( inew,jnew,knew , iold,jold,kold ) is a function that takes
00760         as input 3 indices in [0..nxnew-1,0..nynew-1,0..nznew-1] (the output
00761         image index space) and returns indices in the input image space.
00762       - indices in and out of wf() are floats.
00763       - by default, values outside the original image volume are treated
00764         as zero; this mode can be turned off using mri_warp3D_zerout().
00765 ----------------------------------------------------------------------------*/
00766 
00767 INLINE MRI_IMAGE *mri_warp3D(
00768                     MRI_IMAGE *im, int nxnew, int nynew, int nznew,
00769                     void wf(float,float,float,float *,float *,float *) )
00770 {
00771    switch( wtype ){
00772 
00773      default:
00774      case MRI_CUBIC:
00775        return mri_warp3D_cubic  ( im , nxnew,nynew,nznew , wf ) ;
00776 
00777      case MRI_LINEAR:
00778        return mri_warp3D_linear ( im , nxnew,nynew,nznew , wf ) ;
00779 
00780      case MRI_NN:
00781        return mri_warp3D_NN     ( im , nxnew,nynew,nznew , wf ) ;
00782 
00783      case MRI_QUINTIC:  /* 06 Aug 2003 */
00784        return mri_warp3D_quintic( im , nxnew,nynew,nznew , wf ) ;
00785    }
00786    return NULL ;  /* unreachable */
00787 }
00788 
00789 /*---------------------------------------------------------------------------*/
00790 
00791 static float sx_scale, sy_scale, sz_scale ;  /* global scaler data */
00792 
00793 /*! Internal transform function for resize "warp". */
00794 
00795 static INLINE void w3dMRI_scaler( float  a,float  b,float  c,
00796                                   float *x,float *y,float *z )
00797 {
00798    *x = sx_scale * a; *y = sy_scale * b; *z = sz_scale * c;
00799 }
00800 
00801 /*--------------------------------------------------------------------------*/
00802 /*! Special image warp just for resizing.
00803 ----------------------------------------------------------------------------*/
00804 
00805 MRI_IMAGE *mri_warp3D_resize( MRI_IMAGE *im , int nxnew, int nynew, int nznew )
00806 {
00807    int nx,ny,nz , nnx,nny,nnz ;
00808 
00809    if( im == NULL ) return NULL ;
00810    nx  = im->nx ; ny  = im->ny ; nz  = im->nz ;
00811    nnx = nxnew  ; nny = nynew  ; nnz = nznew  ;
00812 
00813    if( nnx <= 0 && nny <= 0 && nnz <= 0 ) return NULL ;
00814 
00815    sx_scale = (nnx > 0) ? ((float)nx)/nnx : 0.0 ;  /* global variables */
00816    sy_scale = (nny > 0) ? ((float)ny)/nny : 0.0 ;
00817    sz_scale = (nnz > 0) ? ((float)nz)/nnz : 0.0 ;
00818 
00819    if( nnx <= 0 ){
00820      sx_scale = MAX(sy_scale,sz_scale) ;
00821      nnx      = rint(sx_scale*nx) ;
00822    }
00823    if( nny <= 0 ){
00824      sy_scale = MAX(sx_scale,sz_scale) ;
00825      nny      = rint(sy_scale*ny) ;
00826    }
00827    if( nnz <= 0 ){
00828      sz_scale = MAX(sx_scale,sy_scale) ;
00829      nnz      = rint(sz_scale*nz) ;
00830    }
00831 
00832    return mri_warp3D( im , nnx,nny,nnz , w3dMRI_scaler ) ;
00833 }
00834 
00835 /*---------------------------------------------------------------------------*/
00836 
00837 static float a11_aff,a12_aff,a13_aff,a14_aff ,
00838              a21_aff,a22_aff,a23_aff,a24_aff ,
00839              a31_aff,a32_aff,a33_aff,a34_aff  ;
00840 
00841 /*! Internal transform function for affine "warp". */
00842 
00843 static INLINE void w3dMRI_affine( float  a,float  b,float  c,
00844                                   float *x,float *y,float *z )
00845 {
00846    *x = a11_aff*a + a12_aff*b + a13_aff*c + a14_aff ;
00847    *y = a21_aff*a + a22_aff*b + a23_aff*c + a24_aff ;
00848    *z = a31_aff*a + a32_aff*b + a33_aff*c + a34_aff ;
00849 }
00850 
00851 /*--------------------------------------------------------------------------*/
00852 /*! Special image warp for affine transformation of indices.
00853 ----------------------------------------------------------------------------*/
00854 
00855 MRI_IMAGE * mri_warp3D_affine( MRI_IMAGE *im , THD_vecmat aff )
00856 {
00857    if( im == NULL ) return NULL ;
00858 
00859    a11_aff = aff.mm.mat[0][0] ; a12_aff = aff.mm.mat[0][1] ;
00860    a13_aff = aff.mm.mat[0][2] ; a14_aff = aff.vv.xyz[0]    ;
00861 
00862    a21_aff = aff.mm.mat[1][0] ; a22_aff = aff.mm.mat[1][1] ;
00863    a23_aff = aff.mm.mat[1][2] ; a24_aff = aff.vv.xyz[1]    ;
00864 
00865    a31_aff = aff.mm.mat[2][0] ; a32_aff = aff.mm.mat[2][1] ;
00866    a33_aff = aff.mm.mat[2][2] ; a34_aff = aff.vv.xyz[2]    ;
00867 
00868    return mri_warp3D( im , 0,0,0 , w3dMRI_scaler ) ;
00869 }
00870 
00871 /*===========================================================================*/
00872 /*===========================================================================*/
00873 /******** Functions for warping a dataset using coordinate transforms, *******/
00874 /******** rather than the index transformation of the functions above. *******/
00875 /*===========================================================================*/
00876 /*===========================================================================*/
00877 
00878 static THD_vecmat ijk_to_dicom_in, ijk_to_dicom_out, dicom_to_ijk_in ;
00879 
00880 static void (*warp_out_to_in)(float  ,float  ,float ,
00881                               float *,float *,float * ) ;
00882 
00883 /*----------------------------------------------------------------------------*/
00884 /*! This is the function passed to mri_warp3D() for the ijk_out to ijk_in
00885     transformation.  It does the ijk_out to dicom_out transform here,
00886     then calls the coordinate transform function for dicom_out to dicom_in,
00887     then does the dicom_in to ijk_in transform here.
00888 ------------------------------------------------------------------------------*/
00889 
00890 static INLINE void w3d_warp_func( float xout, float yout, float zout,
00891                                   float *xin, float *yin, float *zin )
00892 {
00893    THD_fvec3 xxx,yyy ; float xi,yi,zi ;
00894 
00895    LOAD_FVEC3(xxx,xout,yout,zout) ;
00896    yyy = VECMAT_VEC( ijk_to_dicom_out , xxx ) ;       /* ijk_out to dicom_out */
00897    warp_out_to_in( yyy.xyz[0],yyy.xyz[1],yyy.xyz[2], /* dicom_out to dicom_in */
00898                   &(xxx.xyz[0]) , &(xxx.xyz[1]), &(xxx.xyz[2]) ) ;
00899    yyy = VECMAT_VEC( dicom_to_ijk_in , xxx ) ;          /* dicom_in to ijk_in */
00900    *xin = yyy.xyz[0] ; *yin = yyy.xyz[1] ; *zin = yyy.xyz[2] ;
00901 }
00902 
00903 /*--------------------------------------------------------------------------*/
00904 /*! - Find the 8 corners of the input dataset (voxel edges, not centers).
00905     - Warp each one using the provided wfunc().
00906     - Return the min and max (x,y,z) coordinates of these warped points.
00907 ----------------------------------------------------------------------------*/
00908 
00909 static void warp_corners( THD_3dim_dataset *inset,
00910                           void wfunc(float,float,float,float *,float *,float *),
00911                           float *xb , float *xt ,
00912                           float *yb , float *yt , float *zb , float *zt )
00913 {
00914    THD_dataxes *daxes = inset->daxes ;
00915    THD_fvec3 corner , wcorn ;
00916    float nx0 = -0.5          , ny0 = -0.5          , nz0 = -0.5           ;
00917    float nx1 = daxes->nxx-0.5, ny1 = daxes->nyy-0.5, nz1 = daxes->nzz-0.5 ;
00918    float xx,yy,zz , xbot,ybot,zbot , xtop,ytop,ztop ;
00919 
00920    LOAD_FVEC3( corner , nx0,ny0,nz0 ) ;
00921    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00922    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00923    xbot = xtop = xx ;
00924    ybot = ytop = yy ;
00925    zbot = ztop = zz ;
00926 
00927    LOAD_FVEC3( corner , nx1,ny0,nz0 ) ;
00928    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00929    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00930    xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
00931    ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
00932    zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
00933 
00934    LOAD_FVEC3( corner , nx0,ny1,nz0 ) ;
00935    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00936    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00937    xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
00938    ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
00939    zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
00940 
00941    LOAD_FVEC3( corner , nx1,ny1,nz0 ) ;
00942    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00943    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00944    xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
00945    ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
00946    zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
00947 
00948    LOAD_FVEC3( corner , nx0,ny0,nz1 ) ;
00949    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00950    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00951    xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
00952    ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
00953    zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
00954 
00955    LOAD_FVEC3( corner , nx1,ny0,nz1 ) ;
00956    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00957    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00958    xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
00959    ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
00960    zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
00961 
00962    LOAD_FVEC3( corner , nx0,ny1,nz1 ) ;
00963    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00964    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00965    xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
00966    ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
00967    zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
00968 
00969    LOAD_FVEC3( corner , nx1,ny1,nz1 ) ;
00970    wcorn = VECMAT_VEC( ijk_to_dicom_in , corner ) ;
00971    wfunc( wcorn.xyz[0],wcorn.xyz[1],wcorn.xyz[2] , &xx,&yy,&zz ) ;
00972    xbot = MIN(xx,xbot); xtop = MAX(xx,xtop);
00973    ybot = MIN(yy,ybot); ytop = MAX(yy,ytop);
00974    zbot = MIN(zz,zbot); ztop = MAX(zz,ztop);
00975 
00976    *xb = xbot; *xt = xtop;
00977    *yb = ybot; *yt = ytop; *zb = zbot; *zt = ztop; return;
00978 }
00979 
00980 /*--------------------------------------------------------------------------*/
00981 /*! Geometrically transform a 3D dataset, producing a new dataset.
00982      - w_in2out transforms DICOM coords from input grid to output grid
00983        (only needed if newggg != NULL and flag == WARP3D_NEWGRID;
00984         otherwise, w_in2out can be NULL)
00985      - w_out2in is the inverse of w_in2out (cannot be NULL)
00986      - newggg is a pointer to something that determines the new
00987        grid on which the output dataset will be computed
00988         - NULL means compute on the inset grid
00989         - otherwise, flag determines how newggg is interpreted:
00990         - WARP3D_NEWGRID => newggg is a "float *" to a new grid size
00991         - WARP3D_NEWDSET => newggg is a "THD_3dim_dataset *" to a dataset
00992                             whose grid will be used for the output grid
00993      - prefix = new dataset prefix (if NULL, then "warped")
00994      - zpad   = number of planes to zeropad on each face of inset (>= 0)
00995      - flag   = see above
00996      - Interpolation method can be set using mri_warp3D_method().
00997      - At end, input dataset is unloaded from memory.
00998      - If input is 3D+time, the output dataset won't have slice-wise time
00999        offsets, even if the input did.
01000 ----------------------------------------------------------------------------*/
01001 
01002 INLINE THD_3dim_dataset * THD_warp3D(
01003                      THD_3dim_dataset *inset ,
01004                      void w_in2out(float,float,float,float *,float *,float *),
01005                      void w_out2in(float,float,float,float *,float *,float *),
01006                      void *newggg , char *prefix , int zpad , int flag        )
01007 {
01008    THD_3dim_dataset *outset , *qset , *newgset=NULL ;
01009    int nxin,nyin,nzin,nvals , ival ;
01010    int nxout,nyout,nzout ;
01011    float xbot,xtop , ybot,ytop , zbot,ztop ;
01012    int use_newgrid , use_oldgrid ;
01013    float ddd_newgrid , fac ;
01014    MRI_IMAGE *inim , *outim , *wim ;
01015 
01016 ENTRY("THD_warp3D") ;
01017 
01018    /*-- decide which grid the output will be computed on --*/
01019 
01020    use_newgrid = ((flag & WARP3D_GRIDMASK) == WARP3D_NEWGRID) && (newggg != NULL) ;
01021    if( use_newgrid ){
01022      float *gg = (float *)newggg ;
01023      if( thd_floatscan(1,gg) == 0 && *gg > 0.0 ) ddd_newgrid = *gg ;
01024      else                                        use_newgrid = 0 ;
01025    } else if( newggg != NULL && (flag & WARP3D_GRIDMASK) == WARP3D_NEWDSET ){
01026      newgset = (THD_3dim_dataset *) newggg ;
01027      if( !ISVALID_DSET(newgset) ) newgset = NULL ;
01028    }
01029    use_oldgrid = !use_newgrid && (newgset == NULL) ;
01030 
01031    /*-- see if inputs are valid --*/
01032 
01033    if( !ISVALID_DSET(inset)             ||
01034        w_out2in == NULL                 ||
01035        (w_in2out == NULL && use_newgrid)  ){
01036 
01037      fprintf(stderr,"** ERROR: THD_warp3D has bad inputs!\n") ;
01038      RETURN(NULL);
01039    }
01040 
01041    /*-- zeropad and replace input dataset, if desired --*/
01042 
01043    if( zpad > 0 ){
01044      qset = THD_zeropad( inset , zpad,zpad,zpad,zpad,zpad,zpad ,
01045                          "Quetzalcoatl" , ZPAD_PURGE ) ;
01046      DSET_unload(inset) ;
01047      if( qset == NULL ){
01048        fprintf(stderr,"** ERROR: THD_warp3D can't zeropad!\n"); RETURN(NULL);
01049      }
01050    } else {
01051      qset = inset ;
01052    }
01053 
01054    /*-- read input data from disk, if not already present --*/
01055 
01056    DSET_load(qset) ;
01057    if( !DSET_LOADED(qset) ){
01058      fprintf(stderr,"** ERROR: THD_warp3D can't load input dataset!\n") ;
01059      if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
01060      RETURN(NULL);
01061    }
01062 
01063    /*-- compute mapping from input dataset (i,j,k) to DICOM coords --*/
01064 
01065    { THD_vecmat ijk_to_xyz , xyz_to_dicom ;
01066 
01067      LOAD_DIAG_MAT( ijk_to_xyz.mm , qset->daxes->xxdel,
01068                                     qset->daxes->yydel, qset->daxes->zzdel );
01069      LOAD_FVEC3   ( ijk_to_xyz.vv , qset->daxes->xxorg,
01070                                     qset->daxes->yyorg, qset->daxes->zzorg );
01071 
01072      xyz_to_dicom.mm = qset->daxes->to_dicomm ;
01073      LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ;
01074 
01075      ijk_to_dicom_in = MUL_VECMAT( xyz_to_dicom , ijk_to_xyz ) ;
01076      dicom_to_ijk_in = INV_VECMAT( ijk_to_dicom_in ) ;
01077    }
01078 
01079    /*-- make empty output dataset --*/
01080 
01081    nxin  = DSET_NX(qset) ;
01082    nyin  = DSET_NY(qset) ;
01083    nzin  = DSET_NZ(qset) ;
01084    nvals = DSET_NVALS(qset) ;
01085 
01086    if( nxin*nyin*nzin <= 1 ){
01087      fprintf(stderr,"** ERROR: THD_warp3D has nxin=%d nyin=%d nzin=%d!\n",
01088              nxin,nyin,nzin ) ;
01089      if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
01090      RETURN(NULL) ;
01091    }
01092 
01093    if( use_oldgrid ){                /*-- output is on same grid as input --*/
01094 
01095      outset = EDIT_empty_copy( qset ) ;
01096 
01097    } else if( newgset != NULL ){     /*-- output on same grid as newgset --*/
01098 
01099      outset = EDIT_empty_copy( newgset ) ;
01100 
01101      EDIT_dset_items( outset ,
01102                         ADN_nvals       , nvals ,
01103                         ADN_type        , qset->type ,
01104                         ADN_view_type   , qset->view_type ,
01105                         ADN_func_type   , qset->func_type ,
01106                         ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
01107                       ADN_none ) ;
01108 
01109      if( DSET_NUM_TIMES(qset) > 1 )     /* and some time structure? */
01110        EDIT_dset_items( outset ,
01111                           ADN_ntt       , nvals ,
01112                           ADN_tunits    , DSET_TIMEUNITS(qset) ,
01113                           ADN_ttorg     , DSET_TIMEORIGIN(qset) ,
01114                           ADN_ttdel     , DSET_TR(qset) ,
01115                           ADN_ttdur     , DSET_TIMEDURATION(qset) ,
01116                         ADN_none ) ;
01117 
01118      THD_copy_datablock_auxdata( qset->dblk , outset->dblk ) ;  /* 08 Jun 2004 */
01119      INIT_STAT_AUX( qset , MAX_STAT_AUX , outset->stat_aux ) ;
01120 
01121    } else {                          /*-- output is on new grid --*/
01122 
01123      float xmid,ymid,zmid ;
01124      THD_ivec3 nxyz , orixyz ;
01125      THD_fvec3 dxyz , orgxyz ;
01126 
01127      /* compute DICOM coordinates of warped corners */
01128 
01129      warp_corners( qset, w_in2out, &xbot,&xtop, &ybot,&ytop, &zbot,&ztop ) ;
01130 
01131      nxout = (int)( (xtop-xbot)/ddd_newgrid+0.999 ); if( nxout < 1 ) nxout = 1;
01132      nyout = (int)( (ytop-ybot)/ddd_newgrid+0.999 ); if( nyout < 1 ) nyout = 1;
01133      nzout = (int)( (ztop-zbot)/ddd_newgrid+0.999 ); if( nzout < 1 ) nzout = 1;
01134 
01135      xmid = 0.5*(xbot+xtop); ymid = 0.5*(ybot+ytop); zmid = 0.5*(zbot+ztop);
01136      xbot = xmid-0.5*(nxout-1)*ddd_newgrid; xtop = xbot+(nxout-1)*ddd_newgrid;
01137      ybot = ymid-0.5*(nyout-1)*ddd_newgrid; ytop = ybot+(nyout-1)*ddd_newgrid;
01138      zbot = zmid-0.5*(nzout-1)*ddd_newgrid; ztop = zbot+(nzout-1)*ddd_newgrid;
01139 
01140 #if 0
01141      if( verb )
01142        fprintf(stderr,"++ Transformed grid:\n"
01143                       "++   xbot = %10.4g  xtop = %10.4g  nx = %d\n"
01144                       "++   ybot = %10.4g  ytop = %10.4g  ny = %d\n"
01145                       "++   zbot = %10.4g  ztop = %10.4g  nz = %d\n" ,
01146                xbot,xtop,nxout , ybot,ytop,nyout , zbot,ztop,nzout    ) ;
01147 #endif
01148 
01149      if( nxout*nyout*nzout <= 1 ){
01150        fprintf(stderr,"** ERROR: THD_warp3D has nxout=%d nyout=%d nzout=%d!\n",
01151                nxout,nyout,nzout ) ;
01152        if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
01153        RETURN(NULL) ;
01154      }
01155 
01156      nxyz.ijk[0] = nxout ; dxyz.xyz[0] = ddd_newgrid ;  /* setup axes */
01157      nxyz.ijk[1] = nyout ; dxyz.xyz[1] = ddd_newgrid ;
01158      nxyz.ijk[2] = nzout ; dxyz.xyz[2] = ddd_newgrid ;
01159 
01160      orixyz.ijk[0] = ORI_R2L_TYPE ; orgxyz.xyz[0] = xbot ;
01161      orixyz.ijk[1] = ORI_A2P_TYPE ; orgxyz.xyz[1] = ybot ;
01162      orixyz.ijk[2] = ORI_I2S_TYPE ; orgxyz.xyz[2] = zbot ;
01163 
01164      /** create dataset and mangle it into the desired shape **/
01165 
01166      outset = EDIT_empty_copy( NULL ) ;  /* void and formless dataset */
01167 
01168      EDIT_dset_items( outset ,           /* give it some structure! */
01169                         ADN_nxyz        , nxyz ,
01170                         ADN_xyzdel      , dxyz ,
01171                         ADN_xyzorg      , orgxyz ,
01172                         ADN_xyzorient   , orixyz ,
01173                         ADN_malloc_type , DATABLOCK_MEM_MALLOC ,
01174                         ADN_nvals       , nvals ,
01175                         ADN_type        , qset->type ,
01176                         ADN_view_type   , qset->view_type ,
01177                         ADN_func_type   , qset->func_type ,
01178                       ADN_none ) ;
01179 
01180      if( DSET_NUM_TIMES(qset) > 1 )     /* and some time structure? */
01181        EDIT_dset_items( outset ,
01182                           ADN_ntt       , nvals ,
01183                           ADN_tunits    , DSET_TIMEUNITS(qset) ,
01184                           ADN_ttorg     , DSET_TIMEORIGIN(qset) ,
01185                           ADN_ttdel     , DSET_TR(qset) ,
01186                           ADN_ttdur     , DSET_TIMEDURATION(qset) ,
01187                         ADN_none ) ;
01188 
01189      THD_copy_datablock_auxdata( qset->dblk , outset->dblk ) ;  /* 08 Jun 2004 */
01190      INIT_STAT_AUX( qset , MAX_STAT_AUX , outset->stat_aux ) ;
01191 
01192    } /*-- end of warping to new grid --*/
01193 
01194    nxout = DSET_NX(outset) ;
01195    nyout = DSET_NY(outset) ;
01196    nzout = DSET_NZ(outset) ;
01197 
01198    /*-- compute mapping from output dataset (i,j,k) to DICOM coords --*/
01199 
01200    { THD_vecmat ijk_to_xyz , xyz_to_dicom ;
01201 
01202      LOAD_DIAG_MAT( ijk_to_xyz.mm, outset->daxes->xxdel,
01203                                    outset->daxes->yydel, outset->daxes->zzdel );
01204      LOAD_FVEC3   ( ijk_to_xyz.vv, outset->daxes->xxorg,
01205                                    outset->daxes->yyorg, outset->daxes->zzorg );
01206 
01207      xyz_to_dicom.mm = outset->daxes->to_dicomm ;
01208      LOAD_FVEC3( xyz_to_dicom.vv , 0.0,0.0,0.0 ) ;
01209 
01210      ijk_to_dicom_out = MUL_VECMAT( xyz_to_dicom , ijk_to_xyz ) ;
01211    }
01212 
01213    /*-- add prefix to new dataset --*/
01214 
01215    if( !THD_filename_ok(prefix) ) prefix = "warped" ;
01216    EDIT_dset_items( outset , ADN_prefix,prefix , ADN_none ) ;
01217 
01218    /*-- loop over bricks and warp them --*/
01219 
01220    warp_out_to_in = w_out2in ;  /* for use in w3d_warp_func(), supra */
01221 
01222    for( ival=0 ; ival < nvals ; ival++ ){
01223      inim  = DSET_BRICK(qset,ival) ;
01224      fac   = DSET_BRICK_FACTOR(qset,ival) ;
01225      if( fac > 0.0 && fac != 1.0 ) wim = mri_scale_to_float( fac , inim ) ;
01226      else                          wim = inim ;
01227      outim = mri_warp3D( wim , nxout,nyout,nzout , w3d_warp_func ) ;
01228      if( outim == NULL ){
01229        fprintf(stderr,"** ERROR: THD_warp3D fails at ival=%d\n",ival);
01230        DSET_delete(outset);
01231        if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
01232        RETURN(NULL) ;
01233      }
01234 
01235      if( wim != inim ){  /* if we scaled input */
01236        float gtop , fimfac ;
01237        mri_free(wim) ;
01238 
01239        /* 20 Oct 2003: rescale if input was an integer type */
01240 
01241        if( outim->kind == MRI_float && MRI_IS_INT_TYPE(inim->kind) ){
01242          fimfac = fac ;
01243          wim = mri_scalize( outim , inim->kind , &fimfac ) ;
01244          EDIT_BRICK_FACTOR( outset , ival , fimfac ) ;
01245          mri_free(outim) ; outim = wim ;
01246        } else {                                     /* input not an integer: */
01247          EDIT_BRICK_FACTOR( outset , ival , 0.0 ) ; /* output=unscaled float */
01248        }
01249 
01250      }
01251      EDIT_substitute_brick( outset, ival,outim->kind,mri_data_pointer(outim) );
01252      DSET_unload_one( qset , ival ) ;
01253    }
01254 
01255    /*-- done!!! --*/
01256 
01257    if( qset != inset ) DSET_delete(qset) ; else DSET_unload(qset) ;
01258 
01259    THD_load_statistics( outset ) ;
01260    RETURN( outset ) ;
01261 }
01262 
01263 /*--------------------------------------------------------------------------*/
01264 
01265 static THD_vecmat aff_out2in , aff_in2out ; /* data for afo2i() and afi2o() */
01266 
01267 /*! Internal transform for THD_warp3D_affine(). ----------------------------*/
01268 
01269 static INLINE void afo2i( float  a,float  b,float  c,
01270                           float *x,float *y,float *z )
01271 {
01272    THD_fvec3 xxx , yyy ;
01273    LOAD_FVEC3( xxx , a,b,c ) ;
01274    yyy = VECMAT_VEC( aff_out2in , xxx ) ;
01275    UNLOAD_FVEC3( yyy , *x , *y , *z ) ;
01276 }
01277 
01278 /*! Internal transform for THD_warp3D_affine(). ----------------------------*/
01279 
01280 static INLINE void afi2o( float  a,float  b,float  c,
01281                           float *x,float *y,float *z )
01282 {
01283    THD_fvec3 xxx , yyy ;
01284    LOAD_FVEC3( xxx , a,b,c ) ;
01285    yyy = VECMAT_VEC( aff_in2out , xxx ) ;
01286    UNLOAD_FVEC3( yyy , *x , *y , *z ) ;
01287 }
01288 
01289 /*--------------------------------------------------------------------------*/
01290 /*! A version of THD_warp3D() that uses an affine map, specified in the
01291     out2in parameter.  The main advantage over THD_warp3D() is that the
01292     transformation functions internal -- afo2i() and afi2o() -- which
01293     makes it simpler on the caller and also somewhat faster.
01294 ----------------------------------------------------------------------------*/
01295 
01296 THD_3dim_dataset * THD_warp3D_affine(
01297                      THD_3dim_dataset *inset ,
01298                      THD_vecmat out2in ,
01299                      void *newggg , char *prefix , int zpad , int flag )
01300 {
01301    aff_out2in = out2in ;
01302    aff_in2out = INV_VECMAT(aff_out2in) ;
01303 
01304    return THD_warp3D( inset , afi2o,afo2i , newggg,prefix,zpad,flag ) ;
01305 }
01306 
01307 /*--------------------------------------------------------------------------*/
01308 /*! Internal transform functions for TT <-> MNI coords.  Both are RAI,
01309     as per AFNI internal logic.  11 Mar 2004 - RW Cox [Jury Duty Day]
01310     cf. http://www.mrc-cbu.cam.ac.uk/Imaging/Common/mnispace.shtml
01311 ----------------------------------------------------------------------------*/
01312 
01313 static INLINE void w3d_mni2tta( float mx , float my , float mz ,
01314                                 float *tx, float *ty, float *tz )
01315 {
01316    *tx = 0.99 * mx ;
01317 
01318    if( mz > 0.0 ){
01319      *ty =  0.9688 * my + 0.0460 * mz ;
01320      *tz = -0.0485 * my + 0.9189 * mz ;
01321    } else {
01322      *ty =  0.9688 * my + 0.0420 * mz ;
01323      *tz = -0.0485 * my + 0.8390 * mz ;
01324    }
01325 }
01326 
01327 /*-------  Should be the inverse of the above! -----------------------------*/
01328 
01329 static INLINE void w3d_tta2mni( float tx , float ty , float tz ,
01330                                 float *mx, float *my, float *mz )
01331 {
01332    *mx = 1.01010 * tx ;
01333    *my = 1.02962 * ty - 0.05154 * tz ;
01334    *mz = 0.05434 * ty + 1.08554 * tz ;
01335    if( *mz < 0.0 ) *mz *= 1.09523 ;
01336 }
01337 
01338 /*--------------------------------------------------------------------------*/
01339 
01340 THD_3dim_dataset * THD_warp3D_tta2mni( THD_3dim_dataset *inset , void *newggg ,
01341                                        char *prefix , int zpad , int flag )
01342 {
01343    return THD_warp3D( inset , w3d_tta2mni,w3d_mni2tta , NULL,prefix,0,0 ) ;
01344 }
01345 
01346 /*-------  Should be the inverse of the above! -----------------------------*/
01347 
01348 THD_3dim_dataset * THD_warp3D_mni2tta( THD_3dim_dataset *inset , void *newggg ,
01349                                        char *prefix , int zpad , int flag )
01350 {
01351    return THD_warp3D( inset , w3d_mni2tta,w3d_tta2mni , newggg,prefix,zpad,flag ) ;
01352 }