00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006    
00007 
00008 /*===========================================================================
00009   Routines to rotate/shift a 3D volume of bytes using a 4 way shear
00010   decomposition and "two-step" interpolation -- RWCox - Oct 2000.
00011 =============================================================================*/
00012 
00013 #include "thd_shear3d.h"
00014 
00015 #define CACHE 7168 /* good for Pentium processors */
00016 
00017 #define TSBOT 0.3  /* the "optimal" breakpoints for ts_shift */
00018 #define TSTOP 0.7
00019 
00020 #define NNBOT 0.5
00021 
00022 /*---------------------------------------------------------------------------*/
00023 
00024 static int icode  = MRI_TSSHIFT ;
00025 static float sbot = TSBOT ;
00026 
00027 void THD_rota_byte_mode( int code )
00028 {
00029    if( code == MRI_NN ){
00030       icode = MRI_NN      ; sbot = NNBOT ;
00031    } else {
00032       icode = MRI_TSSHIFT ; sbot = TSBOT ;
00033    }
00034 }
00035 
00036 /*---------------------------------------------------------------------------*/
00037 
00038 #undef  DTYPE
00039 #define DTYPE byte  /* AKA unsigned char */
00040 #undef  DSIZE
00041 #define DSIZE 1     /* sizeof(DTYPE) */
00042 
00043 static int    nlcbuf = 0 ;     /* workspace */
00044 static DTYPE * lcbuf = NULL ;
00045 
00046 /*---------------------------------------------------------------------------*/
00047 
00048 static int nn_shift_byte( int n , float af , DTYPE * f )
00049 {
00050    register int ii , ia ;
00051    float aa ;
00052    int ibot,itop ;
00053 
00054    if( fabs(af) < NNBOT ) return 0 ; /* do nothing if shift is too small */
00055 
00056    for( ii=0 ; ii < n && f[ii] == 0 ; ii++ ) ; /* nada */
00057    if( ii == n ) return 0 ;          /* do nothing if data all zero */
00058 
00059    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00060    aa = af - ia ;
00061 
00062    if( n > nlcbuf ){
00063       if( lcbuf != NULL ) free(lcbuf) ;
00064       lcbuf  = (DTYPE *) malloc( DSIZE * n ) ;
00065       nlcbuf = n ;
00066    }
00067 
00068    ibot = -ia  ;   if( ibot < 0   ) ibot = 0 ;
00069    itop = n-2-ia ; if( itop > n-1 ) itop = n-1 ;
00070 
00071 #if 1
00072    memset(lcbuf,0,n*DSIZE) ;   /* seems to be faster */
00073 #else
00074    memset(lcbuf,0,ibot*DSIZE) ;
00075    memset(lcbuf+(itop+1),0,(n-(itop+1))*DSIZE) ;
00076 #endif
00077 
00078    if( aa < NNBOT ){         /* NN to bottom */
00079 
00080       memcpy( lcbuf+ibot, f+(ibot+ia)  , (itop+1-ibot)*DSIZE ) ;
00081 
00082    } else {                  /* NN to top */
00083 
00084       memcpy( lcbuf+ibot, f+(ibot+1+ia), (itop+1-ibot)*DSIZE ) ;
00085 
00086    }
00087    memcpy( f , lcbuf , DSIZE*n ) ;
00088    return 1 ;
00089 }
00090 
00091 /*---------------------------------------------------------------------------
00092    Two-step interpolation and shifting
00093 -----------------------------------------------------------------------------*/
00094 
00095 static int ts_shift_byte( int n , float af , DTYPE * f )
00096 {
00097    register int ii , ia , ix ;
00098    float aa ;
00099    int ibot,itop ;
00100 
00101    if( fabs(af) < TSBOT ) return 0 ; /* do nothing if shift is too small */
00102 
00103    for( ii=0 ; ii < n && f[ii] == 0 ; ii++ ) ; /* nada */
00104    if( ii == n ) return 0 ;          /* do nothing if data all zero */
00105 
00106    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00107    aa = af - ia ;
00108 
00109    if( n > nlcbuf ){
00110       if( lcbuf != NULL ) free(lcbuf) ;
00111       lcbuf  = (DTYPE *) malloc( DSIZE * n ) ;
00112       nlcbuf = n ;
00113    }
00114 
00115    ibot = -ia  ;   if( ibot < 0   ) ibot = 0 ;
00116    itop = n-2-ia ; if( itop > n-1 ) itop = n-1 ;
00117 
00118 #if 1
00119    memset(lcbuf,0,n*DSIZE) ;   /* seems to be faster */
00120 #else
00121    memset(lcbuf,0,ibot*DSIZE) ;
00122    memset(lcbuf+(itop+1),0,(n-(itop+1))*DSIZE) ;
00123 #endif
00124 
00125    if( aa < TSBOT ){         /* NN to bottom */
00126 
00127       memcpy( lcbuf+ibot, f+(ibot+ia)  , (itop+1-ibot)*DSIZE ) ;
00128 
00129    } else if( aa > TSTOP ){  /* NN to top */
00130 
00131       memcpy( lcbuf+ibot, f+(ibot+1+ia), (itop+1-ibot)*DSIZE ) ;
00132 
00133    } else {                  /* average bottom and top */
00134 
00135       for( ii=ibot ; ii <= itop ; ii++ ){
00136          ix = ii + ia ; lcbuf[ii] = ( f[ix] + f[ix+1] ) >> 1 ;
00137       }
00138 
00139    }
00140    memcpy( f , lcbuf , DSIZE*n ) ;
00141    return 1 ;
00142 }
00143 
00144 /*---------------------------------------------------------------------------
00145    Flip a 3D array about the (x,y) axes:
00146     i <--> nx-1-i    j <--> ny-1-j
00147 -----------------------------------------------------------------------------*/
00148 
00149 #define VV(i,j,k) v[(i)+(j)*nx+(k)*nxy]
00150 #define SX(i)     (nx1-(i))
00151 #define SY(j)     (ny1-(j))
00152 #define SZ(k)     (nz1-(k))
00153 
00154 static void flip_xy( int nx , int ny , int nz , DTYPE * v , Tmask * tm )
00155 {
00156    int ii,jj,kk ;
00157    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00158    DTYPE * r1 ;
00159 
00160    r1 = (DTYPE *) malloc(DSIZE*nx) ;  /* save 1 row */
00161 
00162    for( kk=0 ; kk < nz ; kk++ ){              /* for each slice */
00163       for( jj=0 ; jj < ny2 ; jj++ ){          /* first 1/2 of rows */
00164 
00165          /* swap rows jj and ny1-jj, flipping them in ii as well */
00166 
00167          if( TM_XLINE(tm,jj+kk*ny) || TM_XLINE(tm,SY(jj)+kk*ny) ){
00168             for( ii=0; ii < nx; ii++ ) r1[ii]           = VV(SX(ii),SY(jj),kk) ;
00169             for( ii=0; ii < nx; ii++ ) VV(ii,SY(jj),kk) = VV(SX(ii),jj    ,kk) ;
00170             for( ii=0; ii < nx; ii++ ) VV(ii,jj    ,kk) = r1[ii] ;
00171          }
00172       }
00173       if( ny%2 == 1 && TM_XLINE(tm,jj+kk*ny) ){  /* central row? */
00174          for( ii=0; ii < nx; ii++ ) r1[ii]       = VV(SX(ii),jj,kk); /* flip it */
00175          for( ii=0; ii < nx; ii++ ) VV(ii,jj,kk) = r1[ii] ;          /* restore */
00176       }
00177    }
00178 
00179    free(r1) ; return ;
00180 }
00181 
00182 /*---------------------------------------------------------------------------
00183    Flip a 3D array about the (y,z) axes:
00184      j <--> ny-1-j   k <--> nz-1-k
00185 -----------------------------------------------------------------------------*/
00186 
00187 static void flip_yz( int nx , int ny , int nz , DTYPE * v , Tmask * tm )
00188 {
00189    int ii,jj,kk ;
00190    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00191    DTYPE * r1 ;
00192 
00193    r1 = (DTYPE *) malloc(DSIZE*ny) ;
00194 
00195    for( ii=0 ; ii < nx ; ii++ ){
00196       for( kk=0 ; kk < nz2 ; kk++ ){
00197          if( TM_YLINE(tm,kk+ii*nz) || TM_YLINE(tm,SZ(kk)+ii*nz) ){
00198             for( jj=0; jj < ny; jj++ ) r1[jj]           = VV(ii,SY(jj),SZ(kk)) ;
00199             for( jj=0; jj < ny; jj++ ) VV(ii,jj,SZ(kk)) = VV(ii,SY(jj),kk    ) ;
00200             for( jj=0; jj < ny; jj++ ) VV(ii,jj,kk    ) = r1[jj] ;
00201          }
00202       }
00203       if( nz%2 == 1 && TM_YLINE(tm,kk+ii*nz) ){
00204          for( jj=0; jj < ny; jj++ ) r1[jj]       = VV(ii,SY(jj),kk) ;
00205          for( jj=0; jj < ny; jj++ ) VV(ii,jj,kk) = r1[jj] ;
00206       }
00207    }
00208 
00209    free(r1) ; return ;
00210 }
00211 
00212 /*---------------------------------------------------------------------------
00213    Flip a 3D array about the (x,z) axes:
00214      i <--> nx-1-i   k <--> nz-1-k
00215 -----------------------------------------------------------------------------*/
00216 
00217 static void flip_xz( int nx , int ny , int nz , DTYPE * v , Tmask * tm )
00218 {
00219    int ii,jj,kk ;
00220    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00221    DTYPE * r1 ;
00222 
00223    r1 = (DTYPE *) malloc(DSIZE*nx) ;
00224 
00225    for( jj=0 ; jj < ny ; jj++ ){
00226       for( kk=0 ; kk < nz2 ; kk++ ){
00227          if( TM_XLINE(tm,jj+kk*ny) || TM_XLINE(tm,jj+SZ(kk)*ny) ){
00228             for( ii=0; ii < nx; ii++ ) r1[ii]           = VV(SX(ii),jj,SZ(kk)) ;
00229             for( ii=0; ii < nx; ii++ ) VV(ii,jj,SZ(kk)) = VV(SX(ii),jj,kk    ) ;
00230             for( ii=0; ii < nx; ii++ ) VV(ii,jj,kk    ) = r1[ii] ;
00231          }
00232       }
00233       if( nz%2 == 1 && TM_XLINE(tm,jj+kk*ny) ){
00234          for( ii=0; ii < nx; ii++ ) r1[ii]       = VV(SX(ii),jj,kk) ;
00235          for( ii=0; ii < nx; ii++ ) VV(ii,jj,kk) = r1[ii] ;
00236       }
00237    }
00238 
00239    free(r1) ; return ;
00240 }
00241 
00242 /*---------------------------------------------------------------------------
00243    Apply an x-axis shear to a 3D array: x -> x + a*y + b*z + s
00244    (dilation factor "f" assumed to be 1.0)
00245 -----------------------------------------------------------------------------*/
00246 
00247 static void apply_xshear( float a , float b , float s ,
00248                           int nx , int ny , int nz , DTYPE * v , Tmask * tm )
00249 {
00250    DTYPE * fj0 ;
00251    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00252    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00253    int ii,jj,kk ;
00254    float st ;
00255 
00256    /* don't do anything if shift is too small */
00257 
00258    st = fabs(a)*ny2 + fabs(b)*nz2 + fabs(s); if( st < sbot ) return ;
00259 
00260    switch( icode ){
00261       default:
00262       case MRI_TSSHIFT:
00263         for( kk=0 ; kk < nz ; kk++ ){
00264           for( jj=0 ; jj < ny ; jj++ )
00265             if( TM_XLINE(tm,jj+kk*ny) )
00266               ts_shift_byte( nx, a*(jj-ny2)+b*(kk-nz2)+s, v+(jj*nx+kk*nxy) );
00267         }
00268       break ;
00269 
00270       case MRI_NN:
00271         for( kk=0 ; kk < nz ; kk++ ){
00272           for( jj=0 ; jj < ny ; jj++ )
00273             if( TM_XLINE(tm,jj+kk*ny) )
00274               nn_shift_byte( nx, a*(jj-ny2)+b*(kk-nz2)+s, v+(jj*nx+kk*nxy) );
00275         }
00276       break ;
00277    }
00278 
00279    return ;
00280 }
00281 
00282 /*---------------------------------------------------------------------------
00283    Apply a y-axis shear to a 3D array: y -> y + a*x + b*z + s
00284 -----------------------------------------------------------------------------*/
00285 
00286 static void apply_yshear( float a , float b , float s ,
00287                           int nx , int ny , int nz , DTYPE * v , Tmask * tm )
00288 {
00289    DTYPE * fj0 ;
00290    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00291    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00292    int ii,jj,kk ;
00293    float st ;
00294    int xnum , xx , xtop , *wk ;
00295 
00296    /* don't do anything if shift is too small */
00297 
00298    st = fabs(a)*nx2 + fabs(b)*nz2 + fabs(s) ; if( st < sbot ) return ;
00299 
00300    xnum = CACHE / (ny*DSIZE) ; if( xnum < 1 ) xnum = 1 ;
00301    fj0 = (DTYPE *) malloc( DSIZE * xnum*ny ) ;
00302    wk  = (int *)   malloc( sizeof(int)*xnum ) ;
00303 
00304    switch( icode ){
00305       default:
00306       case MRI_TSSHIFT:
00307         for( kk=0 ; kk < nz ; kk++ ){
00308           for( ii=0 ; ii < nx ; ii+=xnum ){
00309             xtop = MIN(nx-ii,xnum) ;
00310             for( xx=0 ; xx < xtop ; xx++ )
00311               wk[xx] = fabs(a*(ii+xx-nx2)+b*(kk-nz2)+s) > TSBOT
00312                        && TM_YLINE(tm,kk+(ii+xx)*nz) ;
00313             for( jj=0; jj < ny; jj++ )
00314               for( xx=0 ; xx < xtop ; xx++ )
00315                 if( wk[xx] ) fj0[jj+xx*ny] = VV(ii+xx,jj,kk) ;
00316             for( xx=0 ; xx < xtop ; xx++ )
00317               if( wk[xx] )
00318                wk[xx] = ts_shift_byte(ny, a*(ii+xx-nx2)+b*(kk-nz2)+s, fj0+xx*ny);
00319             for( jj=0; jj < ny; jj++ )
00320               for( xx=0 ; xx < xtop ; xx++ )
00321                 if( wk[xx] ) VV(ii+xx,jj,kk) = fj0[jj+xx*ny] ;
00322           }
00323         }
00324       break ;
00325 
00326       case MRI_NN:
00327         for( kk=0 ; kk < nz ; kk++ ){
00328           for( ii=0 ; ii < nx ; ii+=xnum ){
00329             xtop = MIN(nx-ii,xnum) ;
00330             for( xx=0 ; xx < xtop ; xx++ )
00331               wk[xx] = fabs(a*(ii+xx-nx2)+b*(kk-nz2)+s) > NNBOT
00332                        && TM_YLINE(tm,kk+(ii+xx)*nz) ;
00333             for( jj=0; jj < ny; jj++ )
00334               for( xx=0 ; xx < xtop ; xx++ )
00335                 if( wk[xx] ) fj0[jj+xx*ny] = VV(ii+xx,jj,kk) ;
00336             for( xx=0 ; xx < xtop ; xx++ )
00337               if( wk[xx] )
00338                wk[xx] = nn_shift_byte(ny, a*(ii+xx-nx2)+b*(kk-nz2)+s, fj0+xx*ny);
00339             for( jj=0; jj < ny; jj++ )
00340               for( xx=0 ; xx < xtop ; xx++ )
00341                 if( wk[xx] ) VV(ii+xx,jj,kk) = fj0[jj+xx*ny] ;
00342           }
00343         }
00344       break ;
00345    }
00346 
00347    free(wk) ; free(fj0) ; return ;
00348 }
00349 
00350 /*---------------------------------------------------------------------------
00351    Apply a z-axis shear to a 3D array: z -> z + a*x + b*y + s
00352 -----------------------------------------------------------------------------*/
00353 
00354 static void apply_zshear( float a , float b , float s ,
00355                           int nx , int ny , int nz , DTYPE * v , Tmask * tm )
00356 {
00357    DTYPE * fj0 ;
00358    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00359    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00360    int ii,jj,kk ;
00361    float st ;
00362    int xnum , xx , xtop , *wk ;
00363 
00364    /* don't do anything if shift is too small */
00365 
00366    st = fabs(a)*nx2 + fabs(b)*ny2 + fabs(s) ; if( st < sbot ) return ;
00367 
00368    xnum = CACHE / (nz*DSIZE) ; if( xnum < 1 ) xnum = 1 ;
00369    fj0 = (DTYPE *) malloc( DSIZE * xnum*nz ) ;
00370    wk  = (int *)   malloc( sizeof(int)*xnum ) ;
00371 
00372    switch( icode ){
00373       default:
00374       case MRI_TSSHIFT:
00375         for( jj=0 ; jj < ny ; jj++ ){
00376           for( ii=0 ; ii < nx ; ii+=xnum ){
00377             xtop = MIN(nx-ii,xnum) ;
00378             for( xx=0 ; xx < xtop ; xx++ )
00379                wk[xx] = fabs(a*(ii+xx-nx2)+b*(jj-ny2)+s) > TSBOT
00380                         && TM_ZLINE(tm,ii+jj*nx+xx) ;
00381             for( kk=0; kk < nz; kk++ )
00382               for( xx=0 ; xx < xtop ; xx++ )
00383                 if( wk[xx] ) fj0[kk+xx*nz] = VV(ii+xx,jj,kk) ;
00384             for( xx=0 ; xx < xtop ; xx++ )
00385               if( wk[xx] )
00386                wk[xx] = ts_shift_byte(nz, a*(ii+xx-nx2)+b*(jj-ny2)+s, fj0+xx*nz);
00387             for( kk=0; kk < nz; kk++ )
00388               for( xx=0 ; xx < xtop ; xx++ )
00389                 if( wk[xx] ) VV(ii+xx,jj,kk) = fj0[kk+xx*nz] ;
00390           }
00391         }
00392       break ;
00393 
00394       case MRI_NN:
00395         for( jj=0 ; jj < ny ; jj++ ){
00396           for( ii=0 ; ii < nx ; ii+=xnum ){
00397             xtop = MIN(nx-ii,xnum) ;
00398             for( xx=0 ; xx < xtop ; xx++ )
00399                wk[xx] = fabs(a*(ii+xx-nx2)+b*(jj-ny2)+s) > NNBOT
00400                         && TM_ZLINE(tm,ii+jj*nx+xx) ;
00401             for( kk=0; kk < nz; kk++ )
00402               for( xx=0 ; xx < xtop ; xx++ )
00403                 if( wk[xx] ) fj0[kk+xx*nz] = VV(ii+xx,jj,kk) ;
00404             for( xx=0 ; xx < xtop ; xx++ )
00405               if( wk[xx] )
00406                wk[xx] = nn_shift_byte(nz, a*(ii+xx-nx2)+b*(jj-ny2)+s, fj0+xx*nz);
00407             for( kk=0; kk < nz; kk++ )
00408               for( xx=0 ; xx < xtop ; xx++ )
00409                 if( wk[xx] ) VV(ii+xx,jj,kk) = fj0[kk+xx*nz] ;
00410           }
00411         }
00412       break ;
00413    }
00414 
00415    free(wk) ; free(fj0) ; return ;
00416 }
00417 
00418 /*---------------------------------------------------------------------------
00419    Apply a set of shears to a 3D array of bytes.
00420    Note that we assume that the dilation factors ("f") are all 1.
00421 -----------------------------------------------------------------------------*/
00422 
00423 static void apply_3shear( MCW_3shear shr ,
00424                           int nx, int ny, int nz, DTYPE * vol , Tmask * tm )
00425 {
00426    int qq ;
00427    float a , b , s ;
00428 
00429    if( ! ISVALID_3SHEAR(shr) ) return ;
00430 
00431    /* carry out a preliminary 180 flippo ? */
00432 
00433    if( shr.flip0 >= 0 ){
00434       switch( shr.flip0 + shr.flip1 ){
00435          case 1: flip_xy( nx,ny,nz,vol,tm ) ; break ;
00436          case 2: flip_xz( nx,ny,nz,vol,tm ) ; break ;
00437          case 3: flip_yz( nx,ny,nz,vol,tm ) ; break ;
00438       }
00439    }
00440 
00441    /* apply each shear */
00442 
00443    for( qq=0 ; qq < 4 ; qq++ ){
00444       switch( shr.ax[qq] ){
00445          case 0:
00446             a = shr.scl[qq][1] ;
00447             b = shr.scl[qq][2] ;
00448             s = shr.sft[qq]    ;
00449             apply_xshear( a,b,s , nx,ny,nz , vol , (qq==0)? tm : NULL ) ;
00450          break ;
00451 
00452          case 1:
00453             a = shr.scl[qq][0] ;
00454             b = shr.scl[qq][2] ;
00455             s = shr.sft[qq]    ;
00456             apply_yshear( a,b,s , nx,ny,nz , vol , (qq==0)? tm : NULL ) ;
00457          break ;
00458 
00459          case 2:
00460             a = shr.scl[qq][0] ;
00461             b = shr.scl[qq][1] ;
00462             s = shr.sft[qq]    ;
00463             apply_zshear( a,b,s , nx,ny,nz , vol , (qq==0)? tm : NULL ) ;
00464          break ;
00465       }
00466    }
00467 
00468    return ;
00469 }
00470 
00471 /*---------------------------------------------------------------------------
00472   Rotate and translate a 3D volume.
00473 -----------------------------------------------------------------------------*/
00474 
00475 void THD_rota_vol_byte( int   nx   , int   ny   , int   nz   ,
00476                         float xdel , float ydel , float zdel , DTYPE * vol ,
00477                         int ax1,float th1, int ax2,float th2, int ax3,float th3,
00478                         int dcode , float dx , float dy , float dz , Tmask * tm )
00479 {
00480    MCW_3shear shr ;
00481 
00482    if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) return ;
00483 
00484    if( xdel == 0.0 ) xdel = 1.0 ;
00485    if( ydel == 0.0 ) ydel = 1.0 ;
00486    if( zdel == 0.0 ) zdel = 1.0 ;
00487 
00488    if( th1 == 0.0 && th2 == 0.0 && th3 == 0.0 ){  /* nudge rotation */
00489       th1 = 1.e-6 ; th2 = 1.1e-6 ; th3 = 0.9e-6 ;
00490    }
00491 
00492    shr = rot_to_shear( ax1,-th1 , ax2,-th2 , ax3,-th3 ,
00493                        dcode,dx,dy,dz , xdel,ydel,zdel ) ;
00494 
00495    if( ! ISVALID_3SHEAR(shr) ){
00496       fprintf(stderr,"*** THD_rota_vol_byte: can't compute shear transformation!\n") ;
00497       return ;
00498    }
00499 
00500    /*****************************************/
00501 
00502    apply_3shear( shr , nx,ny,nz , vol , tm ) ;
00503 
00504    /*****************************************/
00505 
00506    return ;
00507 }
00508 
00509 #if 0
00510 /****************************************************************************
00511   Alternative entries, with rotation specified via a 3x3 matrix
00512   and shift as a 3-vector -- RWCox - 16 July 2000
00513 *****************************************************************************/
00514 
00515 /*---------------------------------------------------------------------------
00516   Rotate and translate a 3D volume
00517 -----------------------------------------------------------------------------*/
00518 
00519 #undef CLIPIT
00520 
00521 void THD_rota_vol_matvec_byte( int   nx   , int   ny   , int   nz   ,
00522                                float xdel , float ydel , float zdel , DTYPE * vol ,
00523                                THD_mat33 rmat , THD_fvec3 tvec , Tmask * tm )
00524 {
00525    MCW_3shear shr ;
00526    int dcode ;
00527 
00528    if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) return ;
00529 
00530    if( xdel == 0.0 ) xdel = 1.0 ;
00531    if( ydel == 0.0 ) ydel = 1.0 ;
00532    if( zdel == 0.0 ) zdel = 1.0 ;
00533 
00534    shr = rot_to_shear_matvec( rmat , tvec , xdel,ydel,zdel ) ;
00535 
00536    if( ! ISVALID_3SHEAR(shr) ){
00537       fprintf(stderr,"*** THD_rota_vol_byte: can't compute shear transformation!\n") ;
00538       return ;
00539    }
00540 
00541    /************************************/
00542 
00543    apply_3shear( shr , nx,ny,nz , vol , tm ) ;
00544 
00545    /************************************/
00546 
00547    return ;
00548 }
00549 #endif