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 byte pairs 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 8192  /* good for Pentium processors */
00016 
00017 /*---------------------------------------------------------------------------
00018    Two-step interpolation and shifting
00019 -----------------------------------------------------------------------------*/
00020 
00021 typedef unsigned short ushort ;
00022 
00023 #undef  DTYPE
00024 #define DTYPE ushort
00025 #undef  DSIZE
00026 #define DSIZE 2       /* sizeof(DTYPE) */
00027 
00028 static int    nlcbuf = 0 ;     /* workspace */
00029 static DTYPE * lcbuf = NULL ;
00030 
00031 #define TSBOT 0.3  /* the "optimal" breakpoints */
00032 #define TSTOP 0.7
00033 
00034 static void ts_shift_2byte( int n , float af , DTYPE * f )
00035 {
00036    float aa ;
00037    int ibot,itop , ia ;
00038 
00039    if( fabs(af) < TSBOT ) return ;  /* don't do anything if shift is too small */
00040 
00041    af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;  /* ia = floor */
00042    aa = af - ia ;
00043 
00044    if( n > nlcbuf ){
00045       if( lcbuf != NULL ) free(lcbuf) ;
00046       lcbuf  = (DTYPE *) malloc( DSIZE * n ) ;
00047       nlcbuf = n ;
00048    }
00049 
00050    ibot = -ia  ;   if( ibot < 0   ) ibot = 0 ;
00051    itop = n-2-ia ; if( itop > n-1 ) itop = n-1 ;
00052 
00053 #if 1
00054    memset(lcbuf,0,n*DSIZE) ;   /* seems to be faster */
00055 #else
00056    memset(lcbuf,0,ibot*DSIZE) ;
00057    memset(lcbuf+(itop+1),0,(n-(itop+1))*DSIZE) ;
00058 #endif
00059 
00060    if( aa < TSBOT ){         /* NN to bottom */
00061 
00062       memcpy( lcbuf+ibot, f+(ibot+ia)  , (itop+1-ibot)*DSIZE ) ;
00063 
00064    } else if( aa > TSTOP ){  /* NN to top */
00065 
00066       memcpy( lcbuf+ibot, f+(ibot+1+ia), (itop+1-ibot)*DSIZE ) ;
00067 
00068    } else {                  /* average bottom and top */
00069       register byte *fb=(byte *)(f+(ibot+ia)) , lb=(byte *)(lcbuf+ibot) ;
00070       register int ii ;
00071 
00072       for( ii=ibot ; ii <= itop ; ii++ ){
00073          *lb = ( *fb + *(fb+2) ) >> 1 ; fb++ ; lb++ ;
00074          *lb = ( *fb + *(fb+2) ) >> 1 ; fb++ ; lb++ ;
00075       }
00076 
00077    }
00078    memcpy( f , lcbuf , DSIZE*n ) ;
00079    return ;
00080 }
00081 
00082 /*---------------------------------------------------------------------------
00083    Flip a 3D array about the (x,y) axes:
00084     i <--> nx-1-i    j <--> ny-1-j
00085 -----------------------------------------------------------------------------*/
00086 
00087 #define VV(i,j,k) v[(i)+(j)*nx+(k)*nxy]
00088 #define SX(i)     (nx1-(i))
00089 #define SY(j)     (ny1-(j))
00090 #define SZ(k)     (nz1-(k))
00091 
00092 static void flip_xy( int nx , int ny , int nz , DTYPE * v )
00093 {
00094    int ii,jj,kk ;
00095    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00096    DTYPE * r1 ;
00097 
00098    r1 = (DTYPE *) malloc(DSIZE*nx) ;  /* save 1 row */
00099 
00100    for( kk=0 ; kk < nz ; kk++ ){              /* for each slice */
00101       for( jj=0 ; jj < ny2 ; jj++ ){          /* first 1/2 of rows */
00102 
00103          /* swap rows jj and ny1-jj, flipping them in ii as well */
00104 
00105          for( ii=0 ; ii < nx ; ii++ ) r1[ii]           = VV(SX(ii),SY(jj),kk) ;
00106          for( ii=0 ; ii < nx ; ii++ ) VV(ii,SY(jj),kk) = VV(SX(ii),jj    ,kk) ;
00107          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj    ,kk) = r1[ii] ;
00108       }
00109       if( ny%2 == 1 ){                                             /* central row? */
00110          for( ii=0 ; ii < nx ; ii++ ) r1[ii]       = VV(SX(ii),jj,kk) ; /* flip it */
00111          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk) = r1[ii] ;           /* restore */
00112       }
00113    }
00114 
00115    free(r1) ; return ;
00116 }
00117 
00118 /*---------------------------------------------------------------------------
00119    Flip a 3D array about the (y,z) axes:
00120      j <--> ny-1-j   k <--> nz-1-k
00121 -----------------------------------------------------------------------------*/
00122 
00123 static void flip_yz( int nx , int ny , int nz , DTYPE * v )
00124 {
00125    int ii,jj,kk ;
00126    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00127    DTYPE * r1 ;
00128 
00129    r1 = (DTYPE *) malloc(DSIZE*ny) ;
00130 
00131    for( ii=0 ; ii < nx ; ii++ ){
00132       for( kk=0 ; kk < nz2 ; kk++ ){
00133          for( jj=0 ; jj < ny ; jj++ ) r1[jj]           = VV(ii,SY(jj),SZ(kk)) ;
00134          for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,SZ(kk)) = VV(ii,SY(jj),kk    ) ;
00135          for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,kk    ) = r1[jj] ;
00136       }
00137       if( nz%2 == 1 ){
00138          for( jj=0 ; jj < ny ; jj++ ) r1[jj]       = VV(ii,SY(jj),kk) ;
00139          for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,kk) = r1[jj] ;
00140       }
00141    }
00142 
00143    free(r1) ; return ;
00144 }
00145 
00146 /*---------------------------------------------------------------------------
00147    Flip a 3D array about the (x,z) axes:
00148      i <--> nx-1-i   k <--> nz-1-k
00149 -----------------------------------------------------------------------------*/
00150 
00151 static void flip_xz( int nx , int ny , int nz , DTYPE * v )
00152 {
00153    int ii,jj,kk ;
00154    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00155    DTYPE * r1 ;
00156 
00157    r1 = (DTYPE *) malloc(DSIZE*nx) ;
00158 
00159    for( jj=0 ; jj < ny ; jj++ ){
00160       for( kk=0 ; kk < nz2 ; kk++ ){
00161          for( ii=0 ; ii < nx ; ii++ ) r1[ii]           = VV(SX(ii),jj,SZ(kk)) ;
00162          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,SZ(kk)) = VV(SX(ii),jj,kk    ) ;
00163          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk    ) = r1[ii] ;
00164       }
00165       if( nz%2 == 1 ){
00166          for( ii=0 ; ii < nx ; ii++ ) r1[ii]       = VV(SX(ii),jj,kk) ;
00167          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk) = r1[ii] ;
00168       }
00169    }
00170 
00171    free(r1) ; return ;
00172 }
00173 
00174 /*---------------------------------------------------------------------------
00175    Apply an x-axis shear to a 3D array: x -> x + a*y + b*z + s
00176    (dilation factor "f" assumed to be 1.0)
00177 -----------------------------------------------------------------------------*/
00178 
00179 static void apply_xshear( float a , float b , float s ,
00180                           int nx , int ny , int nz , DTYPE * v )
00181 {
00182    DTYPE * fj0 , * fj1 ;
00183    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00184    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00185    int ii,jj,kk , nup,nst ;
00186    float a0 , a1 , st ;
00187 
00188    /* don't do anything if shift is too small */
00189 
00190    st = fabs(a)*ny2 + fabs(b)*nz2 + fabs(s) ; if( st < TSBOT ) return ;
00191 
00192    for( kk=0 ; kk < nz ; kk++ ){
00193       for( jj=0 ; jj < ny ; jj++ )
00194          ts_shift_2byte( nx, a*(jj-ny2) + b*(kk-nz2) + s, v + (jj*nx + kk*nxy) ) ;
00195    }
00196 
00197    return ;
00198 }
00199 
00200 /*---------------------------------------------------------------------------
00201    Apply a y-axis shear to a 3D array: y -> y + a*x + b*z + s
00202 -----------------------------------------------------------------------------*/
00203 
00204 static void apply_yshear( float a , float b , float s ,
00205                           int nx , int ny , int nz , DTYPE * v )
00206 {
00207    DTYPE * fj0 , * fj1 ;
00208    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00209    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00210    int ii,jj,kk , nup,nst ;
00211    float a0 , a1 , st ;
00212    int xnum , xx , xtop ;
00213 
00214    /* don't do anything if shift is too small */
00215 
00216    st = fabs(a)*nx2 + fabs(b)*nz2 + fabs(s) ; if( st < TSBOT ) return ;
00217 
00218    xnum = CACHE / (ny*DSIZE) ; if( xnum < 1 ) xnum = 1 ;
00219    fj0 = (DTYPE *) malloc( DSIZE * xnum*ny ) ;
00220 
00221    for( kk=0 ; kk < nz ; kk++ ){
00222       for( ii=0 ; ii < nx ; ii+=xnum ){
00223          xtop = MIN(nx-ii,xnum) ;
00224          for( jj=0; jj < ny; jj++ )
00225             for( xx=0 ; xx < xtop ; xx++ )
00226                fj0[jj+xx*ny] = VV(ii+xx,jj,kk) ;
00227          for( xx=0 ; xx < xtop ; xx++ )
00228             ts_shift_2byte( ny, a*(ii+xx-nx2) + b*(kk-nz2) + s, fj0+xx*ny ) ;
00229          for( jj=0; jj < ny; jj++ )
00230             for( xx=0 ; xx < xtop ; xx++ )
00231                VV(ii+xx,jj,kk) = fj0[jj+xx*ny] ;
00232       }
00233    }
00234 
00235    free(fj0) ; return ;
00236 }
00237 
00238 /*---------------------------------------------------------------------------
00239    Apply a z-axis shear to a 3D array: z -> z + a*x + b*y + s
00240 -----------------------------------------------------------------------------*/
00241 
00242 static void apply_zshear( float a , float b , float s ,
00243                           int nx , int ny , int nz , DTYPE * v )
00244 {
00245    DTYPE * fj0 , * fj1 ;
00246    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00247    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00248    int ii,jj,kk , nup,nst ;
00249    float a0 , a1 , st ;
00250    int xnum , xx , xtop ;
00251 
00252    /* don't do anything if shift is too small */
00253 
00254    st = fabs(a)*nx2 + fabs(b)*ny2 + fabs(s) ; if( st < TSBOT ) return ;
00255 
00256    xnum = CACHE / (nz*DSIZE) ; if( xnum < 1 ) xnum = 1 ;
00257    fj0 = (DTYPE *) malloc( DSIZE * xnum*nz ) ;
00258 
00259    for( jj=0 ; jj < ny ; jj++ ){
00260       for( ii=0 ; ii < nx ; ii+=xnum ){
00261          xtop = MIN(nx-ii,xnum) ;
00262          for( kk=0; kk < nz; kk++ )
00263             for( xx=0 ; xx < xtop ; xx++ )
00264                fj0[kk+xx*nz] = VV(ii+xx,jj,kk) ;
00265          for( xx=0 ; xx < xtop ; xx++ )
00266             ts_shift_2byte( nz, a*(ii+xx-nx2) + b*(jj-ny2) + s, fj0+xx*nz ) ;
00267          for( kk=0; kk < nz; kk++ )
00268             for( xx=0 ; xx < xtop ; xx++ )
00269                VV(ii+xx,jj,kk) = fj0[kk+xx*nz] ;
00270       }
00271    }
00272 
00273    free(fj0) ; return ;
00274 }
00275 
00276 /*---------------------------------------------------------------------------
00277    Apply a set of shears to a 3D array of byte pairs.
00278    Note that we assume that the dilation factors ("f") are all 1.
00279 -----------------------------------------------------------------------------*/
00280 
00281 static void apply_3shear( MCW_3shear shr ,
00282                           int   nx   , int   ny   , int   nz   , DTYPE * vol )
00283 {
00284    int qq ;
00285    float a , b , s ;
00286 
00287    if( ! ISVALID_3SHEAR(shr) ) return ;
00288 
00289    /* carry out a preliminary 180 flippo ? */
00290 
00291    if( shr.flip0 >= 0 ){
00292       switch( shr.flip0 + shr.flip1 ){
00293          case 1: flip_xy( nx,ny,nz,vol ) ; break ;
00294          case 2: flip_xz( nx,ny,nz,vol ) ; break ;
00295          case 3: flip_yz( nx,ny,nz,vol ) ; break ;
00296         default:                           return ;  /* should not occur */
00297       }
00298    }
00299 
00300    /* apply each shear */
00301 
00302    for( qq=0 ; qq < 4 ; qq++ ){
00303       switch( shr.ax[qq] ){
00304          case 0:
00305             a = shr.scl[qq][1] ;
00306             b = shr.scl[qq][2] ;
00307             s = shr.sft[qq]    ;
00308             apply_xshear( a,b,s , nx,ny,nz , vol ) ;
00309          break ;
00310 
00311          case 1:
00312             a = shr.scl[qq][0] ;
00313             b = shr.scl[qq][2] ;
00314             s = shr.sft[qq]    ;
00315             apply_yshear( a,b,s , nx,ny,nz , vol ) ;
00316          break ;
00317 
00318          case 2:
00319             a = shr.scl[qq][0] ;
00320             b = shr.scl[qq][1] ;
00321             s = shr.sft[qq]    ;
00322             apply_zshear( a,b,s , nx,ny,nz , vol ) ;
00323          break ;
00324       }
00325    }
00326 
00327    return ;
00328 }
00329 
00330 /*---------------------------------------------------------------------------
00331   Rotate and translate a 3D volume.
00332 -----------------------------------------------------------------------------*/
00333 
00334 void THD_rota_vol_2byte( int   nx   , int   ny   , int   nz   ,
00335                          float xdel , float ydel , float zdel , DTYPE * vol ,
00336                          int ax1,float th1, int ax2,float th2, int ax3,float th3,
00337                          int dcode , float dx , float dy , float dz )
00338 {
00339    MCW_3shear shr ;
00340 
00341    if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) return ;
00342 
00343    if( xdel == 0.0 ) xdel = 1.0 ;
00344    if( ydel == 0.0 ) ydel = 1.0 ;
00345    if( zdel == 0.0 ) zdel = 1.0 ;
00346 
00347    if( th1 == 0.0 && th2 == 0.0 && th3 == 0.0 ){  /* nudge rotation */
00348       th1 = 1.e-6 ; th2 = 1.1e-6 ; th3 = 0.9e-6 ;
00349    }
00350 
00351    shr = rot_to_shear( ax1,-th1 , ax2,-th2 , ax3,-th3 ,
00352                        dcode,dx,dy,dz , xdel,ydel,zdel ) ;
00353 
00354    if( ! ISVALID_3SHEAR(shr) ){
00355       fprintf(stderr,"*** THD_rota_vol_2byte: can't compute shear transformation!\n") ;
00356       return ;
00357    }
00358 
00359    /************************************/
00360 
00361    apply_3shear( shr , nx,ny,nz , vol ) ;
00362 
00363    /************************************/
00364 
00365    return ;
00366 }
00367 
00368 /****************************************************************************
00369   Alternative entries, with rotation specified via a 3x3 matrix
00370   and shift as a 3-vector -- RWCox - 16 July 2000
00371 *****************************************************************************/
00372 
00373 /*---------------------------------------------------------------------------
00374   Rotate and translate a 3D volume
00375 -----------------------------------------------------------------------------*/
00376 
00377 #undef CLIPIT
00378 
00379 void THD_rota_vol_matvec_2byte( int   nx   , int   ny   , int   nz   ,
00380                                 float xdel, float ydel, float zdel, DTYPE * vol,
00381                                 THD_mat33 rmat , THD_fvec3 tvec )
00382 {
00383    MCW_3shear shr ;
00384    int dcode ;
00385 
00386    if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) return ;
00387 
00388    if( xdel == 0.0 ) xdel = 1.0 ;
00389    if( ydel == 0.0 ) ydel = 1.0 ;
00390    if( zdel == 0.0 ) zdel = 1.0 ;
00391 
00392    shr = rot_to_shear_matvec( rmat , tvec , xdel,ydel,zdel ) ;
00393 
00394    if( ! ISVALID_3SHEAR(shr) ){
00395       fprintf(stderr,"*** THD_rota_vol_2byte: can't compute shear transformation!\n") ;
00396       return ;
00397    }
00398 
00399    /************************************/
00400 
00401    apply_3shear( shr , nx,ny,nz , vol ) ;
00402 
00403    /************************************/
00404 
00405    return ;
00406 }