(  ( (i)<0 || (j)<0 || (i)>=nx || (j)>=ny ) \
                     ? 0.0 : far[(i)+(j)*nx] )

00299 {
00300    static int nupold=0 ;
00301    static complex * row=NULL , * cf=NULL , * cg=NULL ;
00302 
00303    int ii , nby2=nup/2 , n21=nby2+1 ;
00304    complex fac , gac ;
00305    float sf , sg , dk ;
00306 
00307    /* make new memory for row storage? */
00308 
00309    if( nup > nupold ){
00310       if( row != NULL ){ free(row) ; free(cf) ; free(cg) ; }
00311       row = (complex *) malloc( sizeof(complex) * nup ) ;
00312       cf  = (complex *) malloc( sizeof(complex) * n21 ) ;
00313       cg  = (complex *) malloc( sizeof(complex) * n21 ) ;
00314       nupold = nup ;
00315    }
00316 
00317    /* FFT the pair of rows */
00318 
00319    for( ii=0 ; ii < n   ; ii++ ){ row[ii].r = f[ii] ; row[ii].i = g[ii] ; }
00320    for(      ; ii < nup ; ii++ ){ row[ii].r = row[ii].i = 0.0 ; }
00321 
00322    csfft_cox( -1 , nup , row ) ;
00323 
00324    /* untangle FFT coefficients from row into cf,cg */
00325 
00326    cf[0].r = 2.0 * row[0].r ; cf[0].i = 0.0 ;  /* twice too big */
00327    cg[0].r = 2.0 * row[0].i ; cg[0].i = 0.0 ;
00328    for( ii=1 ; ii < nby2 ; ii++ ){
00329       cf[ii].r =  row[ii].r + row[nup-ii].r ;
00330       cf[ii].i =  row[ii].i - row[nup-ii].i ;
00331       cg[ii].r =  row[ii].i + row[nup-ii].i ;
00332       cg[ii].i = -row[ii].r + row[nup-ii].r ;
00333    }
00334    cf[nby2].r = 2.0 * row[nby2].r ; cf[nby2].i = 0.0 ;
00335    cg[nby2].r = 2.0 * row[nby2].i ; cg[nby2].i = 0.0 ;
00336 
00337    /* phase shift both rows (cf,cg) */
00338 
00339    dk = (2.0*PI) / nup ;
00340    sf = -af * dk ; sg = -ag * dk ;
00341    for( ii=1 ; ii <= nby2 ; ii++ ){
00342       fac = CEXPIT(ii*sf) ; cf[ii] = CMULT( fac , cf[ii] ) ;
00343       gac = CEXPIT(ii*sg) ; cg[ii] = CMULT( gac , cg[ii] ) ;
00344    }
00345    cf[nby2].i = 0.0 ; cg[nby2].i = 0.0 ;
00346 
00347    /* retangle the coefficients from 2 rows */
00348 
00349    row[0].r = cf[0].r ; row[0].i = cg[0].r ;
00350    for( ii=1 ; ii < nby2 ; ii++ ){
00351       row[ii].r     =  cf[ii].r - cg[ii].i ;
00352       row[ii].i     =  cf[ii].i + cg[ii].r ;
00353       row[nup-ii].r =  cf[ii].r + cg[ii].i ;
00354       row[nup-ii].i = -cf[ii].i + cg[ii].r ;
00355    }
00356    row[nby2].r = cf[nby2].r ;
00357    row[nby2].i = cg[nby2].r ;
00358 
00359    /* inverse FFT and store back in output arrays */
00360 
00361    csfft_cox( 1 , nup , row ) ;
00362 
00363    sf = 0.5 / nup ;              /* 0.5 to allow for twice too big above */
00364    for( ii=0 ; ii < n ; ii++ ){
00365       f[ii] = sf * row[ii].r ; g[ii] = sf * row[ii].i ;
00366    }
00367 
00368    return ;
00369 }

00374 {
00375    int jj , nxup ;
00376    float * fj0 , * fj1 , * zz=NULL ;
00377    float a0 , a1 ;
00378 
00379    if( a == 0.0 && b == 0.0 ) return ;          /* nothing to do */
00380    if( nx < 2 || ny < 1 || f == NULL ) return ; /* nothing to operate on */
00381 
00382    if( ny%2 == 1 ){                               /* we work in pairs, so */
00383       zz = (float *) malloc( sizeof(float)*nx ) ; /* if not an even number */
00384       for( jj=0 ; jj < nx ; jj++ ) zz[jj] = 0.0 ; /* of rows, make an extra */
00385    }
00386 
00387    nxup = nx ;                                 /* min FFT length */
00388    jj   = 2 ; while( jj < nxup ){ jj *= 2 ; }  /* next power of 2 larger */
00389    nxup = jj ;
00390 
00391    for( jj=0 ; jj < ny ; jj+=2 ){              /* shear rows in pairs */
00392       fj0 = f + (jj*nx) ;                      /* row 0 */
00393       fj1 = (jj < ny-1) ? (fj0 + nx) : zz ;    /* row 1 */
00394       a0  = a*(jj-0.5*ny) + b ;                /* phase ramp for row 0 */
00395       a1  = a0 + a ;                           /* phase ramp for row 1 */
00396       ft_shift2( nx , nxup , a0 , fj0 , a1 , fj1 ) ;
00397    }
00398 
00399    if( zz != NULL ) free(zz) ; /* toss the trash */
00400    return ;
00401 }

00406 {
00407    int jj , nyup , ii ;
00408    float * fj0 , * fj1 ;
00409    float a0 , a1 ;
00410 
00411    if( a == 0.0 && b == 0.0 ) return ;          /* nothing to do */
00412    if( ny < 2 || nx < 1 || f == NULL ) return ; /* nothing to operate on */
00413 
00414    /* make memory for a pair of columns */
00415 
00416    fj0 = (float *) malloc( sizeof(float) * 2*ny ) ; fj1 = fj0 + ny ;
00417 
00418    nyup = ny ;                                 /* min FFT length */
00419    jj   = 2 ; while( jj < nyup ){ jj *= 2 ; }  /* next power of 2 larger */
00420    nyup = jj ;
00421 
00422    for( jj=0 ; jj < nx ; jj+=2 ){              /* shear rows in pairs */
00423 
00424       if( jj < nx-1 ){
00425          for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = f[jj+1+ii*nx]; }
00426       } else {
00427          for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = 0.0; }
00428       }
00429 
00430       a0  = a*(jj-0.5*nx) + b ;                /* phase ramp for row 0 */
00431       a1  = a0 + a ;                           /* phase ramp for row 1 */
00432       ft_shift2( ny , nyup , a0 , fj0 , a1 , fj1 ) ;
00433 
00434       if( jj < nx-1 ){
00435          for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; f[jj+1+ii*nx] = fj1[ii]; }
00436       } else {
00437          for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; }
00438       }
00439    }
00440 
00441    free(fj0) ; return ;
00442 }

00033 {
00034    float rot_dx , rot_dy , rot_cph , rot_sph , top,bot,val ;
00035    MRI_IMAGE *imfl , *newImg ;
00036    MRI_IMARR *impair ;
00037    float *far , *nar ;
00038    float xx,yy , fx,fy ;
00039    int ii,jj, nx,ny , ix,jy , ifx,jfy ;
00040    float f_jm1,f_j00,f_jp1,f_jp2 , wt_m1,wt_00,wt_p1,wt_p2 ;
00041 
00042 ENTRY("mri_affine_bicubic") ;
00043 
00044    if( im == NULL || ! MRI_IS_2D(im) ){
00045      fprintf(stderr,"*** mri_affine only works on 2D images!\n"); RETURN(NULL):
00046    }
00047 
00048    /* if complex image: break into pairs, do separately, reassemble */
00049 
00050    if( im->kind == MRI_complex ){
00051       MRI_IMARR *impair ;
00052       MRI_IMAGE * rim , * iim , * tim ;
00053       impair = mri_complex_to_pair( im ) ;
00054       if( impair == NULL ){
00055          fprintf(stderr,"*** mri_complex_to_pair fails in mri_affine!\n");EXIT(1);
00056       }
00057       rim = IMAGE_IN_IMARR(impair,0) ;
00058       iim = IMAGE_IN_IMARR(impair,1) ;  FREE_IMARR(impair) ;
00059 
00060       tim = mri_affine_bicubic( rim , a11,a12,a21,a22,xshift,yshift);
00061       mri_free(rim) ; rim = tim ;
00062 
00063       tim = mri_affine_bicubic( iim , a11,a12,a21,a22,xshift,yshift);
00064       mri_free(iim) ; iim = tim ;
00065 
00066       newImg = mri_pair_to_complex( rim , iim ) ;
00067       mri_free( rim ) ; mri_free( iim ) ;
00068       MRI_COPY_AUX(newImg,im) ;
00069       RETURN( newImg );
00070    }
00071 
00072    /** movement params **/
00073 
00074    det = a11*a22 - a21*a12 ;
00075    if( fabs(det) < 1.e-5*(fabs(a11)+fabs(a12)+fabs(a21)+fabs(a22)) ){
00076       fprintf(stderr,"*** input determinant=0 in mri_affine!\n");EXIT(1);
00077    }
00078    det = 1.0 / det ;
00079 
00080    rot_cph = cos(phi) ; rot_sph = sin(phi) ;
00081 
00082    rot_dx  = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00083             -(0.5 * im->ny) * rot_sph ;
00084 
00085    rot_dy  = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00086             +(0.5 * im->ny) * (1.0-rot_cph) ;
00087 
00088    /** other initialization **/
00089 
00090    nx = im->nx ;  /* image dimensions */
00091    ny = im->ny ;
00092 
00093    if( im->kind == MRI_float ) imfl = im ;
00094    else                        imfl = mri_to_float( im ) ;
00095 
00096    far = MRI_FLOAT_PTR(imfl) ;              /* access to float data */
00097    newImg = mri_new( nx , nx , MRI_float ) ;   /* output image */
00098    nar = MRI_FLOAT_PTR(newImg) ;               /* output image data */
00099 
00100    bot = top = far[0] ;
00101    for( ii=0 ; ii < nx*ny ; ii++ )
00102            if( far[ii] < bot ) bot = far[ii] ;
00103       else if( far[ii] > top ) top = far[ii] ;
00104 
00105    /*** loop over output points and warp to them ***/
00106 
00107    for( jj=0 ; jj < nx ; jj++ ){
00108       xx = rot_sph * jj + rot_dx - rot_cph ;
00109       yy = rot_cph * jj + rot_dy + rot_sph ;
00110       for( ii=0 ; ii < nx ; ii++ ){
00111 
00112          xx += rot_cph ;  /* get x,y in original image */
00113          yy -= rot_sph ;
00114 
00115          ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ;  /* floor */
00116          jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
00117 
00118          fx    = xx-ix ;
00119          wt_m1 = P_M1(fx) ; wt_00 = P_00(fx) ;
00120          wt_p1 = P_P1(fx) ; wt_p2 = P_P2(fx) ;
00121 
00122          if( ix > 0 && ix < nx-2 && jy > 0 && jy < ny-2 ){
00123             float * fym1, *fy00 , *fyp1 , *fyp2 ;
00124 
00125             fym1 = far + (ix-1 + (jy-1)*nx) ;
00126             fy00 = fym1 + nx ;
00127             fyp1 = fy00 + nx ;
00128             fyp2 = fyp1 + nx ;
00129 
00130             f_jm1 =  wt_m1 * fym1[0] + wt_00 * fym1[1]
00131                    + wt_p1 * fym1[2] + wt_p2 * fym1[3] ;
00132 
00133             f_j00 =  wt_m1 * fy00[0] + wt_00 * fy00[1]
00134                    + wt_p1 * fy00[2] + wt_p2 * fy00[3] ;
00135 
00136             f_jp1 =  wt_m1 * fyp1[0] + wt_00 * fyp1[1]
00137                    + wt_p1 * fyp1[2] + wt_p2 * fyp1[3] ;
00138 
00139             f_jp2 =  wt_m1 * fyp2[0] + wt_00 * fyp2[1]
00140                    + wt_p1 * fyp2[2] + wt_p2 * fyp2[3] ;
00141 
00142          } else {
00143 
00144             f_jm1 =  wt_m1 * FINS(ix-1,jy-1)
00145                    + wt_00 * FINS(ix  ,jy-1)
00146                    + wt_p1 * FINS(ix+1,jy-1)
00147                    + wt_p2 * FINS(ix+2,jy-1) ;
00148 
00149             f_j00 =   wt_m1 * FINS(ix-1,jy)
00150                     + wt_00 * FINS(ix  ,jy)
00151                     + wt_p1 * FINS(ix+1,jy)
00152                     + wt_p2 * FINS(ix+2,jy) ;
00153 
00154             f_jp1 =   wt_m1 * FINS(ix-1,jy+1)
00155                     + wt_00 * FINS(ix  ,jy+1)
00156                     + wt_p1 * FINS(ix+1,jy+1)
00157                     + wt_p2 * FINS(ix+2,jy+1) ;
00158 
00159             f_jp2 =   wt_m1 * FINS(ix-1,jy+2)
00160                     + wt_00 * FINS(ix  ,jy+2)
00161                     + wt_p1 * FINS(ix+1,jy+2)
00162                     + wt_p2 * FINS(ix+2,jy+2) ;
00163          }
00164 
00165 #define THIRTYSIX 2.7777778e-2  /* 1./36.0, actually */
00166 
00167          fy  = yy-jy ;
00168          val = (  P_M1(fy) * f_jm1 + P_00(fy) * f_j00
00169                 + P_P1(fy) * f_jp1 + P_P2(fy) * f_jp2 ) * THIRTYSIX ;
00170 
00171               if( val < bot ) nar[ii+jj*nx] = bot ;  /* too small! */
00172          else if( val > top ) nar[ii+jj*nx] = top ;  /* too big!   */
00173          else                 nar[ii+jj*nx] = val ;  /* just right */
00174 
00175       }
00176    }
00177 
00178    /*** cleanup and return ***/
00179 
00180    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00181    MRI_COPY_AUX(newImg,im) ;
00182    RETURN( newImg );
00183 }

00196 {
00197    float rot_dx , rot_dy , rot_cph , rot_sph ;
00198    MRI_IMAGE *imfl , *newImg ;
00199    MRI_IMARR *impair ;
00200    float *far , *nar ;
00201    float xx,yy , fx,fy ;
00202    int ii,jj, nx,ny , ix,jy ;
00203    float f_j00,f_jp1 , wt_00,wt_p1 ;
00204 
00205 ENTRY("mri_rota_bilinear") ;
00206 
00207    if( im == NULL || ! MRI_IS_2D(im) ){
00208      fprintf(stderr,"*** mri_rota_bilinear only works on 2D images!\n"); RETURN(NULL);
00209    }
00210 
00211    /** if complex image, break into pairs, do each separately, put back together **/
00212 
00213    if( im->kind == MRI_complex ){
00214       MRI_IMARR *impair ;
00215       MRI_IMAGE * rim , * iim , * tim ;
00216       impair = mri_complex_to_pair( im ) ;
00217       if( impair == NULL ){
00218          fprintf(stderr,"*** mri_complex_to_pair fails in mri_rota!\n") ; EXIT(1) ;
00219       }
00220       rim = IMAGE_IN_IMARR(impair,0) ;
00221       iim = IMAGE_IN_IMARR(impair,1) ;  FREE_IMARR(impair) ;
00222       tim = mri_rota_bilinear( rim , aa,bb,phi ) ; mri_free( rim ) ; rim = tim ;
00223       tim = mri_rota_bilinear( iim , aa,bb,phi ) ; mri_free( iim ) ; iim = tim ;
00224       newImg = mri_pair_to_complex( rim , iim ) ;
00225       mri_free( rim ) ; mri_free( iim ) ;
00226       MRI_COPY_AUX(newImg,im) ;
00227       RETURN( newImg );
00228    }
00229 
00230    /** rotation params **/
00231 
00232    rot_cph = cos(phi) ; rot_sph = sin(phi) ;
00233 
00234    rot_dx  = (0.5 * im->nx) * (1.0-rot_cph) - aa*rot_cph - bb*rot_sph
00235             -(0.5 * im->ny) * rot_sph ;
00236 
00237    rot_dy  = (0.5 * im->nx) * rot_sph + aa*rot_sph - bb*rot_cph
00238             +(0.5 * im->ny) * (1.0-rot_cph) ;
00239 
00240    /** other initialization **/
00241 
00242    nx = im->nx ;  /* image dimensions */
00243    ny = im->ny ;
00244 
00245    if( im->kind == MRI_float ) imfl = im ;
00246    else                        imfl = mri_to_float( im ) ;
00247 
00248    far = MRI_FLOAT_PTR(imfl) ;              /* access to float data */
00249    newImg = mri_new( nx , nx , MRI_float ) ;   /* output image */
00250    nar = MRI_FLOAT_PTR(newImg) ;               /* output image data */
00251 
00252    /*** loop over output points and warp to them ***/
00253 
00254    for( jj=0 ; jj < nx ; jj++ ){
00255       xx = rot_sph * jj + rot_dx - rot_cph ;
00256       yy = rot_cph * jj + rot_dy + rot_sph ;
00257       for( ii=0 ; ii < nx ; ii++ ){
00258 
00259          xx += rot_cph ;  /* get x,y in original image */
00260          yy -= rot_sph ;
00261 
00262          ix = (xx >= 0.0) ? ((int) xx) : ((int) xx)-1 ;  /* floor */
00263          jy = (yy >= 0.0) ? ((int) yy) : ((int) yy)-1 ;
00264 
00265          fx = xx-ix ; wt_00 = 1.0 - fx ; wt_p1 = fx ;
00266 
00267          if( ix >= 0 && ix < nx-1 && jy >= 0 && jy < ny-1 ){
00268             float *fy00 , *fyp1 ;
00269 
00270             fy00 = far + (ix + jy*nx) ; fyp1 = fy00 + nx ;
00271 
00272             f_j00 = wt_00 * fy00[0] + wt_p1 * fy00[1] ;
00273             f_jp1 = wt_00 * fyp1[0] + wt_p1 * fyp1[1] ;
00274 
00275          } else {
00276             f_j00 = wt_00 * FINS(ix,jy  ) + wt_p1 * FINS(ix+1,jy  ) ;
00277             f_jp1 = wt_00 * FINS(ix,jy+1) + wt_p1 * FINS(ix+1,jy+1) ;
00278          }
00279 
00280          fy  = yy-jy ; nar[ii+jj*nx] = (1.0-fy) * f_j00 + fy * f_jp1 ;
00281 
00282       }
00283    }
00284 
00285    /*** cleanup and return ***/
00286 
00287    if( im != imfl ) mri_free(imfl) ;  /* throw away unneeded workspace */
00288    MRI_COPY_AUX(newImg,im) ;
00289    RETURN( newImg );
00290 }

00447 {
00448    double cph , sph ;
00449    float a , b , bot,top ;
00450    MRI_IMAGE *flim ;
00451    float *flar ;
00452    int ii , nxy ;
00453 
00454 ENTRY("mri_rota_shear") ;
00455 
00456    if( im == NULL || ! MRI_IS_2D(im) ){
00457      fprintf(stderr,"*** mri_rota_shear only works on 2D images!\n"); RETURN(NULL);
00458    }
00459 
00460    /** if complex image, break into pairs, do each separately, put back together **/
00461 
00462    if( im->kind == MRI_complex ){
00463       MRI_IMARR *impair ;
00464       MRI_IMAGE * rim , * iim , * tim ;
00465       impair = mri_complex_to_pair( im ) ;
00466       if( impair == NULL ){
00467          fprintf(stderr,"*** mri_complex_to_pair fails in mri_rota!\n") ; EXIT(1) ;
00468       }
00469       rim  = IMAGE_IN_IMARR(impair,0) ;
00470       iim  = IMAGE_IN_IMARR(impair,1) ;  FREE_IMARR(impair) ;
00471       tim  = mri_rota_shear( rim , aa,bb,phi ) ; mri_free( rim ) ; rim = tim ;
00472       tim  = mri_rota_shear( iim , aa,bb,phi ) ; mri_free( iim ) ; iim = tim ;
00473       flim = mri_pair_to_complex( rim , iim ) ;
00474       mri_free( rim ) ; mri_free( iim ) ;
00475       MRI_COPY_AUX(flim,im) ;
00476       RETURN( flim );
00477    }
00478 
00479    /** copy input to output **/
00480 
00481    flim = mri_to_float( im ) ;
00482    flar = MRI_FLOAT_PTR( flim ) ;
00483 
00484    /* find range of image data */
00485 
00486    bot = top = flar[0] ; nxy = im->nx * im->ny ;
00487    for( ii=1 ; ii < nxy ; ii++ )
00488            if( flar[ii] < bot ) bot = flar[ii] ;
00489       else if( flar[ii] > top ) top = flar[ii] ;
00490 
00491    /** rotation params **/
00492 
00493    cph = cos(phi) ; sph = sin(phi) ;
00494 
00495    /* More than 90 degrees?
00496       Must be reduced to less than 90 degrees by a 180 degree flip. */
00497 
00498    if( cph < 0.0 ){
00499       int ii , jj , top , nx=flim->nx , ny=flim->ny ;
00500       float val ;
00501 
00502       top = (nx+1)/2 ;
00503       for( jj=0 ; jj < ny ; jj++ ){
00504          for( ii=1 ; ii < top ; ii++ ){
00505             val               = flar[jj*nx+ii] ;
00506             flar[jj*nx+ii]    = flar[jj*nx+nx-ii] ;
00507             flar[jj*nx+nx-ii] = val ;
00508          }
00509       }
00510 
00511       top = (ny+1)/2 ;
00512       for( ii=0 ; ii < nx ; ii++ ){
00513          for( jj=1 ; jj < top ; jj++ ){
00514             val                 = flar[ii+jj*nx] ;
00515             flar[ii+jj*nx]      = flar[ii+(ny-jj)*nx] ;
00516             flar[ii+(ny-jj)*nx] = val ;
00517          }
00518       }
00519 
00520       cph = -cph ; sph = -sph ;
00521    }
00522 
00523    /* compute shear factors for each direction */
00524 
00525    b = sph ;
00526    a = (b != 0.0 ) ? ((cph - 1.0) / b) : (0.0) ;
00527 
00528    /* shear thrice */
00529 
00530    ft_xshear( a , 0.0       , im->nx , im->ny , flar ) ;
00531    ft_yshear( b , bb        , im->nx , im->ny , flar ) ;
00532    ft_xshear( a , aa - a*bb , im->nx , im->ny , flar ) ;
00533 
00534    /* make sure data does not go out of original range */
00535 
00536    for( ii=0 ; ii < nxy ; ii++ )
00537            if( flar[ii] < bot ) flar[ii] = bot ;
00538       else if( flar[ii] > top ) flar[ii] = top ;
00539 
00540    RETURN( flim );
00541 }

00544 {
00545    switch(mode){
00546 
00547       default:
00548       case MRI_BICUBIC:  return mri_rota( im,aa,bb,phi ) ;
00549 
00550       case MRI_BILINEAR: return mri_rota_bilinear( im,aa,bb,phi ) ;
00551 
00552       case MRI_FOURIER:  return mri_rota_shear( im,aa,bb,phi ) ;
00553    }
00554 }