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

00320 {
00321    static int nupold=0 ;
00322    static complex * row=NULL , * cf=NULL , * cg=NULL ;
00323 
00324    int ii , nby2=nup/2 , n21=nby2+1 ;
00325    complex fac , gac ;
00326    float sf , sg , dk ;
00327 
00328    /* make new memory for row storage? */
00329 
00330    if( nup > nupold ){
00331       if( row != NULL ){ free(row) ; free(cf) ; free(cg) ; }
00332       row = (complex *) malloc( sizeof(complex) * nup ) ;
00333       cf  = (complex *) malloc( sizeof(complex) * n21 ) ;
00334       cg  = (complex *) malloc( sizeof(complex) * n21 ) ;
00335       nupold = nup ;
00336    }
00337 
00338    /* FFT the pair of rows */
00339 
00340    for( ii=0 ; ii < n   ; ii++ ){ row[ii].r = f[ii] ; row[ii].i = g[ii] ; }
00341    for(      ; ii < nup ; ii++ ){ row[ii].r = row[ii].i = 0.0 ; }
00342 
00343    csfft_cox( -1 , nup , row ) ;
00344 
00345    /* untangle FFT coefficients from row into cf,cg */
00346 
00347    cf[0].r = 2.0 * row[0].r ; cf[0].i = 0.0 ;  /* twice too big */
00348    cg[0].r = 2.0 * row[0].i ; cg[0].i = 0.0 ;
00349    for( ii=1 ; ii < nby2 ; ii++ ){
00350       cf[ii].r =  row[ii].r + row[nup-ii].r ;
00351       cf[ii].i =  row[ii].i - row[nup-ii].i ;
00352       cg[ii].r =  row[ii].i + row[nup-ii].i ;
00353       cg[ii].i = -row[ii].r + row[nup-ii].r ;
00354    }
00355    cf[nby2].r = 2.0 * row[nby2].r ; cf[nby2].i = 0.0 ;
00356    cg[nby2].r = 2.0 * row[nby2].i ; cg[nby2].i = 0.0 ;
00357 
00358    /* phase shift both rows (cf,cg) */
00359 
00360    dk = (2.0*PI) / nup ;
00361    sf = -af * dk ; sg = -ag * dk ;
00362    for( ii=1 ; ii <= nby2 ; ii++ ){
00363       fac = CEXPIT(ii*sf) ; cf[ii] = CMULT( fac , cf[ii] ) ;
00364       gac = CEXPIT(ii*sg) ; cg[ii] = CMULT( gac , cg[ii] ) ;
00365    }
00366    cf[nby2].i = 0.0 ; cg[nby2].i = 0.0 ;
00367 
00368    /* retangle the coefficients from 2 rows */
00369 
00370    row[0].r = cf[0].r ; row[0].i = cg[0].r ;
00371    for( ii=1 ; ii < nby2 ; ii++ ){
00372       row[ii].r     =  cf[ii].r - cg[ii].i ;
00373       row[ii].i     =  cf[ii].i + cg[ii].r ;
00374       row[nup-ii].r =  cf[ii].r + cg[ii].i ;
00375       row[nup-ii].i = -cf[ii].i + cg[ii].r ;
00376    }
00377    row[nby2].r = cf[nby2].r ;
00378    row[nby2].i = cg[nby2].r ;
00379 
00380    /* inverse FFT and store back in output arrays */
00381 
00382    csfft_cox( 1 , nup , row ) ;
00383 
00384    sf = 0.5 / nup ;              /* 0.5 to allow for twice too big above */
00385    for( ii=0 ; ii < n ; ii++ ){
00386       f[ii] = sf * row[ii].r ; g[ii] = sf * row[ii].i ;
00387    }
00388 
00389    return ;
00390 }

00395 {
00396    int jj , nxup ;
00397    float * fj0 , * fj1 , * zz=NULL ;
00398    float a0 , a1 ;
00399 
00400    if( a == 0.0 && b == 0.0 ) return ;          /* nothing to do */
00401    if( nx < 2 || ny < 1 || f == NULL ) return ; /* nothing to operate on */
00402 
00403    if( ny%2 == 1 ){                               /* we work in pairs, so */
00404       zz = (float *) malloc( sizeof(float)*nx ) ; /* if not an even number */
00405       for( jj=0 ; jj < nx ; jj++ ) zz[jj] = 0.0 ; /* of rows, make an extra */
00406    }
00407 
00408    nxup = nx ;                                 /* min FFT length */
00409    jj   = 2 ; while( jj < nxup ){ jj *= 2 ; }  /* next power of 2 larger */
00410    nxup = jj ;
00411 
00412    for( jj=0 ; jj < ny ; jj+=2 ){              /* shear rows in pairs */
00413       fj0 = f + (jj*nx) ;                      /* row 0 */
00414       fj1 = (jj < ny-1) ? (fj0 + nx) : zz ;    /* row 1 */
00415       a0  = a*(jj-0.5*ny) + b ;                /* phase ramp for row 0 */
00416       a1  = a0 + a ;                           /* phase ramp for row 1 */
00417       ft_shift2( nx , nxup , a0 , fj0 , a1 , fj1 ) ;
00418    }
00419 
00420    if( zz != NULL ) free(zz) ; /* toss the trash */
00421    return ;
00422 }

00427 {
00428    int jj , nyup , ii ;
00429    float * fj0 , * fj1 ;
00430    float a0 , a1 ;
00431 
00432    if( a == 0.0 && b == 0.0 ) return ;          /* nothing to do */
00433    if( ny < 2 || nx < 1 || f == NULL ) return ; /* nothing to operate on */
00434 
00435    /* make memory for a pair of columns */
00436 
00437    fj0 = (float *) malloc( sizeof(float) * 2*ny ) ; fj1 = fj0 + ny ;
00438 
00439    nyup = ny ;                                 /* min FFT length */
00440    jj   = 2 ; while( jj < nyup ){ jj *= 2 ; }  /* next power of 2 larger */
00441    nyup = jj ;
00442 
00443    for( jj=0 ; jj < nx ; jj+=2 ){              /* shear rows in pairs */
00444 
00445       if( jj < nx-1 ){
00446          for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = f[jj+1+ii*nx]; }
00447       } else {
00448          for( ii=0; ii < ny; ii++ ){ fj0[ii] = f[jj+ii*nx]; fj1[ii] = 0.0; }
00449       }
00450 
00451       a0  = a*(jj-0.5*nx) + b ;                /* phase ramp for row 0 */
00452       a1  = a0 + a ;                           /* phase ramp for row 1 */
00453       ft_shift2( ny , nyup , a0 , fj0 , a1 , fj1 ) ;
00454 
00455       if( jj < nx-1 ){
00456          for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; f[jj+1+ii*nx] = fj1[ii]; }
00457       } else {
00458          for( ii=0; ii < ny; ii++ ){ f[jj+ii*nx] = fj0[ii]; }
00459       }
00460    }
00461 
00462    free(fj0) ; return ;
00463 }

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

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

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

00563 {
00564    switch(mode){
00565 
00566       default:
00567       case MRI_BICUBIC:  return mri_rota( im,aa,bb,phi ) ;
00568 
00569       case MRI_BILINEAR: return mri_rota_bilinear( im,aa,bb,phi ) ;
00570 
00571       case MRI_FOURIER:  return mri_rota_shear( im,aa,bb,phi ) ;
00572    }
00573 }