00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 #define FINS(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00012 ? 0.0 : far[(i)+(j)*nx] )
00013
00014
00015
00016 #undef USE_CGRID
00017
00018 #ifdef USE_CGRID
00019 # define CGRID 8192
00020 static int p_first = 1 ;
00021 static float p_m1[CGRID+1], p_00[CGRID+1],
00022 p_p1[CGRID+1], p_p2[CGRID+1] ;
00023 #endif
00024
00025
00026
00027 #define P_M1(x) ((x)*(1.0-(x))*((x)-2.0))
00028 #define P_00(x) (3.0*((x)+1.0)*((x)-1.0)*((x)-2.0))
00029 #define P_P1(x) (3.0*(x)*((x)+1.0)*(2.0-(x)))
00030 #define P_P2(x) ((x)*((x)+1.0)*((x)-1.0))
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042 MRI_IMAGE *mri_rota( MRI_IMAGE *im, float aa, float bb, float phi )
00043 {
00044 float rot_dx , rot_dy , rot_cph , rot_sph , top,bot,val ;
00045 MRI_IMAGE *imfl , *newImg ;
00046 MRI_IMARR *impair ;
00047 float *far , *nar ;
00048 float xx,yy , fx,fy ;
00049 int ii,jj, nx,ny , ix,jy , ifx,jfy ;
00050 float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
00051
00052 #ifdef USE_CGRID
00053 if( p_first ){
00054 p_first = 0 ;
00055 xx = 1.0 / CGRID ;
00056 for( ii=0 ; ii <= CGRID ; ii++ ){
00057 yy = ii * xx ;
00058 p_m1[ii] = P_M1(yy) ;
00059 p_00[ii] = P_00(yy) ;
00060 p_p1[ii] = P_P1(yy) ;
00061 p_p2[ii] = P_P2(yy) ;
00062 }
00063 }
00064 #endif
00065
00066 if( im == NULL || ! MRI_IS_2D(im) ){
00067 fprintf(stderr,"*** mri_rota only works on 2D images!\n") ; EXIT(1) ;
00068 }
00069
00070
00071
00072 if( im->kind == MRI_complex ){
00073 MRI_IMARR *impair ;
00074 MRI_IMAGE * rim , * iim , * tim ;
00075 impair = mri_complex_to_pair( im ) ;
00076 if( impair == NULL ){
00077 fprintf(stderr,"*** mri_complex_to_pair fails in mri_rota!\n") ; EXIT(1) ;
00078 }
00079 rim = IMAGE_IN_IMARR(impair,0) ;
00080 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
00081 tim = mri_rota( rim , aa,bb,phi ) ; mri_free( rim ) ; rim = tim ;
00082 tim = mri_rota( iim , aa,bb,phi ) ; mri_free( iim ) ; iim = tim ;
00083 newImg = mri_pair_to_complex( rim , iim ) ;
00084 mri_free( rim ) ; mri_free( iim ) ;
00085 MRI_COPY_AUX(newImg,im) ;
00086 return newImg ;
00087 }
00088
00089
00090
00091 rot_cph = cos(phi) ; rot_sph = sin(phi) ;
00092
00093 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00094 -(0.5 * im->ny) * rot_sph ;
00095
00096 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00097 +(0.5 * im->ny) * (1.0-rot_cph) ;
00098
00099
00100
00101 nx = im->nx ;
00102 ny = im->ny ;
00103
00104 if( im->kind == MRI_float ) imfl = im ;
00105 else imfl = mri_to_float( im ) ;
00106
00107 far = MRI_FLOAT_PTR(imfl) ;
00108 newImg = mri_new( nx , nx , MRI_float ) ;
00109 nar = MRI_FLOAT_PTR(newImg) ;
00110
00111 bot = top = far[0] ;
00112 for( ii=0 ; ii < nx*ny ; ii++ )
00113 if( far[ii] < bot ) bot = far[ii] ;
00114 else if( far[ii] > top ) top = far[ii] ;
00115
00116
00117
00118 for( jj=0 ; jj < nx ; jj++ ){
00119 xx = rot_sph * jj + rot_dx - rot_cph ;
00120 yy = rot_cph * jj + rot_dy + rot_sph ;
00121 for( ii=0 ; ii < nx ; ii++ ){
00122
00123 xx += rot_cph ;
00124 yy -= rot_sph ;
00125
00126 ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ;
00127 jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
00128
00129 #ifdef USE_CGRID
00130 ifx = (xx-ix)*CGRID + 0.499 ;
00131 wt_m1 = p_m1[ifx] ; wt_00 = p_00[ifx] ;
00132 wt_p1 = p_p1[ifx] ; wt_p2 = p_p2[ifx] ;
00133 #else
00134 fx = xx-ix ;
00135 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ;
00136 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
00137 #endif
00138
00139 if( ix > 0 && ix < nx-2 && jy > 0 && jy < ny-2 ){
00140 float * fym1, *fy00 , *fyp1 , *fyp2 ;
00141
00142 fym1 = far + (ix-1 + (jy-1)*nx) ;
00143 fy00 = fym1 + nx ;
00144 fyp1 = fy00 + nx ;
00145 fyp2 = fyp1 + nx ;
00146
00147 f_jm1 = wt_m1 * fym1[0] + wt_00 * fym1[1]
00148 + wt_p1 * fym1[2] + wt_p2 * fym1[3] ;
00149
00150 f_j00 = wt_m1 * fy00[0] + wt_00 * fy00[1]
00151 + wt_p1 * fy00[2] + wt_p2 * fy00[3] ;
00152
00153 f_jp1 = wt_m1 * fyp1[0] + wt_00 * fyp1[1]
00154 + wt_p1 * fyp1[2] + wt_p2 * fyp1[3] ;
00155
00156 f_jp2 = wt_m1 * fyp2[0] + wt_00 * fyp2[1]
00157 + wt_p1 * fyp2[2] + wt_p2 * fyp2[3] ;
00158
00159 } else {
00160
00161 f_jm1 = wt_m1 * FINS(ix-1,jy-1)
00162 + wt_00 * FINS(ix ,jy-1)
00163 + wt_p1 * FINS(ix+1,jy-1)
00164 + wt_p2 * FINS(ix+2,jy-1) ;
00165
00166 f_j00 = wt_m1 * FINS(ix-1,jy)
00167 + wt_00 * FINS(ix ,jy)
00168 + wt_p1 * FINS(ix+1,jy)
00169 + wt_p2 * FINS(ix+2,jy) ;
00170
00171 f_jp1 = wt_m1 * FINS(ix-1,jy+1)
00172 + wt_00 * FINS(ix ,jy+1)
00173 + wt_p1 * FINS(ix+1,jy+1)
00174 + wt_p2 * FINS(ix+2,jy+1) ;
00175
00176 f_jp2 = wt_m1 * FINS(ix-1,jy+2)
00177 + wt_00 * FINS(ix ,jy+2)
00178 + wt_p1 * FINS(ix+1,jy+2)
00179 + wt_p2 * FINS(ix+2,jy+2) ;
00180 }
00181
00182 #define THIRTYSIX 2.7777778e-2
00183
00184 #ifdef USE_CGRID
00185 jfy = (yy-jy)*CGRID + 0.499 ;
00186 val = ( p_m1[jfy] * f_jm1 + p_00[jfy] * f_j00
00187 + p_p1[jfy] * f_jp1 + p_p2[jfy] * f_jp2 ) * THIRTYSIX ;
00188 #else
00189 fy = yy-jy ;
00190 val = ( P_M1(fy) * f_jm1 + P_00(fy) * f_j00
00191 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) * THIRTYSIX ;
00192 #endif
00193
00194 if( val < bot ) nar[ii+jj*nx] = bot ;
00195 else if( val > top ) nar[ii+jj*nx] = top ;
00196 else nar[ii+jj*nx] = val ;
00197
00198 }
00199 }
00200
00201
00202
00203 if( im != imfl ) mri_free(imfl) ;
00204 MRI_COPY_AUX(newImg,im) ;
00205 return newImg ;
00206 }
00207
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218 MRI_IMAGE *mri_rota_bilinear( MRI_IMAGE *im, float aa, float bb, float phi )
00219 {
00220 float rot_dx , rot_dy , rot_cph , rot_sph ;
00221 MRI_IMAGE *imfl , *newImg ;
00222 MRI_IMARR *impair ;
00223 float *far , *nar ;
00224 float xx,yy , fx,fy ;
00225 int ii,jj, nx,ny , ix,jy ;
00226 float f_j00,f_jp1 , wt_00,wt_p1 ;
00227
00228 if( im == NULL || ! MRI_IS_2D(im) ){
00229 fprintf(stderr,"*** mri_rota_bilinear only works on 2D images!\n") ; EXIT(1) ;
00230 }
00231
00232
00233
00234 if( im->kind == MRI_complex ){
00235 MRI_IMARR *impair ;
00236 MRI_IMAGE * rim , * iim , * tim ;
00237 impair = mri_complex_to_pair( im ) ;
00238 if( impair == NULL ){
00239 fprintf(stderr,"*** mri_complex_to_pair fails in mri_rota!\n") ; EXIT(1) ;
00240 }
00241 rim = IMAGE_IN_IMARR(impair,0) ;
00242 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
00243 tim = mri_rota_bilinear( rim , aa,bb,phi ) ; mri_free( rim ) ; rim = tim ;
00244 tim = mri_rota_bilinear( iim , aa,bb,phi ) ; mri_free( iim ) ; iim = tim ;
00245 newImg = mri_pair_to_complex( rim , iim ) ;
00246 mri_free( rim ) ; mri_free( iim ) ;
00247 MRI_COPY_AUX(newImg,im) ;
00248 return newImg ;
00249 }
00250
00251
00252
00253 rot_cph = cos(phi) ; rot_sph = sin(phi) ;
00254
00255 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00256 -(0.5 * im->ny) * rot_sph ;
00257
00258 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00259 +(0.5 * im->ny) * (1.0-rot_cph) ;
00260
00261
00262
00263 nx = im->nx ;
00264 ny = im->ny ;
00265
00266 if( im->kind == MRI_float ) imfl = im ;
00267 else imfl = mri_to_float( im ) ;
00268
00269 far = MRI_FLOAT_PTR(imfl) ;
00270 newImg = mri_new( nx , nx , MRI_float ) ;
00271 nar = MRI_FLOAT_PTR(newImg) ;
00272
00273
00274
00275 for( jj=0 ; jj < nx ; jj++ ){
00276 xx = rot_sph * jj + rot_dx - rot_cph ;
00277 yy = rot_cph * jj + rot_dy + rot_sph ;
00278 for( ii=0 ; ii < nx ; ii++ ){
00279
00280 xx += rot_cph ;
00281 yy -= rot_sph ;
00282
00283 ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ;
00284 jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
00285
00286 fx = xx-ix ; wt_00 = 1.0 - fx ; wt_p1 = fx ;
00287
00288 if( ix >= 0 && ix < nx-1 && jy >= 0 && jy < ny-1 ){
00289 float *fy00 , *fyp1 ;
00290
00291 fy00 = far + (ix + jy*nx) ; fyp1 = fy00 + nx ;
00292
00293 f_j00 = wt_00 * fy00[0] + wt_p1 * fy00[1] ;
00294 f_jp1 = wt_00 * fyp1[0] + wt_p1 * fyp1[1] ;
00295
00296 } else {
00297 f_j00 = wt_00 * FINS(ix,jy ) + wt_p1 * FINS(ix+1,jy ) ;
00298 f_jp1 = wt_00 * FINS(ix,jy+1) + wt_p1 * FINS(ix+1,jy+1) ;
00299 }
00300
00301 fy = yy-jy ; nar[ii+jj*nx] = (1.0-fy) * f_j00 + fy * f_jp1 ;
00302
00303 }
00304 }
00305
00306
00307
00308 if( im != imfl ) mri_free(imfl) ;
00309 MRI_COPY_AUX(newImg,im) ;
00310 return newImg ;
00311 }
00312
00313
00314
00315
00316
00317
00318
00319 void ft_shift2( int n, int nup, float af, float * f, float ag, float * g )
00320 {
00321 static int nupold=0 ;
00322 static complex * row=NULL , * cf=NULL , * cg=NULL ;
00323
00324 int ii , nby2=nup/2 , n21=nby2+1 ;
00325 complex fac , gac ;
00326 float sf , sg , dk ;
00327
00328
00329
00330 if( nup > nupold ){
00331 if( row != NULL ){ free(row) ; free(cf) ; free(cg) ; }
00332 row = (complex *) malloc( sizeof(complex) * nup ) ;
00333 cf = (complex *) malloc( sizeof(complex) * n21 ) ;
00334 cg = (complex *) malloc( sizeof(complex) * n21 ) ;
00335 nupold = nup ;
00336 }
00337
00338
00339
00340 for( ii=0 ; ii < n ; ii++ ){ row[ii].r = f[ii] ; row[ii].i = g[ii] ; }
00341 for( ; ii < nup ; ii++ ){ row[ii].r = row[ii].i = 0.0 ; }
00342
00343 csfft_cox( -1 , nup , row ) ;
00344
00345
00346
00347 cf[0].r = 2.0 * row[0].r ; cf[0].i = 0.0 ;
00348 cg[0].r = 2.0 * row[0].i ; cg[0].i = 0.0 ;
00349 for( ii=1 ; ii < nby2 ; ii++ ){
00350 cf[ii].r = row[ii].r + row[nup-ii].r ;
00351 cf[ii].i = row[ii].i - row[nup-ii].i ;
00352 cg[ii].r = row[ii].i + row[nup-ii].i ;
00353 cg[ii].i = -row[ii].r + row[nup-ii].r ;
00354 }
00355 cf[nby2].r = 2.0 * row[nby2].r ; cf[nby2].i = 0.0 ;
00356 cg[nby2].r = 2.0 * row[nby2].i ; cg[nby2].i = 0.0 ;
00357
00358
00359
00360 dk = (2.0*PI) / nup ;
00361 sf = -af * dk ; sg = -ag * dk ;
00362 for( ii=1 ; ii <= nby2 ; ii++ ){
00363 fac = CEXPIT(ii*sf) ; cf[ii] = CMULT( fac , cf[ii] ) ;
00364 gac = CEXPIT(ii*sg) ; cg[ii] = CMULT( gac , cg[ii] ) ;
00365 }
00366 cf[nby2].i = 0.0 ; cg[nby2].i = 0.0 ;
00367
00368
00369
00370 row[0].r = cf[0].r ; row[0].i = cg[0].r ;
00371 for( ii=1 ; ii < nby2 ; ii++ ){
00372 row[ii].r = cf[ii].r - cg[ii].i ;
00373 row[ii].i = cf[ii].i + cg[ii].r ;
00374 row[nup-ii].r = cf[ii].r + cg[ii].i ;
00375 row[nup-ii].i = -cf[ii].i + cg[ii].r ;
00376 }
00377 row[nby2].r = cf[nby2].r ;
00378 row[nby2].i = cg[nby2].r ;
00379
00380
00381
00382 csfft_cox( 1 , nup , row ) ;
00383
00384 sf = 0.5 / nup ;
00385 for( ii=0 ; ii < n ; ii++ ){
00386 f[ii] = sf * row[ii].r ; g[ii] = sf * row[ii].i ;
00387 }
00388
00389 return ;
00390 }
00391
00392
00393
00394 void ft_xshear( float a , float b , int nx , int ny , float * f )
00395 {
00396 int jj , nxup ;
00397 float * fj0 , * fj1 , * zz=NULL ;
00398 float a0 , a1 ;
00399
00400 if( a == 0.0 && b == 0.0 ) return ;
00401 if( nx < 2 || ny < 1 || f == NULL ) return ;
00402
00403 if( ny%2 == 1 ){
00404 zz = (float *) malloc( sizeof(float)*nx ) ;
00405 for( jj=0 ; jj < nx ; jj++ ) zz[jj] = 0.0 ;
00406 }
00407
00408 nxup = nx ;
00409 jj = 2 ; while( jj < nxup ){ jj *= 2 ; }
00410 nxup = jj ;
00411
00412 for( jj=0 ; jj < ny ; jj+=2 ){
00413 fj0 = f + (jj*nx) ;
00414 fj1 = (jj < ny-1) ? (fj0 + nx) : zz ;
00415 a0 = a*(jj-0.5*ny) + b ;
00416 a1 = a0 + a ;
00417 ft_shift2( nx , nxup , a0 , fj0 , a1 , fj1 ) ;
00418 }
00419
00420 if( zz != NULL ) free(zz) ;
00421 return ;
00422 }
00423
00424
00425
00426 void ft_yshear( float a , float b , int nx , int ny , float * f )
00427 {
00428 int jj , nyup , ii ;
00429 float * fj0 , * fj1 ;
00430 float a0 , a1 ;
00431
00432 if( a == 0.0 && b == 0.0 ) return ;
00433 if( ny < 2 || nx < 1 || f == NULL ) return ;
00434
00435
00436
00437 fj0 = (float *) malloc( sizeof(float) * 2*ny ) ; fj1 = fj0 + ny ;
00438
00439 nyup = ny ;
00440 jj = 2 ; while( jj < nyup ){ jj *= 2 ; }
00441 nyup = jj ;
00442
00443 for( jj=0 ; jj < nx ; jj+=2 ){
00444
00445 if( jj < nx-1 ){
00446 for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = f[jj+1+ii*nx]; }
00447 } else {
00448 for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = 0.0; }
00449 }
00450
00451 a0 = a*(jj-0.5*nx) + b ;
00452 a1 = a0 + a ;
00453 ft_shift2( ny , nyup , a0 , fj0 , a1 , fj1 ) ;
00454
00455 if( jj < nx-1 ){
00456 for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; f[jj+1+ii*nx] = fj1[ii]; }
00457 } else {
00458 for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; }
00459 }
00460 }
00461
00462 free(fj0) ; return ;
00463 }
00464
00465
00466
00467 MRI_IMAGE * mri_rota_shear( MRI_IMAGE *im, float aa, float bb, float phi )
00468 {
00469 double cph , sph ;
00470 float a , b , bot,top ;
00471 MRI_IMAGE *flim ;
00472 float *flar ;
00473 int ii , nxy ;
00474
00475 if( im == NULL || ! MRI_IS_2D(im) ){
00476 fprintf(stderr,"*** mri_rota_shear only works on 2D images!\n") ; EXIT(1) ;
00477 }
00478
00479
00480
00481 if( im->kind == MRI_complex ){
00482 MRI_IMARR *impair ;
00483 MRI_IMAGE * rim , * iim , * tim ;
00484 impair = mri_complex_to_pair( im ) ;
00485 if( impair == NULL ){
00486 fprintf(stderr,"*** mri_complex_to_pair fails in mri_rota!\n") ; EXIT(1) ;
00487 }
00488 rim = IMAGE_IN_IMARR(impair,0) ;
00489 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
00490 tim = mri_rota_shear( rim , aa,bb,phi ) ; mri_free( rim ) ; rim = tim ;
00491 tim = mri_rota_shear( iim , aa,bb,phi ) ; mri_free( iim ) ; iim = tim ;
00492 flim = mri_pair_to_complex( rim , iim ) ;
00493 mri_free( rim ) ; mri_free( iim ) ;
00494 MRI_COPY_AUX(flim,im) ;
00495 return flim ;
00496 }
00497
00498
00499
00500 flim = mri_to_float( im ) ;
00501 flar = MRI_FLOAT_PTR( flim ) ;
00502
00503
00504
00505 bot = top = flar[0] ; nxy = im->nx * im->ny ;
00506 for( ii=1 ; ii < nxy ; ii++ )
00507 if( flar[ii] < bot ) bot = flar[ii] ;
00508 else if( flar[ii] > top ) top = flar[ii] ;
00509
00510
00511
00512 cph = cos(phi) ; sph = sin(phi) ;
00513
00514
00515
00516
00517 if( cph < 0.0 ){
00518 int ii , jj , top , nx=flim->nx , ny=flim->ny ;
00519 float val ;
00520
00521 top = (nx+1)/2 ;
00522 for( jj=0 ; jj < ny ; jj++ ){
00523 for( ii=1 ; ii < top ; ii++ ){
00524 val = flar[jj*nx+ii] ;
00525 flar[jj*nx+ii] = flar[jj*nx+nx-ii] ;
00526 flar[jj*nx+nx-ii] = val ;
00527 }
00528 }
00529
00530 top = (ny+1)/2 ;
00531 for( ii=0 ; ii < nx ; ii++ ){
00532 for( jj=1 ; jj < top ; jj++ ){
00533 val = flar[ii+jj*nx] ;
00534 flar[ii+jj*nx] = flar[ii+(ny-jj)*nx] ;
00535 flar[ii+(ny-jj)*nx] = val ;
00536 }
00537 }
00538
00539 cph = -cph ; sph = -sph ;
00540 }
00541
00542
00543
00544 b = sph ;
00545 a = (b != 0.0 ) ? ((cph - 1.0) / b) : (0.0) ;
00546
00547
00548
00549 ft_xshear( a , 0.0 , im->nx , im->ny , flar ) ;
00550 ft_yshear( b , bb , im->nx , im->ny , flar ) ;
00551 ft_xshear( a , aa - a*bb , im->nx , im->ny , flar ) ;
00552
00553
00554
00555 for( ii=0 ; ii < nxy ; ii++ )
00556 if( flar[ii] < bot ) flar[ii] = bot ;
00557 else if( flar[ii] > top ) flar[ii] = top ;
00558
00559 return flim ;
00560 }
00561
00562 MRI_IMAGE * mri_rota_variable( int mode, MRI_IMAGE *im, float aa, float bb, float phi )
00563 {
00564 switch(mode){
00565
00566 default:
00567 case MRI_BICUBIC: return mri_rota( im,aa,bb,phi ) ;
00568
00569 case MRI_BILINEAR: return mri_rota_bilinear( im,aa,bb,phi ) ;
00570
00571 case MRI_FOURIER: return mri_rota_shear( im,aa,bb,phi ) ;
00572 }
00573 }