Doxygen Source Code Documentation
mri_affine.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | FINS(i, j) |
#define | P_M1(x) ((x)*(1.0-(x))*((x)-2.0)) |
#define | P_00(x) (3.0*((x)+1.0)*((x)-1.0)*((x)-2.0)) |
#define | P_P1(x) (3.0*(x)*((x)+1.0)*(2.0-(x))) |
#define | P_P2(x) ((x)*((x)+1.0)*((x)-1.0)) |
#define | THIRTYSIX 2.7777778e-2 |
Functions | |
MRI_IMAGE * | mri_affine_bicubic (MRI_IMAGE *im, float a11, float a12, float a21, float a22, float xshift, float yshift) |
MRI_IMAGE * | mri_rota_bilinear (MRI_IMAGE *im, float aa, float bb, float phi) |
void | ft_shift2 (int n, int nup, float af, float *f, float ag, float *g) |
void | ft_xshear (float a, float b, int nx, int ny, float *f) |
void | ft_yshear (float a, float b, int nx, int ny, float *f) |
MRI_IMAGE * | mri_rota_shear (MRI_IMAGE *im, float aa, float bb, float phi) |
MRI_IMAGE * | mri_rota_variable (int mode, MRI_IMAGE *im, float aa, float bb, float phi) |
Define Documentation
|
Value: ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \ ? 0.0 : far[(i)+(j)*nx] ) Definition at line 11 of file mri_affine.c. |
|
Definition at line 17 of file mri_affine.c. Referenced by cub_shift(), mri_affine_bicubic(), mri_rota(), mri_warp3D_cubic(), mri_warp_bicubic(), and shifter(). |
|
Definition at line 16 of file mri_affine.c. Referenced by cub_shift(), mri_affine_bicubic(), mri_rota(), mri_warp3D_cubic(), mri_warp_bicubic(), and shifter(). |
|
Definition at line 18 of file mri_affine.c. Referenced by cub_shift(), mri_affine_bicubic(), mri_rota(), mri_warp3D_cubic(), mri_warp_bicubic(), and shifter(). |
|
Definition at line 19 of file mri_affine.c. Referenced by cub_shift(), mri_affine_bicubic(), mri_rota(), mri_warp3D_cubic(), mri_warp_bicubic(), and shifter(). |
|
|
Function Documentation
|
other initialization * Definition at line 298 of file mri_affine.c.
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 /* make new memory for row storage? */ 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 /* FFT the pair of rows */ 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 /* untangle FFT coefficients from row into cf,cg */ 00325 00326 cf[0].r = 2.0 * row[0].r ; cf[0].i = 0.0 ; /* twice too big */ 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 /* phase shift both rows (cf,cg) */ 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 /* retangle the coefficients from 2 rows */ 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 /* inverse FFT and store back in output arrays */ 00360 00361 csfft_cox( 1 , nup , row ) ; 00362 00363 sf = 0.5 / nup ; /* 0.5 to allow for twice too big above */ 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 } |
|
Definition at line 373 of file mri_affine.c. References a, free, ft_shift2(), and malloc.
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 ; /* nothing to do */ 00380 if( nx < 2 || ny < 1 || f == NULL ) return ; /* nothing to operate on */ 00381 00382 if( ny%2 == 1 ){ /* we work in pairs, so */ 00383 zz = (float *) malloc( sizeof(float)*nx ) ; /* if not an even number */ 00384 for( jj=0 ; jj < nx ; jj++ ) zz[jj] = 0.0 ; /* of rows, make an extra */ 00385 } 00386 00387 nxup = nx ; /* min FFT length */ 00388 jj = 2 ; while( jj < nxup ){ jj *= 2 ; } /* next power of 2 larger */ 00389 nxup = jj ; 00390 00391 for( jj=0 ; jj < ny ; jj+=2 ){ /* shear rows in pairs */ 00392 fj0 = f + (jj*nx) ; /* row 0 */ 00393 fj1 = (jj < ny-1) ? (fj0 + nx) : zz ; /* row 1 */ 00394 a0 = a*(jj-0.5*ny) + b ; /* phase ramp for row 0 */ 00395 a1 = a0 + a ; /* phase ramp for row 1 */ 00396 ft_shift2( nx , nxup , a0 , fj0 , a1 , fj1 ) ; 00397 } 00398 00399 if( zz != NULL ) free(zz) ; /* toss the trash */ 00400 return ; 00401 } |
|
Definition at line 405 of file mri_affine.c. References a, free, ft_shift2(), and malloc.
00406 { 00407 int jj , nyup , ii ; 00408 float * fj0 , * fj1 ; 00409 float a0 , a1 ; 00410 00411 if( a == 0.0 && b == 0.0 ) return ; /* nothing to do */ 00412 if( ny < 2 || nx < 1 || f == NULL ) return ; /* nothing to operate on */ 00413 00414 /* make memory for a pair of columns */ 00415 00416 fj0 = (float *) malloc( sizeof(float) * 2*ny ) ; fj1 = fj0 + ny ; 00417 00418 nyup = ny ; /* min FFT length */ 00419 jj = 2 ; while( jj < nyup ){ jj *= 2 ; } /* next power of 2 larger */ 00420 nyup = jj ; 00421 00422 for( jj=0 ; jj < nx ; jj+=2 ){ /* shear rows in pairs */ 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 ; /* phase ramp for row 0 */ 00431 a1 = a0 + a ; /* phase ramp for row 1 */ 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 } |
|
------------------------------------------------------------------- Carry out a general affine transform on an image, using bicubic interpolation: xnew = a11 * xold + a12 * yold + xshift ynew = a21 * xold + a22 * yold + yshift If the input image is MRI_complex, the output will be also; otherwise, the output will be MRI_float. ----------------------------------------------------------------------* Definition at line 30 of file mri_affine.c. References ENTRY, EXIT, far, FINS, FREE_IMARR, IMAGE_IN_IMARR, MRI_IMAGE::kind, mri_complex_to_pair(), MRI_COPY_AUX, MRI_FLOAT_PTR, mri_free(), MRI_IS_2D, mri_new(), mri_pair_to_complex(), mri_to_float(), MRI_IMAGE::nx, MRI_IMAGE::ny, P_00, P_M1, P_P1, P_P2, RETURN, and top.
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 /* if complex image: break into pairs, do separately, reassemble */ 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 /** movement params **/ 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 /** other initialization **/ 00089 00090 nx = im->nx ; /* image dimensions */ 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) ; /* access to float data */ 00097 newImg = mri_new( nx , nx , MRI_float ) ; /* output image */ 00098 nar = MRI_FLOAT_PTR(newImg) ; /* output image data */ 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 /*** loop over output points and warp to them ***/ 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 ; /* get x,y in original image */ 00113 yy -= rot_sph ; 00114 00115 ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ; /* floor */ 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 /* 1./36.0, actually */ 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 ; /* too small! */ 00172 else if( val > top ) nar[ii+jj*nx] = top ; /* too big! */ 00173 else nar[ii+jj*nx] = val ; /* just right */ 00174 00175 } 00176 } 00177 00178 /*** cleanup and return ***/ 00179 00180 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */ 00181 MRI_COPY_AUX(newImg,im) ; 00182 RETURN( newImg ); 00183 } |
|
------------------------------------------------------------------- Rotate and shift an image, using bilinear interpolation: aa = shift in x bb = shift in y phi = angle in radians Sort of a replacement for mri_rotate in mri_warp.c; supposed to be faster. If the input image is MRI_complex, the output will be also; otherwise, the output will be MRI_float. ----------------------------------------------------------------------* Definition at line 195 of file mri_affine.c.
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 /** if complex image, break into pairs, do each separately, put back together **/ 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 /** rotation params **/ 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 /** other initialization **/ 00241 00242 nx = im->nx ; /* image dimensions */ 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) ; /* access to float data */ 00249 newImg = mri_new( nx , nx , MRI_float ) ; /* output image */ 00250 nar = MRI_FLOAT_PTR(newImg) ; /* output image data */ 00251 00252 /*** loop over output points and warp to them ***/ 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 ; /* get x,y in original image */ 00260 yy -= rot_sph ; 00261 00262 ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ; /* floor */ 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 /*** cleanup and return ***/ 00286 00287 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */ 00288 MRI_COPY_AUX(newImg,im) ; 00289 RETURN( newImg ); 00290 } |
|
Definition at line 446 of file mri_affine.c.
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 /** if complex image, break into pairs, do each separately, put back together **/ 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 /** copy input to output **/ 00480 00481 flim = mri_to_float( im ) ; 00482 flar = MRI_FLOAT_PTR( flim ) ; 00483 00484 /* find range of image data */ 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 /** rotation params **/ 00492 00493 cph = cos(phi) ; sph = sin(phi) ; 00494 00495 /* More than 90 degrees? 00496 Must be reduced to less than 90 degrees by a 180 degree flip. */ 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 /* compute shear factors for each direction */ 00524 00525 b = sph ; 00526 a = (b != 0.0 ) ? ((cph - 1.0) / b) : (0.0) ; 00527 00528 /* shear thrice */ 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 /* make sure data does not go out of original range */ 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 } |
|
rotation params * Definition at line 543 of file mri_affine.c.
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 } |