00014 {
00015    char prefix[64] = "psd." , 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 PSD of image time series.\n"
00031               "Copyright 1994,1995 Medical College of Wisconsin.\n\n"
00032               "Usage: impsd [-prefix prefix] [-dB] image_files\n" ) ;
00033       exit(0) ;
00034    }
00035 
00036    machdep() ;
00037 
00038    narg = 1 ;
00039    kim  = 0 ;
00040 
00041    /*** switches ***/
00042 
00043    while( argv[narg][0] == '-' ){
00044 
00045       if( strncmp(argv[narg],"-prefix",3) == 0 ){
00046          strncpy( prefix , argv[++narg] , 64 ) ;
00047          ii = strlen(prefix) ;
00048          if( prefix[ii] != '.' ){ prefix[ii+1] = '.' ; prefix[ii+2] = '\0' ; }
00049          narg++ ; continue ;
00050       }
00051 
00052       if( strncmp(argv[narg],"-dB",3) == 0 ){
00053          ldecibel = TRUE ;
00054          narg++ ; continue ;
00055       }
00056 
00057       if( strncmp(argv[narg],"-",1) == 0 ){
00058          fprintf( stderr , "*** illegal switch: %s\a\n" , argv[narg] ) ;
00059          exit(1) ;
00060       }
00061    }
00062 
00063    /** read and initialize images **/
00064 
00065    inimar = mri_read_many_files( argc-narg , argv+narg ) ;
00066    if( inimar == NULL ){
00067       fprintf(stderr,"*** no input images read!\a\n") ;
00068       exit(1) ;
00069    } else if( inimar->num < 32 ){
00070       fprintf(stderr,"*** less than 32 input images read!\a\n") ;
00071       exit(1) ;
00072    }
00073 
00074    nimage = inimar->num ;
00075    for( kim=0 ; kim < nimage; kim++ ){
00076 
00077       inim = IMAGE_IN_IMARR(inimar,kim) ;
00078 
00079       if( ! MRI_IS_2D(inim) ){
00080          fprintf(stderr,"*** can only process 2D images!\a\n") ;
00081          exit(1) ;
00082       }
00083 
00084       if( kim == 0 ){                        /* check size for compatibility */
00085          nx = inim->nx ;
00086          ny = inim->ny ;
00087       } else if( inim->nx != nx || inim->ny != ny ){
00088          fprintf( stderr ,
00089                   "image %d size doesn't match 1st image!\n" , kim+1 ) ;
00090          exit(1) ;
00091       }
00092 
00093       if( inim->kind != MRI_float ){  /* convert to floating point */
00094          tempim = mri_to_float( inim ) ;
00095          mri_free( inim ) ;
00096          IMAGE_IN_IMARR(inimar,kim) = inim = tempim ;
00097       }
00098    }
00099 
00100    /*** cut image count back to power of 2 ***/
00101 
00102    for( ii=2 ; ii <= nimage ; ii *= 2 ) ; /* null loop */
00103    kim    = nimage ;
00104    nimage = ii / 2 ;
00105    nfreq  = nimage/2 - 1 ;
00106    if( nimage < kim ){
00107       fprintf( stderr , "cutting image count back to %d from %d\n" ,
00108                nimage , kim ) ;
00109       for( ii=nimage ; ii < kim ; ii++ ){
00110          mri_free( IMAGE_IN_IMARR(inimar,ii) ) ;
00111          IMAGE_IN_IMARR(inimar,ii) = NULL ;
00112       }
00113    }
00114 
00115    /*** load array of pointers to data arrays,
00116         create array of output images ,
00117         load array of pointers to output arrays ,
00118         create taper array                        ***/
00119 
00120    inar  = (float **)  malloc( sizeof(float *) * nimage ) ;
00121    outar = (float **)  malloc( sizeof(float *) * nfreq  ) ;
00122    dat   = (complex *) malloc( sizeof(complex) * nimage ) ;
00123 
00124    if( inar==NULL || outar==NULL || dat==NULL ){
00125       fprintf(stderr,"*** cannot malloc workspace for PSD!\a\n") ;
00126       exit(1) ;
00127    }
00128 
00129    for( kim=0 ; kim < nimage ; kim++ )
00130       inar[kim] = MRI_FLOAT_PTR( IMAGE_IN_IMARR(inimar,kim) ) ;
00131 
00132    INIT_IMARR(outimar) ;
00133    for( kim=0 ; kim < nfreq ; kim++ ){
00134       outim      = mri_new( nx , ny , MRI_float ) ;
00135       outar[kim] = MRI_FLOAT_PTR( outim ) ;
00136       ADDTO_IMARR(outimar,outim) ;
00137    }
00138 
00139    taper = mri_setup_taper( nimage , 0.50 ) ;
00140 
00141    /*** foreach time series:
00142            remove mean, taper, FFT, square, store in output ***/
00143 
00144    npix = nx * ny ;
00145    for( ii=0 ; ii < npix ; ii++ ){
00146       sum = 0.0 ;
00147       for( kim=0 ; kim < nimage ; kim++ ) sum += inar[kim][ii] ;
00148       sum /= nimage ;
00149 
00150       for( kim=0 ; kim < nimage ; kim++ ){
00151          dat[kim].r = taper[kim]*(inar[kim][ii] - sum) ;
00152          dat[kim].i = 0.0 ;
00153       }
00154       csfft_cox( -1 , nimage , dat ) ;
00155 
00156       for( kim=0 ; kim < nfreq ; kim++ )
00157          outar[kim][ii] = CSQR(dat[kim+1]) ;
00158 
00159       if( ldecibel ){
00160          for( kim=0 ; kim < nfreq ; kim++ )
00161             outar[kim][ii] = 10.0 * log10( outar[kim][ii] + 1.e-30 ) ;
00162       }
00163    }
00164 
00165    /*** toss input ***/
00166 
00167    DESTROY_IMARR( inimar ) ;
00168    free( taper ) ;
00169 
00170    /*** scale to shorts and write to disk ***/
00171 
00172    ptop = pbot = outar[0][0] ;
00173    for( kim=0 ; kim < nfreq ; kim++ ){
00174       sum = mri_max( IMAGE_IN_IMARR(outimar,kim) ) ; if( sum > ptop ) ptop = sum ;
00175       sum = mri_min( IMAGE_IN_IMARR(outimar,kim) ) ; if( sum < pbot ) pbot = sum ;
00176    }
00177 
00178    fprintf( stderr , "/n minimum = %e  maximum = %e\n" , pbot,ptop ) ;
00179 
00180    if( ldecibel ){
00181       pbot = ptop - 50.0 ;  /* 50 dB range */
00182    } else {
00183       pbot = 0.0 ;
00184    }
00185 
00186    scale  = 10000.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 }