Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
imfft.c
Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007 #include <string.h>
00008 #include "mrilib.h"
00009
00010 #define IMMAX 1024
00011 #define NFMAX 512
00012
00013 int main( int argc , char *argv[] )
00014 {
00015 char prefix[64] = "fft." , fname[128] ;
00016 int narg , ii , nx,ny,npix , nimage,nfreq , kim ;
00017
00018 MRI_IMARR *inimar , *outimar ;
00019
00020 MRI_IMAGE *tempim , *inim , *outim ;
00021 float **inar , **outar , *taper ;
00022 complex *dat ;
00023 float sum , scale , pbot,ptop ;
00024 short *tempar ;
00025 int ldecibel = FALSE ;
00026
00027
00028
00029 if( argc < 3 || strncmp(argv[1],"-help",5) == 0 ){
00030 printf( "Computation of |FFT| of image time series.\n"
00031 "Usage: imfft [-prefix prefix] image_files\n" ) ;
00032 exit(0) ;
00033 }
00034
00035 machdep() ;
00036
00037 narg = 1 ;
00038 kim = 0 ;
00039
00040
00041
00042 while( argv[narg][0] == '-' ){
00043
00044 if( strncmp(argv[narg],"-prefix",3) == 0 ){
00045 strncpy( prefix , argv[++narg] , 64 ) ;
00046 ii = strlen(prefix) ;
00047 if( prefix[ii] != '.' ){ prefix[ii+1] = '.' ; prefix[ii+2] = '\0' ; }
00048 narg++ ; continue ;
00049 }
00050
00051 if( strncmp(argv[narg],"-",1) == 0 ){
00052 fprintf( stderr , "*** illegal switch: %s\a\n" , argv[narg] ) ;
00053 exit(1) ;
00054 }
00055 }
00056
00057
00058
00059 inimar = mri_read_many_files( argc-narg , argv+narg ) ;
00060 if( inimar == NULL ){
00061 fprintf(stderr,"*** no input images read!\a\n") ;
00062 exit(1) ;
00063 } else if( inimar->num < 32 ){
00064 fprintf(stderr,"*** less than 32 input images read!\a\n") ;
00065 exit(1) ;
00066 }
00067
00068 nimage = inimar->num ;
00069 for( kim=0 ; kim < nimage; kim++ ){
00070
00071 inim = IMAGE_IN_IMARR(inimar,kim) ;
00072
00073 if( ! MRI_IS_2D(inim) ){
00074 fprintf(stderr,"*** can only process 2D images!\a\n") ;
00075 exit(1) ;
00076 }
00077
00078 if( kim == 0 ){
00079 nx = inim->nx ;
00080 ny = inim->ny ;
00081 } else if( inim->nx != nx || inim->ny != ny ){
00082 fprintf( stderr ,
00083 "image %d size doesn't match 1st image!\n" , kim+1 ) ;
00084 exit(1) ;
00085 }
00086
00087 if( inim->kind != MRI_float ){
00088 tempim = mri_to_float( inim ) ;
00089 mri_free( inim ) ;
00090 IMAGE_IN_IMARR(inimar,kim) = inim = tempim ;
00091 }
00092 }
00093
00094
00095
00096 for( ii=2 ; ii <= nimage ; ii *= 2 ) ;
00097 kim = nimage ;
00098 nimage = ii / 2 ;
00099 nfreq = nimage/2 - 1 ;
00100 if( nimage < kim ){
00101 fprintf( stderr , "cutting image count back to %d from %d\n" ,
00102 nimage , kim ) ;
00103 for( ii=nimage ; ii < kim ; ii++ ){
00104 mri_free( IMAGE_IN_IMARR(inimar,ii) ) ;
00105 IMAGE_IN_IMARR(inimar,ii) = NULL ;
00106 }
00107 }
00108
00109
00110
00111
00112
00113
00114 inar = (float **) malloc( sizeof(float *) * nimage ) ;
00115 outar = (float **) malloc( sizeof(float *) * nfreq ) ;
00116 dat = (complex *) malloc( sizeof(complex) * nimage ) ;
00117
00118 if( inar==NULL || outar==NULL || dat==NULL ){
00119 fprintf(stderr,"*** cannot malloc workspace for FFT!\a\n") ;
00120 exit(1) ;
00121 }
00122
00123 for( kim=0 ; kim < nimage ; kim++ )
00124 inar[kim] = MRI_FLOAT_PTR( IMAGE_IN_IMARR(inimar,kim) ) ;
00125
00126 INIT_IMARR(outimar) ;
00127 for( kim=0 ; kim < nfreq ; kim++ ){
00128 outim = mri_new( nx , ny , MRI_float ) ;
00129 outar[kim] = MRI_FLOAT_PTR( outim ) ;
00130 ADDTO_IMARR(outimar,outim) ;
00131 }
00132
00133 taper = mri_setup_taper( nimage , 0.20 ) ;
00134
00135
00136
00137
00138 npix = nx * ny ;
00139 for( ii=0 ; ii < npix ; ii++ ){
00140 sum = 0.0 ;
00141 for( kim=0 ; kim < nimage ; kim++ ) sum += inar[kim][ii] ;
00142 sum /= nimage ;
00143
00144 for( kim=0 ; kim < nimage ; kim++ ){
00145 dat[kim].r = (inar[kim][ii] - sum) * taper[kim] ;
00146 dat[kim].i = 0.0 ;
00147 }
00148 csfft_cox( -1 , nimage , dat ) ;
00149
00150 for( kim=0 ; kim < nfreq ; kim++ )
00151 outar[kim][ii] = CABS(dat[kim+1]) ;
00152
00153 #if 0
00154 if( ldecibel ){
00155 for( kim=0 ; kim < nfreq ; kim++ )
00156 outar[kim][ii] = 10.0 * log10( outar[kim][ii] + 1.e-30 ) ;
00157 }
00158 #endif
00159 }
00160
00161
00162
00163 DESTROY_IMARR( inimar ) ;
00164 free( taper ) ;
00165
00166
00167
00168 ptop = pbot = outar[0][0] ;
00169 for( kim=0 ; kim < nfreq ; kim++ ){
00170 sum = mri_max( IMAGE_IN_IMARR(outimar,kim) ) ; if( sum > ptop ) ptop = sum ;
00171 sum = mri_min( IMAGE_IN_IMARR(outimar,kim) ) ; if( sum < pbot ) pbot = sum ;
00172 }
00173
00174 fprintf( stderr , "\n minimum = %g maximum = %g\n" , pbot,ptop ) ;
00175
00176 #if 0
00177 if( ldecibel ){
00178 pbot = ptop - 50.0 ;
00179 } else {
00180 pbot = 0.0 ;
00181 }
00182 #else
00183 pbot = 0.0 ;
00184 #endif
00185
00186 scale = 30000.0 / (ptop-pbot) ;
00187 tempim = mri_new( nx , ny , MRI_short ) ;
00188 tempar = mri_data_pointer( tempim ) ;
00189
00190 for( kim=0 ; kim < nfreq ; kim++ ){
00191
00192 for( ii=0 ; ii < npix ; ii++ ){
00193 tempar[ii] = (outar[kim][ii] < pbot)
00194 ? 0
00195 : (short)(scale * (outar[kim][ii] - pbot) + 0.499) ;
00196 }
00197
00198 mri_free( IMAGE_IN_IMARR(outimar,kim) ) ;
00199
00200 sprintf( fname , "%s%04d" , prefix , kim ) ;
00201 mri_write( fname , tempim ) ;
00202
00203 }
00204
00205 exit(0) ;
00206 }