00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 #include <string.h>
00008 #include "mrilib.h"
00009 
00010 #define IMMAX 1024
00011 #define NFMAX 512   /* half of IMMAX */
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    /*** deal with command line arguments ***/
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    /*** switches ***/
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    /** read and initialize images **/
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 ){                        /* check size for compatibility */
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 ){  /* convert to floating point */
00088          tempim = mri_to_float( inim ) ;
00089          mri_free( inim ) ;
00090          IMAGE_IN_IMARR(inimar,kim) = inim = tempim ;
00091       }
00092    }
00093 
00094    /*** cut image count back to power of 2 ***/
00095 
00096    for( ii=2 ; ii <= nimage ; ii *= 2 ) ; /* null loop */
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    /*** load array of pointers to data arrays,
00110         create array of output images ,
00111         load array of pointers to output arrays ,
00112         create taper array                        ***/
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    /*** foreach time series:
00136            remove mean, taper, FFT, square, store in output ***/
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    /*** toss input ***/
00162 
00163    DESTROY_IMARR( inimar ) ;
00164    free( taper ) ;
00165 
00166    /*** scale to shorts and write to disk ***/
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 ;  /* 50 dB range */
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 }