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   Routines to rotate and shift a 3D volume using a 4 way shear decomposition,
00009   coupled with an FFT-based shifting of the rows -- RWCox [October 1998].
00010 =============================================================================*/
00011 
00012 #include "thd_shear3d.h"  /* 23 Oct 2000: moved shear funcs to thd_shear3d.c */
00013 
00014 /********** 28 Oct 1999: the shifting routines that were here   **********
00015  **********              have been removed to file thd_shift2.c **********/
00016 
00017 /*---------------------------------------------------------------------------
00018    Set the interpolation method for shifting:
00019    input is one of MRI_NN, MRI_LINEAR, MRI_CUBIC, or MRI_FOURIER.
00020 -----------------------------------------------------------------------------*/
00021 
00022 typedef void (*shift_func)(int,int,float,float *,float,float *) ;
00023 static  shift_func shifter      = fft_shift2 ;
00024 static  int        shift_method = MRI_FOURIER ;
00025 
00026 void THD_rota_method( int mode )
00027 {
00028    shift_method = mode ;
00029    switch( mode ){
00030       case MRI_NN:      shifter = nn_shift2    ; break ;
00031       case MRI_TSSHIFT: shifter = ts_shift2    ; break ;  /* Dec 1999 */
00032 
00033       case MRI_LINEAR:  shifter = lin_shift2   ; break ;
00034       case MRI_FOURIER: shifter = fft_shift2   ; break ;
00035       default:
00036       case MRI_CUBIC:   shifter = cub_shift2   ; break ;
00037       case MRI_QUINTIC: shifter = quint_shift2 ; break ;  /* Nov 1998 */
00038       case MRI_HEPTIC:  shifter = hept_shift2  ; break ;  /* Nov 1998 */
00039 
00040       case MRI_FOURIER_NOPAD: shifter = fft_shift2   ; break ;  /* 13 May 2003 */
00041    }
00042    return ;
00043 }
00044 
00045 /*---------------------------------------------------------------------------
00046    Flip a 3D array about the (x,y) axes:
00047     i <--> nx-1-i    j <--> ny-1-j
00048 -----------------------------------------------------------------------------*/
00049 
00050 #define VV(i,j,k) v[(i)+(j)*nx+(k)*nxy]
00051 #define SX(i)     (nx1-(i))
00052 #define SY(j)     (ny1-(j))
00053 #define SZ(k)     (nz1-(k))
00054 
00055 static void flip_xy( int nx , int ny , int nz , float * v )
00056 {
00057    int ii,jj,kk ;
00058    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00059    float * r1 ;
00060 
00061    r1 = (float *) malloc(sizeof(float)*nx) ;  /* save 1 row */
00062 
00063    for( kk=0 ; kk < nz ; kk++ ){              /* for each slice */
00064       for( jj=0 ; jj < ny2 ; jj++ ){          /* first 1/2 of rows */
00065 
00066          /* swap rows jj and ny1-jj, flipping them in ii as well */
00067 
00068          for( ii=0 ; ii < nx ; ii++ ) r1[ii]           = VV(SX(ii),SY(jj),kk) ;
00069          for( ii=0 ; ii < nx ; ii++ ) VV(ii,SY(jj),kk) = VV(SX(ii),jj    ,kk) ;
00070          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj    ,kk) = r1[ii] ;
00071       }
00072       if( ny%2 == 1 ){                                             /* central row? */
00073          for( ii=0 ; ii < nx ; ii++ ) r1[ii]       = VV(SX(ii),jj,kk) ; /* flip it */
00074          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk) = r1[ii] ;           /* restore */
00075       }
00076    }
00077 
00078    free(r1) ; return ;
00079 }
00080 
00081 /*---------------------------------------------------------------------------
00082    Flip a 3D array about the (y,z) axes:
00083      j <--> ny-1-j   k <--> nz-1-k
00084 -----------------------------------------------------------------------------*/
00085 
00086 static void flip_yz( int nx , int ny , int nz , float * v )
00087 {
00088    int ii,jj,kk ;
00089    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00090    float * r1 ;
00091 
00092    r1 = (float *) malloc(sizeof(float)*ny) ;
00093 
00094    for( ii=0 ; ii < nx ; ii++ ){
00095       for( kk=0 ; kk < nz2 ; kk++ ){
00096          for( jj=0 ; jj < ny ; jj++ ) r1[jj]           = VV(ii,SY(jj),SZ(kk)) ;
00097          for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,SZ(kk)) = VV(ii,SY(jj),kk    ) ;
00098          for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,kk    ) = r1[jj] ;
00099       }
00100       if( nz%2 == 1 ){
00101          for( jj=0 ; jj < ny ; jj++ ) r1[jj]       = VV(ii,SY(jj),kk) ;
00102          for( jj=0 ; jj < ny ; jj++ ) VV(ii,jj,kk) = r1[jj] ;
00103       }
00104    }
00105 
00106    free(r1) ; return ;
00107 }
00108 
00109 /*---------------------------------------------------------------------------
00110    Flip a 3D array about the (x,z) axes:
00111      i <--> nx-1-i   k <--> nz-1-k
00112 -----------------------------------------------------------------------------*/
00113 
00114 static void flip_xz( int nx , int ny , int nz , float * v )
00115 {
00116    int ii,jj,kk ;
00117    int nx1=nx-1,nx2=nx/2, ny1=ny-1,ny2=ny/2, nz1=nz-1,nz2=nz/2, nxy=nx*ny ;
00118    float * r1 ;
00119 
00120    r1 = (float *) malloc(sizeof(float)*nx) ;
00121 
00122    for( jj=0 ; jj < ny ; jj++ ){
00123       for( kk=0 ; kk < nz2 ; kk++ ){
00124          for( ii=0 ; ii < nx ; ii++ ) r1[ii]           = VV(SX(ii),jj,SZ(kk)) ;
00125          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,SZ(kk)) = VV(SX(ii),jj,kk    ) ;
00126          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk    ) = r1[ii] ;
00127       }
00128       if( nz%2 == 1 ){
00129          for( ii=0 ; ii < nx ; ii++ ) r1[ii]       = VV(SX(ii),jj,kk) ;
00130          for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk) = r1[ii] ;
00131       }
00132    }
00133 
00134    free(r1) ; return ;
00135 }
00136 
00137 /*---------------------------------------------------------------------------
00138    Apply an x-axis shear to a 3D array: x -> x + a*y + b*z + s
00139    (dilation factor "f" assumed to be 1.0)
00140 -----------------------------------------------------------------------------*/
00141 
00142 static void apply_xshear( float a , float b , float s ,
00143                           int nx , int ny , int nz , float * v )
00144 {
00145    float * fj0 , * fj1 ;
00146    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00147    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00148    int ii,jj,kk , nup,nst ;
00149    float a0 , a1 , st ;
00150 
00151 ENTRY("apply_xshear") ;
00152 
00153    /* don't do anything if shift is less than 0.001 pixel */
00154 
00155    st = fabs(a)*ny2 + fabs(b)*nz2 + fabs(s) ; if( st < 1.e-3 ) EXRETURN ;
00156 
00157    if( shift_method == MRI_FOURIER ){
00158       nst = nx + 0.5*st ;
00159       nup = csfft_nextup_one35(nst) ;
00160    } else if( shift_method == MRI_FOURIER_NOPAD ){
00161       nup = csfft_nextup_even(nx) ;
00162    }
00163 
00164    for( kk=0 ; kk < nz ; kk++ ){
00165       for( jj=0 ; jj < ny ; jj+=2 ){
00166          fj0 = v + (jj*nx + kk*nxy) ;
00167          fj1 = (jj < ny1) ? (fj0 + nx) : NULL ;   /* allow for odd ny */
00168          a0  = a*(jj-ny2) + b*(kk-nz2) + s ;
00169          a1  = a0 + a ;
00170          shifter( nx , nup , a0 , fj0 , a1 , fj1 ) ;
00171       }
00172    }
00173 
00174    EXRETURN ;
00175 }
00176 
00177 /*---------------------------------------------------------------------------
00178    Apply a y-axis shear to a 3D array: y -> y + a*x + b*z + s
00179 -----------------------------------------------------------------------------*/
00180 
00181 static void apply_yshear( float a , float b , float s ,
00182                           int nx , int ny , int nz , float * v )
00183 {
00184    float * fj0 , * fj1 ;
00185    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00186    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00187    int ii,jj,kk , nup,nst ;
00188    float a0 , a1 , st ;
00189 
00190 ENTRY("apply_yshear") ;
00191 
00192    /* don't do anything if shift is less than 0.001 pixel */
00193 
00194    st = fabs(a)*nx2 + fabs(b)*nz2 + fabs(s) ; if( st < 1.e-3 ) EXRETURN ;
00195 
00196    if( shift_method == MRI_FOURIER ){
00197       nst = ny + 0.5*st ;
00198       nup = csfft_nextup_one35(nst) ;
00199    } else if( shift_method == MRI_FOURIER_NOPAD ){
00200       nup = csfft_nextup_even(ny) ;
00201    }
00202 
00203    fj0 = (float *) malloc( sizeof(float) * 2*ny ) ; fj1 = fj0 + ny ;
00204 
00205    for( kk=0 ; kk < nz ; kk++ ){
00206       for( ii=0 ; ii < nx1 ; ii+=2 ){
00207          for( jj=0; jj < ny; jj++ ){ fj0[jj] = VV(ii,jj,kk) ; fj1[jj] = VV(ii+1,jj,kk) ; }
00208          a0 = a*(ii-nx2) + b*(kk-nz2) + s ;
00209          a1 = a0 + a ;
00210          shifter( ny , nup , a0 , fj0 , a1 , fj1 ) ;
00211          for( jj=0; jj < ny; jj++ ){ VV(ii,jj,kk) = fj0[jj] ; VV(ii+1,jj,kk) = fj1[jj] ; }
00212       }
00213 
00214       if( ii == nx1 ){                                       /* allow for odd nx */
00215          for( jj=0; jj < ny; jj++ ) fj0[jj] = VV(ii,jj,kk) ;
00216          a0 = a*(ii-nx2) + b*(kk-nz2) + s ;
00217          shifter( ny , nup , a0 , fj0 , a1 , NULL ) ;
00218          for( jj=0; jj < ny; jj++ ) VV(ii,jj,kk) = fj0[jj] ;
00219       }
00220    }
00221 
00222    free(fj0) ; EXRETURN ;
00223 }
00224 
00225 /*---------------------------------------------------------------------------
00226    Apply a z-axis shear to a 3D array: z -> z + a*x + b*y + s
00227 -----------------------------------------------------------------------------*/
00228 
00229 static void apply_zshear( float a , float b , float s ,
00230                           int nx , int ny , int nz , float * v )
00231 {
00232    float * fj0 , * fj1 ;
00233    int   nx1=nx-1    , ny1=ny-1    , nz1=nz-1    , nxy=nx*ny ;
00234    float nx2=0.5*nx1 , ny2=0.5*ny1 , nz2=0.5*nz1 ;
00235    int ii,jj,kk , nup,nst ;
00236    float a0 , a1 , st ;
00237 
00238 ENTRY("apply_zshear") ;
00239 
00240    /* don't do anything if shift is less than 0.001 pixel */
00241 
00242    st = fabs(a)*nx2 + fabs(b)*ny2 + fabs(s) ; if( st < 1.e-3 ) EXRETURN ;
00243 
00244    if( shift_method == MRI_FOURIER ){
00245       nst = nz + 0.5*st ;
00246       nup = csfft_nextup_one35(nst) ;
00247    } else if( shift_method == MRI_FOURIER_NOPAD ){
00248       nup = csfft_nextup_even(nz) ;
00249    }
00250 
00251    fj0 = (float *) malloc( sizeof(float) * 2*nz ) ; fj1 = fj0 + nz ;
00252 
00253    for( jj=0 ; jj < ny ; jj++ ){
00254       for( ii=0 ; ii < nx1 ; ii+=2 ){
00255          for( kk=0; kk < nz; kk++ ){ fj0[kk] = VV(ii,jj,kk) ; fj1[kk] = VV(ii+1,jj,kk) ; }
00256          a0 = a*(ii-nx2) + b*(jj-ny2) + s ;
00257          a1 = a0 + a ;
00258          shifter( nz , nup , a0 , fj0 , a1 , fj1 ) ;
00259          for( kk=0; kk < nz; kk++ ){ VV(ii,jj,kk) = fj0[kk] ; VV(ii+1,jj,kk) = fj1[kk] ; }
00260       }
00261 
00262       if( ii == nx1 ){                                       /* allow for odd nx */
00263          for( kk=0; kk < nz; kk++ ) fj0[kk] = VV(ii,jj,kk) ;
00264          a0 = a*(ii-nx2) + b*(jj-ny2) + s ;
00265          shifter( nz , nup , a0 , fj0 , a1 , NULL ) ;
00266          for( kk=0; kk < nz; kk++ ) VV(ii,jj,kk) = fj0[kk] ;
00267       }
00268    }
00269 
00270    free(fj0) ; EXRETURN ;
00271 }
00272 
00273 /*---------------------------------------------------------------------------
00274    Apply a set of shears to a 3D array of floats.
00275    Note that we assume that the dilation factors ("f") are all 1.
00276 -----------------------------------------------------------------------------*/
00277 
00278 static void apply_3shear( MCW_3shear shr ,
00279                           int   nx   , int   ny   , int   nz   , float * vol )
00280 {
00281    int qq ;
00282    float a , b , s ;
00283 
00284 ENTRY("apply_3shear") ;
00285 
00286    if( ! ISVALID_3SHEAR(shr) ) EXRETURN ;
00287 
00288    /* carry out a preliminary 180 flippo ? */
00289 
00290    if( shr.flip0 >= 0 ){
00291       switch( shr.flip0 + shr.flip1 ){
00292          case 1: flip_xy( nx,ny,nz,vol ) ; break ;
00293          case 2: flip_xz( nx,ny,nz,vol ) ; break ;
00294          case 3: flip_yz( nx,ny,nz,vol ) ; break ;
00295         default:                           EXRETURN ;  /* should not occur */
00296       }
00297    }
00298 
00299    /* apply each shear */
00300 
00301    for( qq=0 ; qq < 4 ; qq++ ){
00302       switch( shr.ax[qq] ){
00303          case 0:
00304             a = shr.scl[qq][1] ;
00305             b = shr.scl[qq][2] ;
00306             s = shr.sft[qq]    ;
00307             apply_xshear( a,b,s , nx,ny,nz , vol ) ;
00308          break ;
00309 
00310          case 1:
00311             a = shr.scl[qq][0] ;
00312             b = shr.scl[qq][2] ;
00313             s = shr.sft[qq]    ;
00314             apply_yshear( a,b,s , nx,ny,nz , vol ) ;
00315          break ;
00316 
00317          case 2:
00318             a = shr.scl[qq][0] ;
00319             b = shr.scl[qq][1] ;
00320             s = shr.sft[qq]    ;
00321             apply_zshear( a,b,s , nx,ny,nz , vol ) ;
00322          break ;
00323       }
00324    }
00325 
00326    EXRETURN ;
00327 }
00328 
00329 /*---------------------------------------------------------------------------
00330    Set zero padding size for rotations:
00331    padding is done before rotate, then stripped off afterwards.
00332    02 Feb 2001 -- RWCox
00333 -----------------------------------------------------------------------------*/
00334 
00335 static int rotpx=0 , rotpy=0 , rotpz = 0 ;
00336 static int rotpset=0 ;
00337 
00338 void THD_rota_setpad( int px , int py , int pz )
00339 {
00340    rotpx = (px > 0) ? px : 0 ;
00341    rotpy = (py > 0) ? py : 0 ;
00342    rotpz = (pz > 0) ? pz : 0 ;
00343    rotpset = 1 ;               /* 05 Feb 2001 */
00344    return ;
00345 }
00346 
00347 /*---------------------------------------------------------------------------*/
00348 
00349 void THD_rota_clearpad(void)   /* 05 Feb 2001 */
00350 {
00351    rotpx=rotpy=rotpz=0; rotpset=1; return;
00352 }
00353 
00354 static void THD_rota_envpad(void)
00355 {
00356    char * eee = my_getenv("AFNI_ROTA_ZPAD") ;
00357    int ppp ;
00358 
00359    if( rotpset ) return ;
00360    eee = my_getenv("AFNI_ROTA_ZPAD") ;
00361    if( eee != NULL ){
00362       ppp = strtol( eee , NULL , 10 ) ;
00363       THD_rota_setpad(ppp,ppp,ppp) ;
00364    }
00365    rotpset = 1 ; return ;
00366 }
00367 
00368 /*---------------------------------------------------------------------------
00369   Rotate and translate a 3D volume.
00370 -----------------------------------------------------------------------------*/
00371 
00372 #undef CLIPIT
00373 
00374 void THD_rota_vol( int   nx   , int   ny   , int   nz   ,
00375                    float xdel , float ydel , float zdel , float * vol ,
00376                    int ax1,float th1, int ax2,float th2, int ax3,float th3,
00377                    int dcode , float dx , float dy , float dz )
00378 {
00379    MCW_3shear shr ;
00380 
00381 #ifdef CLIPIT
00382    register float bot , top ;
00383    register int   nxyz=nx*ny*nz , ii ;
00384 #endif
00385 
00386 ENTRY("THD_rota_vol") ;
00387 
00388    if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) EXRETURN ;
00389 
00390    if( xdel == 0.0 ) xdel = 1.0 ;
00391    if( ydel == 0.0 ) ydel = 1.0 ;
00392    if( zdel == 0.0 ) zdel = 1.0 ;
00393 
00394    if( th1 == 0.0 && th2 == 0.0 && th3 == 0.0 ){  /* nudge rotation */
00395       th1 = 1.e-6 ; th2 = 1.1e-6 ; th3 = 0.9e-6 ;
00396    }
00397 
00398 #if 0
00399 fprintf(stderr,"THD_rota_vol:\n") ;
00400 fprintf(stderr,"  th1=%g  th2=%g  th3=%g\n",th1,th2,th3) ;
00401 fprintf(stderr,"  dx=%g  dy=%g  dz=%g\n",dx,dy,dz) ;
00402 fprintf(stderr,"  xdel=%g  ydel=%g  zdel=%g\n",xdel,ydel,zdel) ;
00403 #endif
00404 
00405    shr = rot_to_shear( ax1,-th1 , ax2,-th2 , ax3,-th3 ,
00406                        dcode,dx,dy,dz , xdel,ydel,zdel ) ;
00407 
00408    if( ! ISVALID_3SHEAR(shr) ){
00409       fprintf(stderr,"*** THD_rota_vol: can't compute shear transformation!\n") ;
00410       EXRETURN ;
00411    }
00412 
00413 #if 0
00414    if( MRILIB_verbose )
00415      DUMP_3SHEAR("Computed shear",shr) ;
00416 #endif
00417 
00418 #ifdef CLIPIT
00419    bot = top = vol[0] ;
00420    for( ii=1 ; ii < nxyz ; ii++ ){
00421            if( vol[ii] < bot ) bot = vol[ii] ;
00422       else if( vol[ii] > top ) top = vol[ii] ;
00423    }
00424    if( bot >= top ) EXRETURN ;
00425 #endif
00426 
00427    /********************************/
00428    /* 02 Feb 2001: include padding */
00429 
00430    { float * vvv , *www ;
00431      int nxp , nyp , nzp ;
00432 
00433      THD_rota_envpad() ;  /* 05 Feb 2001 */
00434 
00435      nxp=nx+2*rotpx ; nyp=ny+2*rotpy ; nzp=nz+2*rotpz ;
00436 
00437      if( rotpx > 0 && rotpy > 0 && rotpz > 0 )
00438         vvv = EDIT_volpad_even( rotpx,rotpy,rotpz , nx,ny,nz , MRI_float,vol ) ;
00439      else
00440         vvv = vol ;
00441 
00442      apply_3shear( shr , nxp,nyp,nzp , vvv ) ;  /*-- do the actual rotation! --*/
00443 
00444      if( vvv != vol ){
00445         www = EDIT_volpad_even( -rotpx,-rotpy,-rotpz , nxp,nyp,nzp , MRI_float,vvv ) ;
00446         free(vvv) ;
00447         memcpy( vol , www , sizeof(float)*nx*ny*nz ) ; free(www) ;
00448      }
00449    }
00450 
00451    /********************************/
00452 
00453 #ifdef CLIPIT
00454    for( ii=0 ; ii < nxyz ; ii++ ){
00455            if( vol[ii] < bot ) vol[ii] = bot ;
00456       else if( vol[ii] > top ) vol[ii] = top ;
00457    }
00458 #endif
00459 
00460    EXRETURN ;
00461 }
00462 
00463 /*---------------------------------------------------------------------------
00464    Like the above, but with geometrical information about the volume
00465    given from the image header
00466 -----------------------------------------------------------------------------*/
00467 
00468 MRI_IMAGE * THD_rota3D( MRI_IMAGE * im ,
00469                         int ax1,float th1, int ax2,float th2, int ax3,float th3,
00470                         int dcode , float dx , float dy , float dz )
00471 {
00472    MRI_IMAGE * jm ;
00473    float * jvol ;
00474 
00475    if( ! MRI_IS_3D(im) ){
00476       fprintf(stderr,"\n*** THD_rota3D: non-3D image input!\n") ;
00477       return NULL ;
00478    }
00479 
00480    jm = mri_new_vol( im->nx , im->ny , im->nz , MRI_float ) ;
00481    MRI_COPY_AUX(jm,im) ;
00482    jvol = MRI_FLOAT_PTR(jm) ;
00483 
00484    EDIT_coerce_type( im->nvox ,
00485                      im->kind , mri_data_pointer(im) , MRI_float , jvol ) ;
00486 
00487    THD_rota_vol(      im->nx ,       im->ny ,       im->nz  ,
00488                  fabs(im->dx) , fabs(im->dy) , fabs(im->dz) , jvol ,
00489                  ax1,th1 , ax2,th2 , ax3,th3 , dcode , dx,dy,dz ) ;
00490 
00491    return jm ;
00492 }
00493 
00494 /****************************************************************************
00495   Alternative entries, with rotation specified via a 3x3 matrix
00496   and shift as a 3-vector -- RWCox - 16 July 2000
00497 *****************************************************************************/
00498 
00499 /*---------------------------------------------------------------------------
00500   Rotate and translate a 3D volume
00501 -----------------------------------------------------------------------------*/
00502 
00503 #undef CLIPIT
00504 
00505 void THD_rota_vol_matvec( int   nx   , int   ny   , int   nz   ,
00506                           float xdel , float ydel , float zdel , float * vol ,
00507                           THD_dmat33 rmat , THD_dfvec3 tvec )
00508 {
00509    MCW_3shear shr ;
00510    int dcode ;
00511 
00512 #ifdef CLIPIT
00513    register float bot , top ;
00514    register int   nxyz=nx*ny*nz , ii ;
00515 #endif
00516 
00517    if( nx < 2 || ny < 2 || nz < 2 || vol == NULL ) return ;
00518 
00519    if( xdel == 0.0 ) xdel = 1.0 ;
00520    if( ydel == 0.0 ) ydel = 1.0 ;
00521    if( zdel == 0.0 ) zdel = 1.0 ;
00522 
00523    shr = rot_to_shear_matvec( rmat , tvec , xdel,ydel,zdel ) ;
00524 
00525    if( ! ISVALID_3SHEAR(shr) ){
00526       fprintf(stderr,"*** THD_rota_vol: can't compute shear transformation!\n") ;
00527       return ;
00528    }
00529 
00530 #if 0
00531    if( MRILIB_verbose )
00532      DUMP_3SHEAR("Computed shear",shr) ;
00533 #endif
00534 
00535 #ifdef CLIPIT
00536    bot = top = vol[0] ;
00537    for( ii=1 ; ii < nxyz ; ii++ ){
00538            if( vol[ii] < bot ) bot = vol[ii] ;
00539       else if( vol[ii] > top ) top = vol[ii] ;
00540    }
00541    if( bot >= top ) return ;
00542 #endif
00543 
00544    /********************************/
00545    /* 02 Feb 2001: include padding */
00546 
00547    { float * vvv , *www ;
00548      int nxp , nyp , nzp ;
00549 
00550      THD_rota_envpad() ;  /* 05 Feb 2001 */
00551 
00552      nxp=nx+2*rotpx ; nyp=ny+2*rotpy ; nzp=nz+2*rotpz ;
00553 
00554      if( rotpx > 0 && rotpy > 0 && rotpz > 0 )
00555         vvv = EDIT_volpad_even( rotpx,rotpy,rotpz , nx,ny,nz , MRI_float,vol ) ;
00556      else
00557         vvv = vol ;
00558 
00559      apply_3shear( shr , nxp,nyp,nzp , vvv ) ;  /*-- do the actual rotation! --*/
00560 
00561      if( vvv != vol ){
00562         www = EDIT_volpad_even( -rotpx,-rotpy,-rotpz , nxp,nyp,nzp , MRI_float,vvv ) ;
00563         free(vvv) ;
00564         memcpy( vol , www , sizeof(float)*nx*ny*nz ) ; free(www) ;
00565      }
00566    }
00567 
00568    /********************************/
00569 
00570 #ifdef CLIPIT
00571    for( ii=0 ; ii < nxyz ; ii++ ){
00572            if( vol[ii] < bot ) vol[ii] = bot ;
00573       else if( vol[ii] > top ) vol[ii] = top ;
00574    }
00575 #endif
00576 
00577    return ;
00578 }
00579 
00580 /*------------------------------------------------------------------------------
00581    14 Feb 2001:
00582    Like the above, but with geometrical information about the volume
00583    given from the image header
00584 --------------------------------------------------------------------------------*/
00585 
00586 MRI_IMAGE * THD_rota3D_matvec( MRI_IMAGE * im, THD_dmat33 rmat,THD_dfvec3 tvec )
00587 {
00588    MRI_IMAGE * jm ;
00589    float * jvol ;
00590 
00591    if( ! MRI_IS_3D(im) ){
00592       fprintf(stderr,"\n*** THD_rota3D_matvec: non-3D image input!\n") ;
00593       return NULL ;
00594    }
00595 
00596    jm = mri_new_vol( im->nx , im->ny , im->nz , MRI_float ) ;
00597    MRI_COPY_AUX(jm,im) ;
00598    jvol = MRI_FLOAT_PTR(jm) ;
00599 
00600    EDIT_coerce_type( im->nvox ,
00601                      im->kind , mri_data_pointer(im) , MRI_float , jvol ) ;
00602 
00603    THD_rota_vol_matvec(      im->nx ,       im->ny ,       im->nz  ,
00604                        fabs(im->dx) , fabs(im->dy) , fabs(im->dz) , jvol ,
00605                        rmat , tvec ) ;
00606    return jm ;
00607 }