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 /*--------------------------------------------------------------------------
00012  Routine to filter a 2D image in Fourier space:
00013   float sigma --> if > 0, sigma for Gaussian blur: exp(-k**2 * sigma**2/2)
00014   int   diffx --> if != 0, apply d/dx:             i * k_x
00015   int   diffy --> if != 0, apply d/dy:             i * k_y
00016   int   code  --> bit mapped options to control filter.  currently:
00017                     FILT_FFT_WRAPAROUND  -->  allows wraparound in FFT
00018 
00019  Note that the units of sigma are the same as those of im->dx, im->dy
00020  (defaults are 1.0 = units of pixel dimensions).
00021 
00022  The output is an image in the floating point format.  Input can be
00023  any format (but complex inputs will have CABS taken first).
00024 ----------------------------------------------------------------------------*/
00025 
00026 #define GET_AS_BIG(name,type,dim)                                           \
00027    do{ if( (dim) > name ## _size ){                                         \
00028           if( name != NULL ) free(name) ;                                   \
00029           name = (type *) malloc( sizeof(type) * (dim) ) ;                  \
00030           if( name == NULL ){                                               \
00031              fprintf(stderr,"*** can't malloc mri_filt_fft workspace\n") ;  \
00032              EXIT(1) ; }                                                    \
00033           name ## _size = (dim) ; }                                         \
00034        break ; } while(1)
00035 
00036 MRI_IMAGE * mri_filt_fft( MRI_IMAGE * im ,
00037                           float sigma , int diffx , int diffy , int code )
00038 {
00039    int jj, nby2 , nx,ny ;
00040    float  dk , aa , k , fac , dx,dy ;
00041    register int ii , nup ;
00042    register float * bfim ;
00043 
00044    static int cx_size  = 0 ;     /* workspaces (will hang around between calls) */
00045    static int gg_size  = 0 ;
00046    static complex * cx = NULL ;
00047    static float   * gg = NULL ;
00048 
00049    MRI_IMAGE * flim ;
00050    float     * flar ;
00051 
00052    /*** initialize ***/
00053 
00054    if( im == NULL ){
00055       fprintf(stderr,"*** mri_filt_fft: NULL input image\n") ;
00056       return NULL ;
00057    }
00058 
00059    if( sigma < 0.0 ){
00060       fprintf(stderr,"*** mri_filt_fft: Illegal control parameters input\n");
00061       return NULL ;
00062    }
00063 
00064    if( ! MRI_IS_2D(im) ){
00065       fprintf(stderr,"*** mri_filt_fft: Only works on 2D images\n") ;
00066       EXIT(1) ;
00067    }
00068 
00069    nx = im->nx ; ny = im->ny ;
00070    dx = fabs(im->dx) ; if( dx == 0.0 ) dx = 1.0 ;
00071    dy = fabs(im->dy) ; if( dy == 0.0 ) dy = 1.0 ;
00072 
00073    aa = sigma * sigma * 0.5 ;
00074 
00075    flim = mri_to_float( im ) ;        /* will be output */
00076    flar = mri_data_pointer( flim ) ;
00077 
00078    if( sigma == 0.0 && diffx == 0 && diffy == 0 ) return flim ;  /* no action! */
00079 
00080    /*** do x-direction ***/
00081 
00082    if( (code & FILT_FFT_WRAPAROUND) == 0 ){
00083       nup = nx + (int)(3.0 * sigma / dx) ;      /* min FFT length */
00084    } else {
00085       nup = nx ;
00086    }
00087    ii  = 4 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00088    nup = ii ; nby2 = nup / 2 ;
00089 
00090    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00091 
00092    dk    = (2.0*PI) / (nup * dx) ;
00093    fac   = 1.0 / nup ;
00094    gg[0] = fac ;
00095    if( aa > 0.0 ){
00096       for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00097    } else {
00098       for( ii=1 ; ii < nup ; ii++ ) gg[ii] = fac ;
00099    }
00100 
00101    if( diffx ){
00102       gg[0] = gg[nby2] = 0.0 ;
00103       for( ii=1 ; ii < nby2 ; ii++ ){ k=ii*dk ; gg[ii] *= k ; gg[nup-ii] *= (-k) ; }
00104    }
00105 
00106    /** July 19 1995: double up on FFTs **/
00107 
00108    for( jj=0 ; jj < ny ; jj+=2 ){
00109       bfim = flar + jj*nx ;
00110       if( jj == ny-1 )
00111          for( ii=0 ; ii<nx ; ii++){ cx[ii].r = bfim[ii] ; cx[ii].i = 0.0 ; }  /* copy in */
00112       else
00113          for( ii=0 ; ii<nx ; ii++){ cx[ii].r = bfim[ii] ; cx[ii].i = bfim[ii+nx] ; }
00114       for( ii=nx; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }                /* zero pad */
00115       csfft_cox( -1 , nup , cx ) ;                                            /* FFT */
00116       for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }  /* filter */
00117       if( diffx ){
00118          float tt ;
00119          for( ii=0 ; ii < nup ; ii++ ){
00120             tt = cx[ii].r ; cx[ii].r = -cx[ii].i ; cx[ii].i = tt ;            /* mult by i */
00121          }
00122       }
00123       csfft_cox(  1 , nup , cx ) ;                                            /* inv FFT */
00124       if( jj == ny-1 )
00125          for( ii=0 ; ii<nx ; ii++){ bfim[ii] = cx[ii].r ; }                   /* copy out */
00126       else
00127          for( ii=0 ; ii<nx ; ii++){ bfim[ii] = cx[ii].r ; bfim[ii+nx] = cx[ii].i ; }
00128    }
00129 
00130    /*** do y-direction ***/
00131 
00132    if( (code & FILT_FFT_WRAPAROUND) == 0 ){
00133       nup = ny + (int)(3.0 * sigma / dy) ;      /* min FFT length */
00134    } else {
00135       nup = ny ;
00136    }
00137    ii  = 2 ; while( ii < nup ){ ii *= 2 ; }  /* next power of 2 larger */
00138    nup = ii ; nby2 = nup / 2 ;
00139 
00140    GET_AS_BIG(cx,complex,nup) ; GET_AS_BIG(gg,float,nup) ;
00141 
00142    dk    = (2.0*PI) / (nup * dy) ;
00143    fac   = 1.0 / nup ;
00144    gg[0] = fac ;
00145 
00146    if( aa > 0.0 ){
00147       for( ii=1 ; ii<=nby2 ; ii++ ){ k=ii*dk; gg[nup-ii]=gg[ii]=fac*exp(-aa*k*k); }
00148    } else {
00149       for( ii=1 ; ii < nup ; ii++ ) gg[ii] = fac ;
00150    }
00151 
00152    if( diffy ){
00153       gg[0] = gg[nby2] = 0.0 ;
00154       for( ii=1 ; ii < nby2 ; ii++ ){ k=ii*dk ; gg[ii] *= k ; gg[nup-ii] *= (-k) ; }
00155    }
00156 
00157    for( jj=0 ; jj < nx ; jj+=2 ){
00158       bfim = flar + jj ;
00159       if( jj == nx-1 )
00160          for( ii=0 ; ii<ny ; ii++){ cx[ii].r = bfim[ii*nx] ; cx[ii].i = 0.0 ; }
00161       else
00162          for( ii=0 ; ii<ny ; ii++){ cx[ii].r = bfim[ii*nx] ; cx[ii].i = bfim[ii*nx+1] ; }
00163       for( ii=ny; ii<nup; ii++){ cx[ii].r = cx[ii].i = 0.0 ; }
00164       csfft_cox( -1 , nup , cx ) ;
00165       for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
00166       if( diffy ){
00167          float tt ;
00168          for( ii=0 ; ii < nup ; ii++ ){
00169             tt = cx[ii].r ; cx[ii].r = -cx[ii].i ; cx[ii].i = tt ;  /* multiply by i */
00170          }
00171       }
00172       csfft_cox(  1 , nup , cx ) ;
00173       if( jj == nx-1 )
00174          for( ii=0 ; ii<ny ; ii++){ bfim[ii*nx] = cx[ii].r ; }
00175       else
00176          for( ii=0 ; ii<ny ; ii++){ bfim[ii*nx] = cx[ii].r ; bfim[ii*nx+1] = cx[ii].i ; }
00177    }
00178 
00179    /*** done! ***/
00180 
00181    return flim ;
00182 }