Doxygen Source Code Documentation
Main Page Alphabetical List Data Structures File List Data Fields Globals Search
impsd.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 256
00012
00013 int main( int argc , char *argv[] )
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
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
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
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 ){
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 ){
00094 tempim = mri_to_float( inim ) ;
00095 mri_free( inim ) ;
00096 IMAGE_IN_IMARR(inimar,kim) = inim = tempim ;
00097 }
00098 }
00099
00100
00101
00102 for( ii=2 ; ii <= nimage ; ii *= 2 ) ;
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
00116
00117
00118
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
00142
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
00166
00167 DESTROY_IMARR( inimar ) ;
00168 free( taper ) ;
00169
00170
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 ;
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 }