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 /*** functions to "warp" images -- not very efficient, but quite general ***/
00012 
00013 #ifdef HP
00014 # undef INLINE
00015 # pragma INLINE  xxMRI_scaler
00016 # pragma INLINE  xxMRI_rotfunc
00017 #endif
00018 
00019 #define FINS(i,j) (  ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
00020                      ? 0.0 : far[(i)+(j)*nx] )
00021 
00022 /**************************************************************************/
00023 
00024 static float sx_scale , sy_scale ;  /* global scaler data */
00025 
00026 INLINE void xxMRI_scaler( float xpr, float ypr, float *xx , float *yy )
00027 {
00028    *xx = sx_scale * xpr ;
00029    *yy = sy_scale * ypr ;
00030    return ;
00031 }
00032 
00033 /**************************************************************************/
00034 
00035 MRI_IMAGE *mri_warp( MRI_IMAGE *im , int nxnew , int nynew , int wtype ,
00036                      void wfunc(float,float,float *,float *)            )
00037 {
00038    switch( wtype ){
00039       case MRI_BILINEAR:
00040          return mri_warp_bilinear( im , nxnew , nynew , wfunc ) ;
00041 
00042       case MRI_BICUBIC:
00043          return mri_warp_bicubic( im , nxnew , nynew , wfunc ) ;
00044 
00045       default:
00046          fprintf( stderr , "mri_warp: illegal wtype %d\n" , wtype ) ;
00047          MRI_FATAL_ERROR ;
00048    }
00049    return NULL ;
00050 }
00051 
00052 /****************************************************************************/
00053 
00054 MRI_IMAGE *mri_resize_NN( MRI_IMAGE *im , int nxnew , int nynew )  /* 08 Jun 2004 */
00055 {
00056    int nx,ny , nnx,nny , ii,jj , pp,qq , bb ;
00057    float fx,fy ;
00058    MRI_IMAGE *nim ;
00059    char *nar , *ar ;
00060 
00061    if( im == NULL ) return NULL ;
00062 
00063    nx  = im->nx ;  ny  = im->ny ;
00064    nnx = nxnew  ;  nny = nynew  ;
00065    fx  = ((float)nx) / (float)nnx ;
00066    fy  = ((float)ny) / (float)nny ;
00067 
00068    nim = mri_new( nnx , nny , im->kind ) ;
00069    nar = mri_data_pointer( nim ) ;
00070    ar  = mri_data_pointer( im ) ;
00071    bb  = im->pixel_size ;
00072 
00073    for( jj=0 ; jj < nny ; jj++ ){
00074      qq = (int)( fy*jj ) ;
00075      for( ii=0 ; ii < nnx ; ii++ ){
00076        pp = (int)( fx*ii ) ;
00077        memcpy( nar + (ii+jj*nnx)*bb , ar + (pp+qq*nx)*bb , bb ) ;
00078      }
00079    }
00080 
00081    MRI_COPY_AUX(nim,im) ;
00082    nim->dx *= fx ;
00083    nim->dy *= fy ;
00084    return nim ;
00085 }
00086 
00087 /****************************************************************************/
00088 
00089 MRI_IMAGE *mri_squareaspect( MRI_IMAGE *im )  /* 08 Jun 2004 */
00090 {
00091    int nx,ny , nxnew,nynew ;
00092    float dx,dy , rr ;
00093 
00094    if( im == NULL ) return NULL ;
00095 
00096    dx = fabs(im->dx) ; dy = fabs(im->dy) ;
00097    if( dx == 0.0 || dy == 0.0 ) return NULL ;
00098    rr = dy / dx ; if( rr == 1.0 ) return NULL ;
00099 
00100    nx = im->nx ; ny = im->ny ;
00101 
00102    if( rr < 1.0 ){
00103      nxnew = rint( nx/rr ) ; if( nxnew <= nx ) return NULL ;
00104      nynew = ny ;
00105    } else {
00106      nynew = rint( ny*rr ) ; if( nynew <= ny ) return NULL ;
00107      nxnew = nx ;
00108    }
00109 
00110    return mri_resize_NN( im , nxnew , nynew ) ;
00111 }
00112 
00113 /****************************************************************************/
00114 
00115 MRI_IMAGE *mri_resize( MRI_IMAGE *im , int nxnew , int nynew )
00116 {
00117    int nx,ny , nnx,nny , wtype ;
00118 
00119    nx  = im->nx ;  ny  = im->ny ;
00120    nnx = nxnew  ;  nny = nynew  ;
00121 
00122    if( nnx <= 0 && nny <= 0 ){
00123       fprintf( stderr , "mri_resize: nxnew,nynew = %d %d\n",nxnew,nynew ) ;
00124       MRI_FATAL_ERROR ;
00125    }
00126 
00127    sx_scale = (nnx > 0) ? ((float)nx)/nnx : 0.0 ;
00128    sy_scale = (nny > 0) ? ((float)ny)/nny : 0.0 ;
00129 
00130    if( nnx <= 0 ){
00131       sx_scale = sy_scale ;
00132       nnx      = sx_scale * nx ;
00133    } else if( nny <= 0 ){
00134       sy_scale = sx_scale ;
00135       nny      = sy_scale * ny ;
00136    }
00137 
00138 
00139 #if 0
00140    wtype = MRI_BICUBIC ;
00141 
00142    if( ( ((nnx>=nx) && (nnx%nx==0)) || ((nnx<nx) && (nx%nnx==0)) ) &&
00143        ( ((nny>=ny) && (nny%ny==0)) || ((nny<ny) && (ny%nny==0)) )   ){
00144 
00145       wtype = MRI_BILINEAR ;
00146    } else {
00147       wtype = MRI_BICUBIC ;
00148    }
00149 
00150    return mri_warp( im , nnx,nny , wtype , xxMRI_scaler ) ;
00151 #else
00152    return mri_warp_bicubic( im , nnx,nny , xxMRI_scaler ) ;
00153 #endif
00154 }
00155 
00156 /**************************************************************************/
00157 
00158 MRI_IMAGE *mri_warp_bicubic( MRI_IMAGE *im , int nxnew , int nynew ,
00159                                      void wf( float,float,float *,float *) )
00160 {
00161    MRI_IMAGE *imfl , *newImg ;
00162    float *far , *nar ;
00163    float xpr,ypr , xx,yy , fx,fy ;
00164    int ii,jj, nx,ny , ix,jy ;
00165    float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
00166    float bot,top,val ;  /* 29 Mar 2003 */
00167 
00168    nx = im->nx ;  /* input image dimensions, for convenience */
00169    ny = im->ny ;
00170 
00171    nxnew = (nxnew > 0) ? nxnew : nx ;  /* default output image sizes */
00172    nynew = (nynew > 0) ? nynew : ny ;
00173 
00174    switch( im->kind ){   /* 29 Mar 2003: allow for different input types */
00175                          /*              by doing components 1 at a time */
00176      case MRI_float:
00177        imfl = im ; break ;
00178 
00179      default:
00180        imfl = mri_to_float(im) ; break ;
00181 
00182      case MRI_short:{
00183        imfl = mri_to_float(im) ;
00184        newImg  = mri_warp_bicubic( imfl , nxnew,nynew , wf ) ;
00185        mri_free(imfl) ;
00186        imfl = mri_to_mri(MRI_short,newImg) ;
00187        mri_free(newImg) ; return imfl ;
00188      }
00189 
00190      case MRI_byte:{
00191        imfl = mri_to_float(im) ;
00192        newImg  = mri_warp_bicubic( imfl , nxnew,nynew , wf ) ;
00193        mri_free(imfl) ;
00194        imfl = mri_to_mri(MRI_byte,newImg) ;
00195        mri_free(newImg) ; return imfl ;
00196      }
00197 
00198      case MRI_rgb:{
00199        MRI_IMARR *imar = mri_rgb_to_3float(im) ;
00200        MRI_IMAGE *rim,*gim,*bim ;
00201        rim = mri_warp_bicubic( IMARR_SUBIM(imar,0), nxnew,nynew, wf ) ;
00202        gim = mri_warp_bicubic( IMARR_SUBIM(imar,1), nxnew,nynew, wf ) ;
00203        bim = mri_warp_bicubic( IMARR_SUBIM(imar,2), nxnew,nynew, wf ) ;
00204        DESTROY_IMARR(imar) ;
00205        newImg = mri_3to_rgb( rim,gim,bim ) ;
00206        mri_free(rim); mri_free(gim); mri_free(bim); return newImg;
00207      }
00208 
00209    }
00210 
00211    /* at this point, imfl is in MRI_float format */
00212 
00213    far = mri_data_pointer( imfl ) ;  /* easy access to float data */
00214 
00215    newImg = mri_new( nxnew , nynew , MRI_float ) ;   /* output image */
00216    nar = mri_data_pointer( newImg ) ;                /* output image data */
00217 
00218    bot = top = far[0] ;                        /* 29 Mar 2003: */
00219    for( ii=1 ; ii < imfl->nvox ; ii++ ){       /* clip output data range */
00220           if( far[ii] > top ) top = far[ii] ;
00221      else if( far[ii] < bot ) bot = far[ii] ;
00222    }
00223 
00224    /*** loop over output points and warp to them ***/
00225 
00226    for( jj=0 ; jj < nynew ; jj++ ){
00227       ypr = jj ;
00228       for( ii=0 ; ii < nxnew ; ii++ ){
00229          xpr = ii ;
00230          wf( xpr,ypr , &xx,&yy ) ;  /* get xx,yy in original image */
00231 
00232          ix = floor( xx ) ;  fx = xx - ix ;   /* integer and       */
00233          jy = floor( yy ) ;  fy = yy - jy ;   /* fractional coords */
00234 
00235 /* define cubic interpolation polynomials and data from original grid */
00236 
00237 #define P_M1(x)  (-(x)*((x)-1)*((x)-2))
00238 #define P_00(x)  (3*((x)+1)*((x)-1)*((x)-2))
00239 #define P_P1(x)  (-3*(x)*((x)+1)*((x)-2))
00240 #define P_P2(x)  ((x)*((x)+1)*((x)-1))
00241 
00242          wt_m1 = P_M1(fx) ;  wt_00 = P_00(fx) ;  /* weights for interpolating */
00243          wt_p1 = P_P1(fx) ;  wt_p2 = P_P2(fx) ;  /* in the x-direction        */
00244 
00245          f_jm1 =  wt_m1 * FINS(ix-1,jy-1)        /* interpolate to ix + fx */
00246                 + wt_00 * FINS(ix  ,jy-1)        /* at levels jy-1 .. jy+2 */
00247                 + wt_p1 * FINS(ix+1,jy-1)
00248                 + wt_p2 * FINS(ix+2,jy-1) ;
00249 
00250          f_j00 =   wt_m1 * FINS(ix-1,jy)
00251                  + wt_00 * FINS(ix  ,jy)
00252                  + wt_p1 * FINS(ix+1,jy)
00253                  + wt_p2 * FINS(ix+2,jy) ;
00254 
00255          f_jp1 =   wt_m1 * FINS(ix-1,jy+1)
00256                  + wt_00 * FINS(ix  ,jy+1)
00257                  + wt_p1 * FINS(ix+1,jy+1)
00258                  + wt_p2 * FINS(ix+2,jy+1) ;
00259 
00260          f_jp2 =   wt_m1 * FINS(ix-1,jy+2)
00261                  + wt_00 * FINS(ix  ,jy+2)
00262                  + wt_p1 * FINS(ix+1,jy+2)
00263                  + wt_p2 * FINS(ix+2,jy+2) ;
00264 
00265          /* interpolate between y-levels to jy+fy */
00266 
00267          val = (  P_M1(fy) * f_jm1 + P_00(fy) * f_j00
00268                 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) / 36.0 ;
00269 
00270               if( val > top ) val = top ;  /* 29 Mar 2003 */
00271          else if( val < bot ) val = bot ;
00272 
00273          nar[ii+jj*nxnew] = val ;
00274       }
00275    }
00276 
00277    /*** cleanup and return ***/
00278 
00279    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00280    return newImg ;
00281 }
00282 
00283 /*************************************************************************/
00284 
00285 MRI_IMAGE *mri_warp_bilinear( MRI_IMAGE *im , int nxnew , int nynew ,
00286                                   void wf( float,float,float *,float *) )
00287 {
00288    MRI_IMAGE *imfl , *newImg ;
00289    float *far , *nar ;
00290    float xpr,ypr , xx,yy , fx,fx1,fy,fy1 , f00,f10,f01,f11 ;
00291    int ii,jj, nx,ny , ix,jy ;
00292 
00293    nx = im->nx ;  /* dimensions, for convenience */
00294    ny = im->ny ;
00295 
00296    nxnew = (nxnew > 0) ? nxnew : nx ;  /* default output image sizes */
00297    nynew = (nynew > 0) ? nynew : ny ;
00298 
00299    if( im->kind == MRI_float ){    /* convert input to float, if needed */
00300       imfl = im ;
00301    } else {
00302       imfl = mri_to_float( im ) ;
00303    }
00304    far = mri_data_pointer( imfl ) ;  /* easy access to float data */
00305 
00306    newImg = mri_new( nxnew , nynew , MRI_float ) ;   /* output image */
00307    nar = mri_data_pointer( newImg ) ;                /* output image data */
00308 
00309    /*** loop over output points and warp to them ***/
00310 
00311    for( jj=0 ; jj < nynew ; jj++ ){
00312       ypr = jj ;
00313       for( ii=0 ; ii < nxnew ; ii++ ){
00314          xpr = ii ;
00315          wf( xpr,ypr , &xx,&yy ) ;  /* get xx,yy in original image */
00316 
00317          ix = floor( xx ) ;  fx = xx - ix ;  fx1 = 1.0 - fx ;
00318          jy = floor( yy ) ;  fy = yy - jy ;  fy1 = 1.0 - fy ;
00319 
00320          f00 = FINS(ix  ,jy  ) ;
00321          f01 = FINS(ix+1,jy  ) ;
00322          f10 = FINS(ix  ,jy+1) ;
00323          f11 = FINS(ix+1,jy+1) ;
00324 
00325          nar[ii+jj*nxnew] = fx1 * ( fy1 * f00 + fy * f01 )
00326                            +fx  * ( fy1 * f10 + fy * f11 ) ;
00327 
00328       }
00329    }
00330 
00331    /*** cleanup and return ***/
00332 
00333    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00334    return newImg ;
00335 }
00336 
00337 /**********************************************************************/
00338 
00339 static float rot_dx , rot_dy , rot_cph , rot_sph ;    /* global rotfunc data */
00340 
00341 INLINE void xxMRI_rotfunc( float xpr , float ypr , float *xx , float *yy )
00342 {
00343    *xx =  rot_cph * xpr + rot_sph * ypr + rot_dx ;
00344    *yy = -rot_sph * xpr + rot_cph * ypr + rot_dy ;
00345    return ;
00346 }
00347 
00348 /**--------------------------------------------------------------------------------
00349     Rotate and shift an image, using bicubic interpolation:
00350        aa  = shift in x
00351        bb  = shift in y
00352        phi = angle in radians
00353        scl = size scale factor (0.0 --> leave the same size)
00354 -----------------------------------------------------------------------------------**/
00355 
00356 MRI_IMAGE *mri_rotate( MRI_IMAGE *im, float aa, float bb, float phi, float scl )
00357 {
00358    MRI_IMAGE  *imwarp ;
00359    int nxnew , nynew ;
00360 
00361    rot_cph = cos(phi) ;
00362    rot_sph = sin(phi) ;
00363 
00364    rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00365            -(0.5 * im->ny) * rot_sph ;
00366 
00367    rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00368            +(0.5 * im->ny) * (1.0-rot_cph) ;
00369 
00370    if( scl <= 0.0 ){
00371       nxnew = nynew = 0 ;
00372    } else {
00373       nxnew = scl * im->nx + 0.49 ;
00374       nynew = scl * im->ny + 0.49 ;
00375       rot_cph /= scl ;
00376       rot_sph /= scl ;
00377    }
00378 
00379    return mri_warp_bicubic( im , nxnew,nynew , xxMRI_rotfunc ) ;
00380 }
00381 
00382 MRI_IMAGE *mri_rotate_bilinear( MRI_IMAGE *im, float aa, float bb, float phi, float scl )
00383 {
00384    MRI_IMAGE  *imwarp ;
00385    int nxnew , nynew ;
00386 
00387    rot_cph = cos(phi) ;
00388    rot_sph = sin(phi) ;
00389 
00390    rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00391            -(0.5 * im->ny) * rot_sph ;
00392 
00393    rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00394            +(0.5 * im->ny) * (1.0-rot_cph) ;
00395 
00396    if( scl <= 0.0 ){
00397       nxnew = nynew = 0 ;
00398    } else {
00399       nxnew = scl * im->nx + 0.49 ;
00400       nynew = scl * im->ny + 0.49 ;
00401       rot_cph /= scl ;
00402       rot_sph /= scl ;
00403    }
00404 
00405    return mri_warp_bilinear( im , nxnew,nynew , xxMRI_rotfunc ) ;
00406 }
00407 
00408 #undef WARP_POINT_ROUTINES
00409 #ifdef WARP_POINT_ROUTINES
00410 
00411 /*--------------------------------------------------------------------
00412   Return one point from a warped image
00413 ----------------------------------------------------------------------*/
00414 
00415 #define INSIDE(i,j) ( ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) ? 0.0            : \
00416                     (kk==MRI_byte)    ? (float) MRI_BYTE_PIXEL(im,i,j)          : \
00417                     (kk==MRI_short)   ? (float) MRI_SHORT_PIXEL(im,i,j)         : \
00418                     (kk==MRI_int)     ? (float) MRI_INT_PIXEL(im,i,j)           : \
00419                     (kk==MRI_float)   ? (float) MRI_FLOAT_PIXEL(im,i,j)         : \
00420                     (kk==MRI_double)  ? (float) MRI_DOUBLE_PIXEL(im,i,j)        : \
00421                     (kk==MRI_complex) ? (float) CABS(MRI_COMPLEX_PIXEL(im,i,j)) : 0.0 )
00422 
00423 float mri_warp_bicubic_point( MRI_IMAGE *im , int ii , int jj ,
00424                                 void wf( float,float,float *,float *) )
00425 {
00426    float xx,yy , fx,fy , newPt ;
00427    int nx,ny , ix,jy , kk ;
00428    float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
00429 
00430    nx = im->nx ;  /* input image dimensions, for convenience */
00431    ny = im->ny ;
00432    kk = im->kind ;
00433 
00434    /*** warp to output point ***/
00435 
00436    wf( (float) ii , (float) jj , &xx,&yy ) ;  /* get xx,yy in original image */
00437 
00438    ix = floor( xx ) ;  fx = xx - ix ;   /* integer and       */
00439    jy = floor( yy ) ;  fy = yy - jy ;   /* fractional coords */
00440 
00441    wt_m1 = P_M1(fx) ;  wt_00 = P_00(fx) ;  /* weights for interpolating */
00442    wt_p1 = P_P1(fx) ;  wt_p2 = P_P2(fx) ;  /* in the x-direction        */
00443 
00444    f_jm1 =  wt_m1 * INSIDE(ix-1,jy-1)        /* interpolate to ix + fx */
00445           + wt_00 * INSIDE(ix  ,jy-1)        /* at levels jy-1 .. jy+2 */
00446           + wt_p1 * INSIDE(ix+1,jy-1)
00447           + wt_p2 * INSIDE(ix+2,jy-1) ;
00448 
00449    f_j00 =   wt_m1 * INSIDE(ix-1,jy)
00450            + wt_00 * INSIDE(ix  ,jy)
00451            + wt_p1 * INSIDE(ix+1,jy)
00452            + wt_p2 * INSIDE(ix+2,jy) ;
00453 
00454    f_jp1 =   wt_m1 * INSIDE(ix-1,jy+1)
00455            + wt_00 * INSIDE(ix  ,jy+1)
00456            + wt_p1 * INSIDE(ix+1,jy+1)
00457            + wt_p2 * INSIDE(ix+2,jy+1) ;
00458 
00459    f_jp2 =   wt_m1 * INSIDE(ix-1,jy+2)
00460            + wt_00 * INSIDE(ix  ,jy+2)
00461            + wt_p1 * INSIDE(ix+1,jy+2)
00462            + wt_p2 * INSIDE(ix+2,jy+2) ;
00463 
00464    /* interpolate between y-levels to jy+fy */
00465 
00466    newPt = (  P_M1(fy) * f_jm1 + P_00(fy) * f_j00
00467           + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) / 36.0 ;
00468 
00469    return newPt ;
00470 }
00471 
00472 /*--------------------------------------------------------------------
00473   Return one point from a rotated image
00474 ----------------------------------------------------------------------*/
00475 
00476 float mri_rotate_point( MRI_IMAGE *im, float aa, float bb, float phi, float scl ,
00477                         int ix , int jy )
00478 {
00479    MRI_IMAGE  *imwarp ;
00480    int nxnew , nynew ;
00481 
00482    rot_cph = cos(phi) ;
00483    rot_sph = sin(phi) ;
00484 
00485    rot_dx = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00486            -(0.5 * im->ny) * rot_sph ;
00487 
00488    rot_dy = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00489            +(0.5 * im->ny) * (1.0-rot_cph) ;
00490 
00491    if( scl > 0.0 ){
00492       rot_cph /= scl ;
00493       rot_sph /= scl ;
00494    }
00495 
00496    return mri_warp_bicubic_point( im , ix,jy , xxMRI_rotfunc ) ;
00497 }
00498 #endif /* WARP_POINT_ROUTINES */