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