00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013 #include "thd_shear3d.h"
00014
00015 #define CACHE 8192
00016
00017
00018
00019
00020
00021 typedef unsigned short ushort ;
00022
00023 #undef DTYPE
00024 #define DTYPE ushort
00025 #undef DSIZE
00026 #define DSIZE 2
00027
00028 static int nlcbuf = 0 ;
00029 static DTYPE * lcbuf = NULL ;
00030
00031 #define TSBOT 0.3
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 ;
00040
00041 af = -af ; ia = (int) af ; if( af < 0 ) ia-- ;
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) ;
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 ){
00061
00062 memcpy( lcbuf+ibot, f+(ibot+ia) , (itop+1-ibot)*DSIZE ) ;
00063
00064 } else if( aa > TSTOP ){
00065
00066 memcpy( lcbuf+ibot, f+(ibot+1+ia), (itop+1-ibot)*DSIZE ) ;
00067
00068 } else {
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
00084
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) ;
00099
00100 for( kk=0 ; kk < nz ; kk++ ){
00101 for( jj=0 ; jj < ny2 ; jj++ ){
00102
00103
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 ){
00110 for( ii=0 ; ii < nx ; ii++ ) r1[ii] = VV(SX(ii),jj,kk) ;
00111 for( ii=0 ; ii < nx ; ii++ ) VV(ii,jj,kk) = r1[ii] ;
00112 }
00113 }
00114
00115 free(r1) ; return ;
00116 }
00117
00118
00119
00120
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
00148
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
00176
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
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
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
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
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
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
00278
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
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 ;
00297 }
00298 }
00299
00300
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
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 ){
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
00370
00371
00372
00373
00374
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 }