1 /*****************************************************************************
2 Major portions of this software are copyrighted by the Medical College
3 of Wisconsin, 1994-2000, and are released under the Gnu General Public
4 License, Version 2. See the file README.Copyright for details.
5 ******************************************************************************/
6
7 #include <string.h>
8 #include "mrilib.h"
9
10 #define IMMAX 1024
11 #define NFMAX 512 /* half of IMMAX */
12
main(int argc,char * argv[])13 int main( int argc , char *argv[] )
14 {
15 char prefix[64] = "fft." , fname[128] ;
16 int narg , ii , nx,ny,npix , nimage,nfreq , kim ;
17
18 MRI_IMARR *inimar , *outimar ;
19
20 MRI_IMAGE *tempim , *inim , *outim ;
21 float **inar , **outar , *taper ;
22 complex *dat ;
23 float sum , scale , pbot,ptop ;
24 short *tempar ;
25 int ldecibel = FALSE ;
26
27 /*** deal with command line arguments ***/
28
29 if( argc < 3 || strncmp(argv[1],"-help",5) == 0 ){
30 printf( "Computation of |FFT| of image time series.\n"
31 "Usage: imfft [-prefix prefix] image_files\n" ) ;
32 exit(0) ;
33 }
34
35 machdep() ;
36
37 narg = 1 ;
38 kim = 0 ;
39
40 /*** switches ***/
41
42 while( argv[narg][0] == '-' ){
43
44 if( strncmp(argv[narg],"-prefix",3) == 0 ){
45 strncpy( prefix , argv[++narg] , 64 ) ;
46 ii = strlen(prefix) ;
47 if( prefix[ii] != '.' ){ prefix[ii+1] = '.' ; prefix[ii+2] = '\0' ; }
48 narg++ ; continue ;
49 }
50
51 if( strncmp(argv[narg],"-",1) == 0 ){
52 fprintf( stderr , "*** illegal switch: %s\a\n" , argv[narg] ) ;
53 exit(1) ;
54 }
55 }
56
57 /** read and initialize images **/
58
59 inimar = mri_read_many_files( argc-narg , argv+narg ) ;
60 if( inimar == NULL ){
61 fprintf(stderr,"*** no input images read!\a\n") ;
62 exit(1) ;
63 } else if( inimar->num < 32 ){
64 fprintf(stderr,"*** less than 32 input images read!\a\n") ;
65 exit(1) ;
66 }
67
68 nimage = inimar->num ;
69 for( kim=0 ; kim < nimage; kim++ ){
70
71 inim = IMAGE_IN_IMARR(inimar,kim) ;
72
73 if( ! MRI_IS_2D(inim) ){
74 fprintf(stderr,"*** can only process 2D images!\a\n") ;
75 exit(1) ;
76 }
77
78 if( kim == 0 ){ /* check size for compatibility */
79 nx = inim->nx ;
80 ny = inim->ny ;
81 } else if( inim->nx != nx || inim->ny != ny ){
82 fprintf( stderr ,
83 "image %d size doesn't match 1st image!\n" , kim+1 ) ;
84 exit(1) ;
85 }
86
87 if( inim->kind != MRI_float ){ /* convert to floating point */
88 tempim = mri_to_float( inim ) ;
89 mri_free( inim ) ;
90 IMAGE_IN_IMARR(inimar,kim) = inim = tempim ;
91 }
92 }
93
94 /*** cut image count back to power of 2 ***/
95
96 for( ii=2 ; ii <= nimage ; ii *= 2 ) ; /* null loop */
97 kim = nimage ;
98 nimage = ii / 2 ;
99 nfreq = nimage/2 - 1 ;
100 if( nimage < kim ){
101 fprintf( stderr , "cutting image count back to %d from %d\n" ,
102 nimage , kim ) ;
103 for( ii=nimage ; ii < kim ; ii++ ){
104 mri_free( IMAGE_IN_IMARR(inimar,ii) ) ;
105 IMAGE_IN_IMARR(inimar,ii) = NULL ;
106 }
107 }
108
109 /*** load array of pointers to data arrays,
110 create array of output images ,
111 load array of pointers to output arrays ,
112 create taper array ***/
113
114 inar = (float **) malloc( sizeof(float *) * nimage ) ;
115 outar = (float **) malloc( sizeof(float *) * nfreq ) ;
116 dat = (complex *) malloc( sizeof(complex) * nimage ) ;
117
118 if( inar==NULL || outar==NULL || dat==NULL ){
119 fprintf(stderr,"*** cannot malloc workspace for FFT!\a\n") ;
120 exit(1) ;
121 }
122
123 for( kim=0 ; kim < nimage ; kim++ )
124 inar[kim] = MRI_FLOAT_PTR( IMAGE_IN_IMARR(inimar,kim) ) ;
125
126 INIT_IMARR(outimar) ;
127 for( kim=0 ; kim < nfreq ; kim++ ){
128 outim = mri_new( nx , ny , MRI_float ) ;
129 outar[kim] = MRI_FLOAT_PTR( outim ) ;
130 ADDTO_IMARR(outimar,outim) ;
131 }
132
133 taper = mri_setup_taper( nimage , 0.20 ) ;
134
135 /*** foreach time series:
136 remove mean, taper, FFT, square, store in output ***/
137
138 npix = nx * ny ;
139 for( ii=0 ; ii < npix ; ii++ ){
140 sum = 0.0 ;
141 for( kim=0 ; kim < nimage ; kim++ ) sum += inar[kim][ii] ;
142 sum /= nimage ;
143
144 for( kim=0 ; kim < nimage ; kim++ ){
145 dat[kim].r = (inar[kim][ii] - sum) * taper[kim] ;
146 dat[kim].i = 0.0 ;
147 }
148 csfft_cox( -1 , nimage , dat ) ;
149
150 for( kim=0 ; kim < nfreq ; kim++ )
151 outar[kim][ii] = CABS(dat[kim+1]) ;
152
153 #if 0
154 if( ldecibel ){
155 for( kim=0 ; kim < nfreq ; kim++ )
156 outar[kim][ii] = 10.0 * log10( outar[kim][ii] + 1.e-30 ) ;
157 }
158 #endif
159 }
160
161 /*** toss input ***/
162
163 DESTROY_IMARR( inimar ) ;
164 free( taper ) ;
165
166 /*** scale to shorts and write to disk ***/
167
168 ptop = pbot = outar[0][0] ;
169 for( kim=0 ; kim < nfreq ; kim++ ){
170 sum = mri_max( IMAGE_IN_IMARR(outimar,kim) ) ; if( sum > ptop ) ptop = sum ;
171 sum = mri_min( IMAGE_IN_IMARR(outimar,kim) ) ; if( sum < pbot ) pbot = sum ;
172 }
173
174 fprintf( stderr , "\n minimum = %g maximum = %g\n" , pbot,ptop ) ;
175
176 #if 0
177 if( ldecibel ){
178 pbot = ptop - 50.0 ; /* 50 dB range */
179 } else {
180 pbot = 0.0 ;
181 }
182 #else
183 pbot = 0.0 ;
184 #endif
185
186 scale = 30000.0 / (ptop-pbot) ;
187 tempim = mri_new( nx , ny , MRI_short ) ;
188 tempar = mri_data_pointer( tempim ) ;
189
190 for( kim=0 ; kim < nfreq ; kim++ ){
191
192 for( ii=0 ; ii < npix ; ii++ ){
193 tempar[ii] = (outar[kim][ii] < pbot)
194 ? 0
195 : (short)(scale * (outar[kim][ii] - pbot) + 0.499) ;
196 }
197
198 mri_free( IMAGE_IN_IMARR(outimar,kim) ) ;
199
200 sprintf( fname , "%s%04d" , prefix , kim ) ;
201 mri_write( fname , tempim ) ;
202
203 }
204
205 exit(0) ;
206 }
207