1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include "mrilib.h"
8
9 /*** NOT 7D SAFE ***/
10
11 #define FINS(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
12 ? 0.0 : far[(i)+(j)*nx] )
13
14 /* cubic interpolation polynomials */
15
16 #define P_M1(x) ((x)*(1.0-(x))*((x)-2.0))
17 #define P_00(x) (3.0*((x)+1.0)*((x)-1.0)*((x)-2.0))
18 #define P_P1(x) (3.0*(x)*((x)+1.0)*(2.0-(x)))
19 #define P_P2(x) ((x)*((x)+1.0)*((x)-1.0))
20
21 /**-------------------------------------------------------------------
22 Carry out a general affine transform on an image, using bicubic
23 interpolation:
24 xnew = a11 * xold + a12 * yold + xshift
25 ynew = a21 * xold + a22 * yold + yshift
26 If the input image is MRI_complex, the output will be also;
27 otherwise, the output will be MRI_float.
28 ----------------------------------------------------------------------**/
29
mri_affine_bicubic(MRI_IMAGE * im,float a11,float a12,float a21,float a22,float xshift,float yshift)30 MRI_IMAGE *mri_affine_bicubic( MRI_IMAGE *im,
31 float a11, float a12, float a21, float a22,
32 float xshift , float yshift )
33 {
34 float rot_dx , rot_dy , rot_cph , rot_sph , top,bot,val ;
35 MRI_IMAGE *imfl , *newImg ;
36 MRI_IMARR *impair ;
37 float *far , *nar ;
38 float xx,yy , fx,fy ;
39 int ii,jj, nx,ny , ix,jy , ifx,jfy ;
40 float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
41
42 ENTRY("mri_affine_bicubic") ;
43
44 if( im == NULL || ! MRI_IS_2D(im) ){
45 fprintf(stderr,"*** mri_affine only works on 2D images!\n"); RETURN(NULL):
46 }
47
48 /* if complex image: break into pairs, do separately, reassemble */
49
50 if( im->kind == MRI_complex ){
51 MRI_IMARR *impair ;
52 MRI_IMAGE * rim , * iim , * tim ;
53 impair = mri_complex_to_pair( im ) ;
54 if( impair == NULL ){
55 fprintf(stderr,"*** mri_complex_to_pair fails in mri_affine!\n");EXIT(1);
56 }
57 rim = IMAGE_IN_IMARR(impair,0) ;
58 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
59
60 tim = mri_affine_bicubic( rim , a11,a12,a21,a22,xshift,yshift);
61 mri_free(rim) ; rim = tim ;
62
63 tim = mri_affine_bicubic( iim , a11,a12,a21,a22,xshift,yshift);
64 mri_free(iim) ; iim = tim ;
65
66 newImg = mri_pair_to_complex( rim , iim ) ;
67 mri_free( rim ) ; mri_free( iim ) ;
68 MRI_COPY_AUX(newImg,im) ;
69 RETURN( newImg );
70 }
71
72 /** movement params **/
73
74 det = a11*a22 - a21*a12 ;
75 if( fabs(det) < 1.e-5*(fabs(a11)+fabs(a12)+fabs(a21)+fabs(a22)) ){
76 fprintf(stderr,"*** input determinant=0 in mri_affine!\n");EXIT(1);
77 }
78 det = 1.0 / det ;
79
80 rot_cph = cos(phi) ; rot_sph = sin(phi) ;
81
82 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
83 -(0.5 * im->ny) * rot_sph ;
84
85 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
86 +(0.5 * im->ny) * (1.0-rot_cph) ;
87
88 /** other initialization **/
89
90 nx = im->nx ; /* image dimensions */
91 ny = im->ny ;
92
93 if( im->kind == MRI_float ) imfl = im ;
94 else imfl = mri_to_float( im ) ;
95
96 far = MRI_FLOAT_PTR(imfl) ; /* access to float data */
97 newImg = mri_new( nx , nx , MRI_float ) ; /* output image */
98 nar = MRI_FLOAT_PTR(newImg) ; /* output image data */
99
100 bot = top = far[0] ;
101 for( ii=0 ; ii < nx*ny ; ii++ )
102 if( far[ii] < bot ) bot = far[ii] ;
103 else if( far[ii] > top ) top = far[ii] ;
104
105 /*** loop over output points and warp to them ***/
106
107 for( jj=0 ; jj < nx ; jj++ ){
108 xx = rot_sph * jj + rot_dx - rot_cph ;
109 yy = rot_cph * jj + rot_dy + rot_sph ;
110 for( ii=0 ; ii < nx ; ii++ ){
111
112 xx += rot_cph ; /* get x,y in original image */
113 yy -= rot_sph ;
114
115 ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ; /* floor */
116 jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
117
118 fx = xx-ix ;
119 wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ;
120 wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
121
122 if( ix > 0 && ix < nx-2 && jy > 0 && jy < ny-2 ){
123 float * fym1, *fy00 , *fyp1 , *fyp2 ;
124
125 fym1 = far + (ix-1 + (jy-1)*nx) ;
126 fy00 = fym1 + nx ;
127 fyp1 = fy00 + nx ;
128 fyp2 = fyp1 + nx ;
129
130 f_jm1 = wt_m1 * fym1[0] + wt_00 * fym1[1]
131 + wt_p1 * fym1[2] + wt_p2 * fym1[3] ;
132
133 f_j00 = wt_m1 * fy00[0] + wt_00 * fy00[1]
134 + wt_p1 * fy00[2] + wt_p2 * fy00[3] ;
135
136 f_jp1 = wt_m1 * fyp1[0] + wt_00 * fyp1[1]
137 + wt_p1 * fyp1[2] + wt_p2 * fyp1[3] ;
138
139 f_jp2 = wt_m1 * fyp2[0] + wt_00 * fyp2[1]
140 + wt_p1 * fyp2[2] + wt_p2 * fyp2[3] ;
141
142 } else {
143
144 f_jm1 = wt_m1 * FINS(ix-1,jy-1)
145 + wt_00 * FINS(ix ,jy-1)
146 + wt_p1 * FINS(ix+1,jy-1)
147 + wt_p2 * FINS(ix+2,jy-1) ;
148
149 f_j00 = wt_m1 * FINS(ix-1,jy)
150 + wt_00 * FINS(ix ,jy)
151 + wt_p1 * FINS(ix+1,jy)
152 + wt_p2 * FINS(ix+2,jy) ;
153
154 f_jp1 = wt_m1 * FINS(ix-1,jy+1)
155 + wt_00 * FINS(ix ,jy+1)
156 + wt_p1 * FINS(ix+1,jy+1)
157 + wt_p2 * FINS(ix+2,jy+1) ;
158
159 f_jp2 = wt_m1 * FINS(ix-1,jy+2)
160 + wt_00 * FINS(ix ,jy+2)
161 + wt_p1 * FINS(ix+1,jy+2)
162 + wt_p2 * FINS(ix+2,jy+2) ;
163 }
164
165 #define THIRTYSIX 2.7777778e-2 /* 1./36.0, actually */
166
167 fy = yy-jy ;
168 val = ( P_M1(fy) * f_jm1 + P_00(fy) * f_j00
169 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) * THIRTYSIX ;
170
171 if( val < bot ) nar[ii+jj*nx] = bot ; /* too small! */
172 else if( val > top ) nar[ii+jj*nx] = top ; /* too big! */
173 else nar[ii+jj*nx] = val ; /* just right */
174
175 }
176 }
177
178 /*** cleanup and return ***/
179
180 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
181 MRI_COPY_AUX(newImg,im) ;
182 RETURN( newImg );
183 }
184
185 /**-------------------------------------------------------------------
186 Rotate and shift an image, using bilinear interpolation:
187 aa = shift in x
188 bb = shift in y
189 phi = angle in radians
190 Sort of a replacement for mri_rotate in mri_warp.c; supposed
191 to be faster. If the input image is MRI_complex, the output
192 will be also; otherwise, the output will be MRI_float.
193 ----------------------------------------------------------------------**/
194
mri_rota_bilinear(MRI_IMAGE * im,float aa,float bb,float phi)195 MRI_IMAGE *mri_rota_bilinear( MRI_IMAGE *im, float aa, float bb, float phi )
196 {
197 float rot_dx , rot_dy , rot_cph , rot_sph ;
198 MRI_IMAGE *imfl , *newImg ;
199 MRI_IMARR *impair ;
200 float *far , *nar ;
201 float xx,yy , fx,fy ;
202 int ii,jj, nx,ny , ix,jy ;
203 float f_j00,f_jp1 , wt_00,wt_p1 ;
204
205 ENTRY("mri_rota_bilinear") ;
206
207 if( im == NULL || ! MRI_IS_2D(im) ){
208 fprintf(stderr,"*** mri_rota_bilinear only works on 2D images!\n"); RETURN(NULL);
209 }
210
211 /** if complex image, break into pairs, do each separately, put back together **/
212
213 if( im->kind == MRI_complex ){
214 MRI_IMARR *impair ;
215 MRI_IMAGE * rim , * iim , * tim ;
216 impair = mri_complex_to_pair( im ) ;
217 if( impair == NULL ){
218 fprintf(stderr,"*** mri_complex_to_pair fails in mri_rota!\n") ; EXIT(1) ;
219 }
220 rim = IMAGE_IN_IMARR(impair,0) ;
221 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
222 tim = mri_rota_bilinear( rim , aa,bb,phi ) ; mri_free( rim ) ; rim = tim ;
223 tim = mri_rota_bilinear( iim , aa,bb,phi ) ; mri_free( iim ) ; iim = tim ;
224 newImg = mri_pair_to_complex( rim , iim ) ;
225 mri_free( rim ) ; mri_free( iim ) ;
226 MRI_COPY_AUX(newImg,im) ;
227 RETURN( newImg );
228 }
229
230 /** rotation params **/
231
232 rot_cph = cos(phi) ; rot_sph = sin(phi) ;
233
234 rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
235 -(0.5 * im->ny) * rot_sph ;
236
237 rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
238 +(0.5 * im->ny) * (1.0-rot_cph) ;
239
240 /** other initialization **/
241
242 nx = im->nx ; /* image dimensions */
243 ny = im->ny ;
244
245 if( im->kind == MRI_float ) imfl = im ;
246 else imfl = mri_to_float( im ) ;
247
248 far = MRI_FLOAT_PTR(imfl) ; /* access to float data */
249 newImg = mri_new( nx , nx , MRI_float ) ; /* output image */
250 nar = MRI_FLOAT_PTR(newImg) ; /* output image data */
251
252 /*** loop over output points and warp to them ***/
253
254 for( jj=0 ; jj < nx ; jj++ ){
255 xx = rot_sph * jj + rot_dx - rot_cph ;
256 yy = rot_cph * jj + rot_dy + rot_sph ;
257 for( ii=0 ; ii < nx ; ii++ ){
258
259 xx += rot_cph ; /* get x,y in original image */
260 yy -= rot_sph ;
261
262 ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ; /* floor */
263 jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
264
265 fx = xx-ix ; wt_00 = 1.0 - fx ; wt_p1 = fx ;
266
267 if( ix >= 0 && ix < nx-1 && jy >= 0 && jy < ny-1 ){
268 float *fy00 , *fyp1 ;
269
270 fy00 = far + (ix + jy*nx) ; fyp1 = fy00 + nx ;
271
272 f_j00 = wt_00 * fy00[0] + wt_p1 * fy00[1] ;
273 f_jp1 = wt_00 * fyp1[0] + wt_p1 * fyp1[1] ;
274
275 } else {
276 f_j00 = wt_00 * FINS(ix,jy ) + wt_p1 * FINS(ix+1,jy ) ;
277 f_jp1 = wt_00 * FINS(ix,jy+1) + wt_p1 * FINS(ix+1,jy+1) ;
278 }
279
280 fy = yy-jy ; nar[ii+jj*nx] = (1.0-fy) * f_j00 + fy * f_jp1 ;
281
282 }
283 }
284
285 /*** cleanup and return ***/
286
287 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
288 MRI_COPY_AUX(newImg,im) ;
289 RETURN( newImg );
290 }
291
292 /*--------------------------------------------------------------------------
293 Routines to rotate using FFTs and the shearing transformation
294 ----------------------------------------------------------------------------*/
295
296 /*** Shift 2 rows at a time with the FFT ***/
297
ft_shift2(int n,int nup,float af,float * f,float ag,float * g)298 void ft_shift2( int n, int nup, float af, float * f, float ag, float * g )
299 {
300 static int nupold=0 ;
301 static complex * row=NULL , * cf=NULL , * cg=NULL ;
302
303 int ii , nby2=nup/2 , n21=nby2+1 ;
304 complex fac , gac ;
305 float sf , sg , dk ;
306
307 /* make new memory for row storage? */
308
309 if( nup > nupold ){
310 if( row != NULL ){ free(row) ; free(cf) ; free(cg) ; }
311 row = (complex *) malloc( sizeof(complex) * nup ) ;
312 cf = (complex *) malloc( sizeof(complex) * n21 ) ;
313 cg = (complex *) malloc( sizeof(complex) * n21 ) ;
314 nupold = nup ;
315 }
316
317 /* FFT the pair of rows */
318
319 for( ii=0 ; ii < n ; ii++ ){ row[ii].r = f[ii] ; row[ii].i = g[ii] ; }
320 for( ; ii < nup ; ii++ ){ row[ii].r = row[ii].i = 0.0 ; }
321
322 csfft_cox( -1 , nup , row ) ;
323
324 /* untangle FFT coefficients from row into cf,cg */
325
326 cf[0].r = 2.0 * row[0].r ; cf[0].i = 0.0 ; /* twice too big */
327 cg[0].r = 2.0 * row[0].i ; cg[0].i = 0.0 ;
328 for( ii=1 ; ii < nby2 ; ii++ ){
329 cf[ii].r = row[ii].r + row[nup-ii].r ;
330 cf[ii].i = row[ii].i - row[nup-ii].i ;
331 cg[ii].r = row[ii].i + row[nup-ii].i ;
332 cg[ii].i = -row[ii].r + row[nup-ii].r ;
333 }
334 cf[nby2].r = 2.0 * row[nby2].r ; cf[nby2].i = 0.0 ;
335 cg[nby2].r = 2.0 * row[nby2].i ; cg[nby2].i = 0.0 ;
336
337 /* phase shift both rows (cf,cg) */
338
339 dk = (2.0*PI) / nup ;
340 sf = -af * dk ; sg = -ag * dk ;
341 for( ii=1 ; ii <= nby2 ; ii++ ){
342 fac = CEXPIT(ii*sf) ; cf[ii] = CMULT( fac , cf[ii] ) ;
343 gac = CEXPIT(ii*sg) ; cg[ii] = CMULT( gac , cg[ii] ) ;
344 }
345 cf[nby2].i = 0.0 ; cg[nby2].i = 0.0 ;
346
347 /* retangle the coefficients from 2 rows */
348
349 row[0].r = cf[0].r ; row[0].i = cg[0].r ;
350 for( ii=1 ; ii < nby2 ; ii++ ){
351 row[ii].r = cf[ii].r - cg[ii].i ;
352 row[ii].i = cf[ii].i + cg[ii].r ;
353 row[nup-ii].r = cf[ii].r + cg[ii].i ;
354 row[nup-ii].i = -cf[ii].i + cg[ii].r ;
355 }
356 row[nby2].r = cf[nby2].r ;
357 row[nby2].i = cg[nby2].r ;
358
359 /* inverse FFT and store back in output arrays */
360
361 csfft_cox( 1 , nup , row ) ;
362
363 sf = 0.5 / nup ; /* 0.5 to allow for twice too big above */
364 for( ii=0 ; ii < n ; ii++ ){
365 f[ii] = sf * row[ii].r ; g[ii] = sf * row[ii].i ;
366 }
367
368 return ;
369 }
370
371 /*** Shear in the x-direction ***/
372
ft_xshear(float a,float b,int nx,int ny,float * f)373 void ft_xshear( float a , float b , int nx , int ny , float * f )
374 {
375 int jj , nxup ;
376 float * fj0 , * fj1 , * zz=NULL ;
377 float a0 , a1 ;
378
379 if( a == 0.0 && b == 0.0 ) return ; /* nothing to do */
380 if( nx < 2 || ny < 1 || f == NULL ) return ; /* nothing to operate on */
381
382 if( ny%2 == 1 ){ /* we work in pairs, so */
383 zz = (float *) malloc( sizeof(float)*nx ) ; /* if not an even number */
384 for( jj=0 ; jj < nx ; jj++ ) zz[jj] = 0.0 ; /* of rows, make an extra */
385 }
386
387 nxup = nx ; /* min FFT length */
388 jj = 2 ; while( jj < nxup ){ jj *= 2 ; } /* next power of 2 larger */
389 nxup = jj ;
390
391 for( jj=0 ; jj < ny ; jj+=2 ){ /* shear rows in pairs */
392 fj0 = f + (jj*nx) ; /* row 0 */
393 fj1 = (jj < ny-1) ? (fj0 + nx) : zz ; /* row 1 */
394 a0 = a*(jj-0.5*ny) + b ; /* phase ramp for row 0 */
395 a1 = a0 + a ; /* phase ramp for row 1 */
396 ft_shift2( nx , nxup , a0 , fj0 , a1 , fj1 ) ;
397 }
398
399 if( zz != NULL ) free(zz) ; /* toss the trash */
400 return ;
401 }
402
403 /*** Shear in the y direction ***/
404
ft_yshear(float a,float b,int nx,int ny,float * f)405 void ft_yshear( float a , float b , int nx , int ny , float * f )
406 {
407 int jj , nyup , ii ;
408 float * fj0 , * fj1 ;
409 float a0 , a1 ;
410
411 if( a == 0.0 && b == 0.0 ) return ; /* nothing to do */
412 if( ny < 2 || nx < 1 || f == NULL ) return ; /* nothing to operate on */
413
414 /* make memory for a pair of columns */
415
416 fj0 = (float *) malloc( sizeof(float) * 2*ny ) ; fj1 = fj0 + ny ;
417
418 nyup = ny ; /* min FFT length */
419 jj = 2 ; while( jj < nyup ){ jj *= 2 ; } /* next power of 2 larger */
420 nyup = jj ;
421
422 for( jj=0 ; jj < nx ; jj+=2 ){ /* shear rows in pairs */
423
424 if( jj < nx-1 ){
425 for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = f[jj+1+ii*nx]; }
426 } else {
427 for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = 0.0; }
428 }
429
430 a0 = a*(jj-0.5*nx) + b ; /* phase ramp for row 0 */
431 a1 = a0 + a ; /* phase ramp for row 1 */
432 ft_shift2( ny , nyup , a0 , fj0 , a1 , fj1 ) ;
433
434 if( jj < nx-1 ){
435 for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; f[jj+1+ii*nx] = fj1[ii]; }
436 } else {
437 for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; }
438 }
439 }
440
441 free(fj0) ; return ;
442 }
443
444 /*** Image rotation using 3 shears ***/
445
mri_rota_shear(MRI_IMAGE * im,float aa,float bb,float phi)446 MRI_IMAGE * mri_rota_shear( MRI_IMAGE *im, float aa, float bb, float phi )
447 {
448 double cph , sph ;
449 float a , b , bot,top ;
450 MRI_IMAGE *flim ;
451 float *flar ;
452 int ii , nxy ;
453
454 ENTRY("mri_rota_shear") ;
455
456 if( im == NULL || ! MRI_IS_2D(im) ){
457 fprintf(stderr,"*** mri_rota_shear only works on 2D images!\n"); RETURN(NULL);
458 }
459
460 /** if complex image, break into pairs, do each separately, put back together **/
461
462 if( im->kind == MRI_complex ){
463 MRI_IMARR *impair ;
464 MRI_IMAGE * rim , * iim , * tim ;
465 impair = mri_complex_to_pair( im ) ;
466 if( impair == NULL ){
467 fprintf(stderr,"*** mri_complex_to_pair fails in mri_rota!\n") ; EXIT(1) ;
468 }
469 rim = IMAGE_IN_IMARR(impair,0) ;
470 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
471 tim = mri_rota_shear( rim , aa,bb,phi ) ; mri_free( rim ) ; rim = tim ;
472 tim = mri_rota_shear( iim , aa,bb,phi ) ; mri_free( iim ) ; iim = tim ;
473 flim = mri_pair_to_complex( rim , iim ) ;
474 mri_free( rim ) ; mri_free( iim ) ;
475 MRI_COPY_AUX(flim,im) ;
476 RETURN( flim );
477 }
478
479 /** copy input to output **/
480
481 flim = mri_to_float( im ) ;
482 flar = MRI_FLOAT_PTR( flim ) ;
483
484 /* find range of image data */
485
486 bot = top = flar[0] ; nxy = im->nx * im->ny ;
487 for( ii=1 ; ii < nxy ; ii++ )
488 if( flar[ii] < bot ) bot = flar[ii] ;
489 else if( flar[ii] > top ) top = flar[ii] ;
490
491 /** rotation params **/
492
493 cph = cos(phi) ; sph = sin(phi) ;
494
495 /* More than 90 degrees?
496 Must be reduced to less than 90 degrees by a 180 degree flip. */
497
498 if( cph < 0.0 ){
499 int ii , jj , top , nx=flim->nx , ny=flim->ny ;
500 float val ;
501
502 top = (nx+1)/2 ;
503 for( jj=0 ; jj < ny ; jj++ ){
504 for( ii=1 ; ii < top ; ii++ ){
505 val = flar[jj*nx+ii] ;
506 flar[jj*nx+ii] = flar[jj*nx+nx-ii] ;
507 flar[jj*nx+nx-ii] = val ;
508 }
509 }
510
511 top = (ny+1)/2 ;
512 for( ii=0 ; ii < nx ; ii++ ){
513 for( jj=1 ; jj < top ; jj++ ){
514 val = flar[ii+jj*nx] ;
515 flar[ii+jj*nx] = flar[ii+(ny-jj)*nx] ;
516 flar[ii+(ny-jj)*nx] = val ;
517 }
518 }
519
520 cph = -cph ; sph = -sph ;
521 }
522
523 /* compute shear factors for each direction */
524
525 b = sph ;
526 a = (b != 0.0 ) ? ((cph - 1.0) / b) : (0.0) ;
527
528 /* shear thrice */
529
530 ft_xshear( a , 0.0 , im->nx , im->ny , flar ) ;
531 ft_yshear( b , bb , im->nx , im->ny , flar ) ;
532 ft_xshear( a , aa - a*bb , im->nx , im->ny , flar ) ;
533
534 /* make sure data does not go out of original range */
535
536 for( ii=0 ; ii < nxy ; ii++ )
537 if( flar[ii] < bot ) flar[ii] = bot ;
538 else if( flar[ii] > top ) flar[ii] = top ;
539
540 RETURN( flim );
541 }
542
mri_rota_variable(int mode,MRI_IMAGE * im,float aa,float bb,float phi)543 MRI_IMAGE * mri_rota_variable( int mode, MRI_IMAGE *im, float aa, float bb, float phi )
544 {
545 switch(mode){
546
547 default:
548 case MRI_BICUBIC: return mri_rota( im,aa,bb,phi ) ;
549
550 case MRI_BILINEAR: return mri_rota_bilinear( im,aa,bb,phi ) ;
551
552 case MRI_FOURIER: return mri_rota_shear( im,aa,bb,phi ) ;
553 }
554 }
555