1 #include "mrilib.h"
2
3 /*** NOT 7D SAFE ***/
4
5 #define FINS(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
6 ? 0.0 : far[(i)+(j)*nx] )
7
8 /**-------------------------------------------------------------------
9 Shift an image, using bilinear interpolation:
10 aa = shift in x
11 bb = shift in y
12 Like mri_rota_bilinear, with phi=0 -- RWCox - 11 Sep 2001
13 ----------------------------------------------------------------------**/
14
mri_shift2D_bilinear(MRI_IMAGE * im,float aa,float bb)15 MRI_IMAGE *mri_shift2D_bilinear( MRI_IMAGE *im, float aa, float bb )
16 {
17 float dx , dy ;
18 MRI_IMAGE *imfl , *newImg ;
19 MRI_IMARR *impair ;
20 float *far , *nar ;
21 float xx,yy , fx,fy ;
22 int ii,jj, nx,ny , ix,jy ;
23 float f_j00,f_jp1 , wt_00,wt_p1 ;
24
25 ENTRY("mri_shift2D_bilinear") ;
26
27 if( im == NULL || !MRI_IS_2D(im) ){
28 fprintf(stderr,"*** mri_shift2D_bilinear only works on 2D images!\n") ;
29 EXIT(1) ;
30 }
31
32 /** if complex image, break into pairs, do each separately, put back together **/
33
34 if( im->kind == MRI_complex ){
35 MRI_IMARR *impair ;
36 MRI_IMAGE *rim , *iim , *tim ;
37 impair = mri_complex_to_pair( im ) ;
38 if( impair == NULL ){
39 fprintf(stderr,"*** mri_complex_to_pair fails in mri_shift2D_bilinear!\n");
40 EXIT(1) ;
41 }
42 rim = IMAGE_IN_IMARR(impair,0) ;
43 iim = IMAGE_IN_IMARR(impair,1) ; FREE_IMARR(impair) ;
44 tim = mri_shift2D_bilinear( rim , aa,bb ); mri_free(rim); rim = tim;
45 tim = mri_shift2D_bilinear( iim , aa,bb ); mri_free(iim); iim = tim;
46 newImg = mri_pair_to_complex( rim , iim ) ;
47 mri_free(rim) ; mri_free(iim) ;
48 MRI_COPY_AUX(newImg,im) ;
49 RETURN(newImg) ;
50 }
51
52 /** shift params **/
53
54 dx = - aa ;
55 dy = - bb ;
56
57 /** other initialization **/
58
59 nx = im->nx ; /* image dimensions */
60 ny = im->ny ;
61
62 if( im->kind == MRI_float ) imfl = im ;
63 else imfl = mri_to_float( im ) ;
64
65 far = MRI_FLOAT_PTR(imfl) ; /* access to float data */
66 newImg = mri_new( nx , nx , MRI_float ) ; /* output image */
67 nar = MRI_FLOAT_PTR(newImg) ; /* output image data */
68
69 /*** loop over output points and warp to them ***/
70
71 for( jj=0 ; jj < nx ; jj++ ){
72 xx = dx - 1.0 ;
73 yy = jj + dy ;
74
75 jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
76
77 for( ii=0 ; ii < nx ; ii++ ){
78
79 xx += 1.0 ; /* get x,y in original image */
80
81 ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ; /* floor */
82
83 fx = xx-ix ; wt_00 = 1.0 - fx ; wt_p1 = fx ;
84
85 if( ix >= 0 && ix < nx-1 && jy >= 0 && jy < ny-1 ){
86 float *fy00 , *fyp1 ;
87
88 fy00 = far + (ix + jy*nx) ; fyp1 = fy00 + nx ;
89
90 f_j00 = wt_00 * fy00[0] + wt_p1 * fy00[1] ;
91 f_jp1 = wt_00 * fyp1[0] + wt_p1 * fyp1[1] ;
92
93 } else {
94 f_j00 = wt_00 * FINS(ix,jy ) + wt_p1 * FINS(ix+1,jy ) ;
95 f_jp1 = wt_00 * FINS(ix,jy+1) + wt_p1 * FINS(ix+1,jy+1) ;
96 }
97
98 fy = yy-jy ; nar[ii+jj*nx] = (1.0-fy) * f_j00 + fy * f_jp1 ;
99
100 }
101 }
102
103 /*** cleanup and return ***/
104
105 if( im != imfl ) mri_free(imfl) ; /* throw away unneeded workspace */
106 MRI_COPY_AUX(newImg,im) ;
107 RETURN(newImg) ;
108 }
109