00001 #include "mrilib.h"
00002 
00003 /*** NOT 7D SAFE ***/
00004 
00005 #define FINS(i,j) (  ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00006                      ? 0.0 : far[(i)+(j)*nx] )
00007 
00008 /**-------------------------------------------------------------------
00009     Shift an image, using bilinear interpolation:
00010        aa  = shift in x
00011        bb  = shift in y
00012     Like mri_rota_bilinear, with phi=0 -- RWCox - 11 Sep 2001
00013 ----------------------------------------------------------------------**/
00014 
00015 MRI_IMAGE *mri_shift2D_bilinear( MRI_IMAGE *im, float aa, float bb )
00016 {
00017    float dx , dy ;
00018    MRI_IMAGE *imfl , *newImg ;
00019    MRI_IMARR *impair ;
00020    float *far , *nar ;
00021    float xx,yy , fx,fy ;
00022    int ii,jj, nx,ny , ix,jy ;
00023    float f_j00,f_jp1 , wt_00,wt_p1 ;
00024 
00025 ENTRY("mri_shift2D_bilinear") ;
00026 
00027    if( im == NULL || !MRI_IS_2D(im) ){
00028       fprintf(stderr,"*** mri_shift2D_bilinear only works on 2D images!\n") ;
00029       EXIT(1) ;
00030    }
00031 
00032    /** if complex image, break into pairs, do each separately, put back together **/
00033 
00034    if( im->kind == MRI_complex ){
00035       MRI_IMARR *impair ;
00036       MRI_IMAGE *rim , *iim , *tim ;
00037       impair = mri_complex_to_pair( im ) ;
00038       if( impair == NULL ){
00039          fprintf(stderr,"*** mri_complex_to_pair fails in mri_shift2D_bilinear!\n");
00040          EXIT(1) ;
00041       }
00042       rim = IMAGE_IN_IMARR(impair,0) ;
00043       iim = IMAGE_IN_IMARR(impair,1) ;  FREE_IMARR(impair) ;
00044       tim = mri_shift2D_bilinear( rim , aa,bb ); mri_free(rim); rim = tim;
00045       tim = mri_shift2D_bilinear( iim , aa,bb ); mri_free(iim); iim = tim;
00046       newImg = mri_pair_to_complex( rim , iim ) ;
00047       mri_free(rim) ; mri_free(iim) ;
00048       MRI_COPY_AUX(newImg,im) ;
00049       RETURN(newImg) ;
00050    }
00051 
00052    /** shift params **/
00053 
00054    dx = - aa ;
00055    dy = - bb ;
00056 
00057    /** other initialization **/
00058 
00059    nx = im->nx ;  /* image dimensions */
00060    ny = im->ny ;
00061 
00062    if( im->kind == MRI_float ) imfl = im ;
00063    else                        imfl = mri_to_float( im ) ;
00064 
00065    far = MRI_FLOAT_PTR(imfl) ;              /* access to float data */
00066    newImg = mri_new( nx , nx , MRI_float ) ;   /* output image */
00067    nar = MRI_FLOAT_PTR(newImg) ;               /* output image data */
00068 
00069    /*** loop over output points and warp to them ***/
00070 
00071    for( jj=0 ; jj < nx ; jj++ ){
00072       xx = dx - 1.0 ;
00073       yy = jj + dy ;
00074 
00075       jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
00076 
00077       for( ii=0 ; ii < nx ; ii++ ){
00078 
00079          xx += 1.0 ;  /* get x,y in original image */
00080 
00081          ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ;  /* floor */
00082 
00083          fx = xx-ix ; wt_00 = 1.0 - fx ; wt_p1 = fx ;
00084 
00085          if( ix >= 0 && ix < nx-1 && jy >= 0 && jy < ny-1 ){
00086             float *fy00 , *fyp1 ;
00087 
00088             fy00 = far + (ix + jy*nx) ; fyp1 = fy00 + nx ;
00089 
00090             f_j00 = wt_00 * fy00[0] + wt_p1 * fy00[1] ;
00091             f_jp1 = wt_00 * fyp1[0] + wt_p1 * fyp1[1] ;
00092 
00093          } else {
00094             f_j00 = wt_00 * FINS(ix,jy  ) + wt_p1 * FINS(ix+1,jy  ) ;
00095             f_jp1 = wt_00 * FINS(ix,jy+1) + wt_p1 * FINS(ix+1,jy+1) ;
00096          }
00097 
00098          fy  = yy-jy ; nar[ii+jj*nx] = (1.0-fy) * f_j00 + fy * f_jp1 ;
00099 
00100       }
00101    }
00102 
00103    /*** cleanup and return ***/
00104 
00105    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00106    MRI_COPY_AUX(newImg,im) ;
00107    RETURN(newImg) ;
00108 }