Doxygen Source Code Documentation
mri_filt_fft.c File Reference
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | GET_AS_BIG(name, type, dim) |
Functions | |
MRI_IMAGE * | mri_filt_fft (MRI_IMAGE *im, float sigma, int diffx, int diffy, int code) |
Define Documentation
|
Value: 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) Definition at line 26 of file mri_filt_fft.c. |
Function Documentation
|
Definition at line 36 of file mri_filt_fft.c. References csfft_cox(), MRI_IMAGE::dx, MRI_IMAGE::dy, EXIT, FILT_FFT_WRAPAROUND, GET_AS_BIG, complex::i, mri_data_pointer(), MRI_IS_2D, mri_to_float(), MRI_IMAGE::nx, MRI_IMAGE::ny, complex::r, and tt. Referenced by main(), mri_2dalign_one(), mri_2dalign_setup(), mri_align_crao(), and mri_align_dfspace().
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 } |