Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
mri_filt_fft.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
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 ;
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
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 ) ;
00076 flar = mri_data_pointer( flim ) ;
00077
00078 if( sigma == 0.0 && diffx == 0 && diffy == 0 ) return flim ;
00079
00080
00081
00082 if( (code & FILT_FFT_WRAPAROUND) == 0 ){
00083 nup = nx + (int)(3.0 * sigma / dx) ;
00084 } else {
00085 nup = nx ;
00086 }
00087 ii = 4 ; while( ii < nup ){ ii *= 2 ; }
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
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 ; }
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 ; }
00115 csfft_cox( -1 , nup , cx ) ;
00116 for( ii=0 ; ii<nup; ii++){ cx[ii].r *= gg[ii] ; cx[ii].i *= gg[ii] ; }
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 ;
00121 }
00122 }
00123 csfft_cox( 1 , nup , cx ) ;
00124 if( jj == ny-1 )
00125 for( ii=0 ; ii<nx ; ii++){ bfim[ii] = cx[ii].r ; }
00126 else
00127 for( ii=0 ; ii<nx ; ii++){ bfim[ii] = cx[ii].r ; bfim[ii+nx] = cx[ii].i ; }
00128 }
00129
00130
00131
00132 if( (code & FILT_FFT_WRAPAROUND) == 0 ){
00133 nup = ny + (int)(3.0 * sigma / dy) ;
00134 } else {
00135 nup = ny ;
00136 }
00137 ii = 2 ; while( ii < nup ){ ii *= 2 ; }
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 ;
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
00180
00181 return flim ;
00182 }