do{ if( (dim) > name ## _size ){                                         \
          if( name != NULL ) free(name) ;                                   \
          name = (type *) malloc( sizeof(type) * (dim) ) ;                  \
          if( name == NULL ){                                               \
             fprintf(stderr,"*** can't malloc mri_filt_fft workspace\n") ;  \
             EXIT(1) ; }                                                    \
          name ## _size = (dim) ; }                                         \
       break ; } while(1)

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 }