00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include "mrilib.h"
00008 
00009 /*** NOT 7D SAFE ***/
00010 
00011 /*! Return far[] if (i,j) is inside the image, otherwise return 0. */
00012 
00013 #define FINS(i,j) (  ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00014                      ? 0.0 : far[(i)+(j)*nx] )
00015 
00016 /*-------------------------------------------------------------------*/
00017 /*! Invert a 2D matrix, dude. */
00018 
00019 static void invert2d( float  axx, float  axy, float  ayx, float  ayy ,
00020                       float *bxx, float *bxy, float *byx, float *byy  )
00021 {
00022    float det = axx*ayy - axy*ayx ;
00023    if( det == 0.0 ){ *bxx=*byy=*bxy=*byx = 0.0 ; return ; }
00024    *bxx =  ayy / det ;
00025    *byy =  axx / det ;
00026    *bxy = -axy / det ;
00027    *byx = -ayx / det ;
00028    return ;
00029 }
00030 
00031 /*-------------------------------------------------------------------*/
00032 /*! Affine transform a 2D image, using bilinear interpolation:
00033     If flag == 0
00034        [ xout ] = [ axx axy ] [ xin ]
00035        [ yout ] = [ ayx ayy ] [ yin ]
00036     If flag == 1
00037        [ xin ] = [ axx axy ] [ xout ]
00038        [ yin ] = [ ayx ayy ] [ yout ]
00039     These are index coordinates, not spatial.
00040 ----------------------------------------------------------------------*/
00041 
00042 MRI_IMAGE *mri_aff2d_byte( MRI_IMAGE *im, int flag ,
00043                            float axx, float axy, float ayx, float ayy )
00044 {
00045    float bxx,bxy,byx,byy , xbase,ybase , xx,yy , fx,fy ;
00046    float f_j00,f_jp1 , wt_00,wt_p1 ;
00047    int ii,jj , nx,ny , ix,jy ;
00048    MRI_IMAGE *newImg ;
00049    byte *far , *nar ;
00050 
00051 ENTRY("mri_aff2d_byte") ;
00052 
00053    if( im == NULL || !MRI_IS_2D(im) || im->kind != MRI_byte ){
00054       fprintf(stderr,"*** mri_aff2d_byte only works on 2D byte images!\n");
00055       RETURN( NULL );
00056    }
00057 
00058    if( flag == 0 ){
00059       invert2d( axx,axy,ayx,ayy , &bxx,&bxy,&byx,&byy ) ;
00060    } else {
00061       bxx = axx ; bxy = axy ; byx = ayx ; byy = ayy ;
00062    }
00063    if( (bxx == 0.0 && bxy == 0.0) || (byx == 0.0 && byy == 0.0) ){
00064       fprintf(stderr,"*** mri_aff2d_byte: input matrix is singular!\n") ;
00065       RETURN( NULL );
00066    }
00067 
00068    nx = im->nx ; ny = im->ny ;
00069    xbase = 0.5*nx*(1.0-bxx) - 0.5*ny*bxy ;
00070    ybase = 0.5*ny*(1.0-byy) - 0.5*nx*byx ;
00071 
00072    far = MRI_BYTE_PTR(im) ;                /* input image data */
00073    newImg = mri_new( nx , nx , MRI_byte ) ;   /* output image */
00074    nar = MRI_BYTE_PTR(newImg) ;               /* output image data */
00075 
00076    /*** loop over output points and warp to them ***/
00077 
00078    for( jj=0 ; jj < nx ; jj++ ){
00079       xx = xbase-bxx + bxy * jj ;
00080       yy = ybase-byx + byy * jj ;
00081       for( ii=0 ; ii < nx ; ii++ ){
00082 
00083          xx += bxx ;  /* get x,y in original image */
00084          yy += byx ;
00085 
00086          ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ;  /* floor */
00087          jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
00088 
00089          fx = xx-ix ; wt_00 = 1.0 - fx ; wt_p1 = fx ;
00090 
00091          if( ix >= 0 && ix < nx-1 && jy >= 0 && jy < ny-1 ){
00092             byte *fy00 , *fyp1 ;
00093 
00094             fy00 = far + (ix + jy*nx) ; fyp1 = fy00 + nx ;
00095 
00096             f_j00 = wt_00 * fy00[0] + wt_p1 * fy00[1] ;
00097             f_jp1 = wt_00 * fyp1[0] + wt_p1 * fyp1[1] ;
00098 
00099          } else {
00100             f_j00 = wt_00 * FINS(ix,jy  ) + wt_p1 * FINS(ix+1,jy  ) ;
00101             f_jp1 = wt_00 * FINS(ix,jy+1) + wt_p1 * FINS(ix+1,jy+1) ;
00102          }
00103 
00104          fy  = yy-jy ; nar[ii+jj*nx] = (1.0-fy) * f_j00 + fy * f_jp1 ;
00105 
00106       }
00107    }
00108 
00109    MRI_COPY_AUX(newImg,im) ;
00110    RETURN( newImg ) ;
00111 }
00112 
00113 /*----------------------------------------------------------------------*/
00114 
00115 /*! Same as FINS(), but for the R component of an RGB image. */
00116 
00117 #define RINS(i,j) (  ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00118                      ? 0.0 : far[3*((i)+(j)*nx)] )
00119 
00120 /*! Same as FINS(), but for the G component of an RGB image. */
00121 
00122 #define GINS(i,j) (  ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00123                      ? 0.0 : far[3*((i)+(j)*nx)+1] )
00124 
00125 /*! Same as FINS(), but for the B component of an RGB image. */
00126 
00127 #define BINS(i,j) (  ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00128                      ? 0.0 : far[3*((i)+(j)*nx)+2] )
00129 
00130 /*----------------------------------------------------------------------*/
00131 /*!  Same as mri_aff2d_byte(), but for RGB images [11 Dec 2000].
00132 ------------------------------------------------------------------------*/
00133 
00134 MRI_IMAGE *mri_aff2d_rgb( MRI_IMAGE *im, int flag ,
00135                           float axx, float axy, float ayx, float ayy )
00136 {
00137    float bxx,bxy,byx,byy , xbase,ybase , xx,yy , fx,fy ;
00138    float f_j00r,f_jp1r , f_j00g,f_jp1g , f_j00b,f_jp1b , wt_00,wt_p1 ;
00139    int jj , nx,ny , ix,jy ;
00140    MRI_IMAGE *newImg ;
00141    byte *far , *nar ;
00142    register int ii ;
00143 
00144 ENTRY("mri_aff2d_rgb") ;
00145 
00146    if( im == NULL || !MRI_IS_2D(im) || im->kind != MRI_rgb ){
00147       fprintf(stderr,"*** mri_aff2d_rgb only works on 2D RGB images!\n");
00148       RETURN( NULL );
00149    }
00150 
00151    if( flag == 0 ){
00152       invert2d( axx,axy,ayx,ayy , &bxx,&bxy,&byx,&byy ) ;
00153    } else {
00154       bxx = axx ; bxy = axy ; byx = ayx ; byy = ayy ;
00155    }
00156    if( (bxx == 0.0 && bxy == 0.0) || (byx == 0.0 && byy == 0.0) ){
00157       fprintf(stderr,"*** mri_aff2d_byte: input matrix is singular!\n") ;
00158       RETURN( NULL );
00159    }
00160 
00161    nx = im->nx ; ny = im->ny ;
00162    xbase = 0.5*nx*(1.0-bxx) - 0.5*ny*bxy ;
00163    ybase = 0.5*ny*(1.0-byy) - 0.5*nx*byx ;
00164 
00165    far = MRI_RGB_PTR(im) ;                /* input image data */
00166    newImg = mri_new( nx , nx , MRI_rgb ) ;   /* output image */
00167    nar = MRI_RGB_PTR(newImg) ;               /* output image data */
00168 
00169    /*** loop over output points and warp to them ***/
00170 
00171    for( jj=0 ; jj < nx ; jj++ ){
00172       xx = xbase-bxx + bxy * jj ;
00173       yy = ybase-byx + byy * jj ;
00174       for( ii=0 ; ii < nx ; ii++ ){
00175 
00176          xx += bxx ;  /* get x,y in original image */
00177          yy += byx ;
00178 
00179          ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ;  /* floor */
00180          jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
00181 
00182          fx = xx-ix ; wt_00 = 1.0 - fx ; wt_p1 = fx ;
00183 
00184          if( ix >= 0 && ix < nx-1 && jy >= 0 && jy < ny-1 ){
00185             byte *fy00 , *fyp1 ;
00186 
00187             fy00 = far + 3*(ix+jy*nx) ; fyp1 = fy00 + 3*nx ;
00188 
00189             f_j00r = wt_00 * fy00[0] + wt_p1 * fy00[3] ;
00190             f_j00g = wt_00 * fy00[1] + wt_p1 * fy00[4] ;
00191             f_j00b = wt_00 * fy00[2] + wt_p1 * fy00[5] ;
00192 
00193             f_jp1r = wt_00 * fyp1[0] + wt_p1 * fyp1[3] ;
00194             f_jp1g = wt_00 * fyp1[1] + wt_p1 * fyp1[4] ;
00195             f_jp1b = wt_00 * fyp1[2] + wt_p1 * fyp1[5] ;
00196 
00197          } else {
00198             f_j00r = wt_00 * RINS(ix,jy  ) + wt_p1 * RINS(ix+1,jy  ) ;
00199             f_j00g = wt_00 * GINS(ix,jy  ) + wt_p1 * GINS(ix+1,jy  ) ;
00200             f_j00b = wt_00 * BINS(ix,jy  ) + wt_p1 * BINS(ix+1,jy  ) ;
00201 
00202             f_jp1r = wt_00 * RINS(ix,jy+1) + wt_p1 * RINS(ix+1,jy+1) ;
00203             f_jp1g = wt_00 * GINS(ix,jy+1) + wt_p1 * GINS(ix+1,jy+1) ;
00204             f_jp1b = wt_00 * BINS(ix,jy+1) + wt_p1 * BINS(ix+1,jy+1) ;
00205          }
00206 
00207          fy = yy-jy ;
00208          nar[3*ii+ 3*jj*nx   ] = (1.0-fy) * f_j00r + fy * f_jp1r ;
00209          nar[3*ii+(3*jj*nx+1)] = (1.0-fy) * f_j00g + fy * f_jp1g ;
00210          nar[3*ii+(3*jj*nx+2)] = (1.0-fy) * f_j00b + fy * f_jp1b ;
00211 
00212       }
00213    }
00214 
00215    MRI_COPY_AUX(newImg,im) ;
00216    RETURN( newImg );
00217 }