00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <mrilib.h>
00023 #include <string.h>
00024
00025 #define NMAX 1050
00026
00027
00028
00029
00030
00031
00032
00033 void detrend( int n , float vec[] , float wt[] )
00034 {
00035 register int ii ;
00036 float sum0 , sum1 , sum2 , det , cf,lf ;
00037
00038 static int first = 1 ;
00039 static float cf0,cf1 , lf0,lf1 ;
00040
00041
00042
00043 if( first ){
00044 first = 0 ;
00045 sum0 = sum1 = sum2 = 0.0 ;
00046 for( ii=0 ; ii < n ; ii++ ){
00047 sum0 += wt[ii] ;
00048 sum1 += wt[ii] * ii ;
00049 sum2 += wt[ii] * ii*ii ;
00050 }
00051 det = sum0 * sum2 - sum1 * sum1 ;
00052 cf0 = sum2 / det ;
00053 cf1 = -sum1 / det ;
00054 lf0 = cf1 ;
00055 lf1 = sum0 / det ;
00056 }
00057
00058
00059
00060 sum0 = sum1 = 0.0 ;
00061 for( ii=0 ; ii < n ; ii++ ){
00062 sum0 += wt[ii] * vec[ii] ;
00063 sum1 += wt[ii] * vec[ii] * ii ;
00064 }
00065
00066 cf = cf0 * sum0 + cf1 * sum1 ;
00067 lf = lf0 * sum0 + lf1 * sum1 ;
00068 for( ii=0 ; ii < n ; ii++ ) vec[ii] -= cf + ii*lf ;
00069
00070 return ;
00071 }
00072
00073
00074
00075
00076
00077
00078
00079 void normalize( int n , float vec[] , float wt[] )
00080 {
00081 register int ii ;
00082 float sqsum ;
00083
00084 detrend( n , vec , wt ) ;
00085
00086 sqsum = 0.0 ;
00087 for( ii=0 ; ii < n ; ii++ ) sqsum += wt[ii] * vec[ii] * vec[ii] ;
00088
00089 if( sqsum < 0.000001 ){
00090 for( ii=0 ; ii < n ; ii++ ) vec[ii] = 0.0 ;
00091 } else {
00092 sqsum = 1.0 / sqrt(sqsum) ;
00093 for( ii=0 ; ii < n ; ii++ ) vec[ii] *= sqsum ;
00094 }
00095
00096 return ;
00097 }
00098
00099
00100
00101
00102
00103 void main( int argc , char *argv[] )
00104 {
00105 MRI_IMAGE *inim[NMAX] , *tempim ;
00106 short *inar[NMAX] ;
00107 int ii,jj,joff,ijoff , kk , narg , nimage , nx , ny , nimax=NMAX ;
00108 float cthresh = 0.5 ;
00109 float scale , fmax,fmin , sum ;
00110
00111 FILE *reffile , *covfile , *savefile=NULL ;
00112 float ref[NMAX] , vec[NMAX] , org[NMAX] ;
00113 float wt[NMAX] ;
00114 int lref = 0 , lcov = 0 ;
00115 char buf[128] , cfname[128] , sfname[128] ;
00116 char *cchh ;
00117
00118 double *vsum , *vsqsum ;
00119 int numv , mvlen , lx,my , moff ;
00120 int cvtype = 2 , cvinp ;
00121
00122 int ldetrend = FALSE , lmedfilt = FALSE , lsave = FALSE ;
00123
00124
00125
00126 narg = 1 ;
00127 kk = 0 ;
00128 do {
00129
00130
00131
00132 if( argc < 4 || strncmp(argv[narg],"-help",2) == 0 ){
00133 fprintf( stderr ,
00134 "Usage: %s [-nim number] [-pcnt threshold] [-detrend] [-noquad] \n" ,
00135 argv[0] ) ;
00136 fprintf( stderr ,
00137 " [-medfilt] [-save savefile] cov_file ref_file im_files\n" );
00138 exit(0) ;
00139 }
00140
00141 if( strncmp(argv[narg],"-save",4) == 0 ){
00142 strcpy( sfname , argv[++narg] ) ;
00143 lsave = TRUE ;
00144 continue ;
00145 }
00146
00147 if( strncmp(argv[narg],"-detrend",4) == 0 ){
00148 ldetrend = TRUE ;
00149 continue ;
00150 }
00151
00152 if( strncmp(argv[narg],"-noquad",4) == 0 ){
00153 cvtype = 1 ;
00154 continue ;
00155 }
00156
00157 if( strncmp(argv[narg],"-medfilt",4) == 0 ){
00158 lmedfilt = TRUE ;
00159 continue ;
00160 }
00161
00162 if( strncmp(argv[narg],"-pcnt",4) == 0 ){
00163 cthresh = strtod( argv[++narg] , NULL ) ;
00164 if( cthresh > 1 ) cthresh /= 100.0 ;
00165 continue ;
00166 }
00167
00168 if( strncmp(argv[narg],"-nim",4) == 0 ){
00169 nimax = strtol( argv[++narg] , NULL , 0 ) ;
00170 if( nimax > NMAX || nimax < 9 ) nimax = NMAX ;
00171 continue ;
00172 }
00173
00174 if( strncmp(argv[narg],"-",1) == 0 ){
00175 fprintf( stderr , "unknown switch %s\n" , argv[narg] ) ;
00176 exit(1) ;
00177 }
00178
00179
00180
00181 if( !lcov ){
00182 covfile = fopen( argv[narg] , "r" ) ;
00183 strcpy( cfname , argv[narg] ) ;
00184 if( covfile == NULL ){
00185 lcov = -1 ;
00186 } else {
00187 lcov = 1 ;
00188 }
00189 continue ;
00190 }
00191
00192 if( !lref ){
00193 if( cthresh <= 0.0 ){
00194 fprintf( stderr , "skipping ref file %s\n" , argv[narg] ) ;
00195 reffile == NULL ;
00196 } else {
00197 reffile = fopen( argv[narg] , "r" ) ;
00198 if( reffile == NULL ){
00199 fprintf( stderr , "cannot open ref file %s\n" , argv[narg] ) ;
00200 exit(1) ;
00201 }
00202 }
00203 lref = 1 ;
00204 continue ;
00205 }
00206
00207 tempim = mri_read( argv[narg] ) ;
00208
00209 if( tempim == NULL ) exit(1) ;
00210 if( kk >= nimax ){
00211 fprintf( stderr , "max #files exceeded at file # %d\n" , kk+1 ) ;
00212 exit(1) ;
00213 }
00214
00215 if( kk == 0 ){
00216
00217 nx = tempim->nx ;
00218 ny = tempim->ny ;
00219
00220 fmin = mri_min( tempim ) ;
00221 fmax = mri_max( tempim ) ;
00222 if( fmin >= fmax ){
00223 fprintf( stderr , "1st image is constant!\n" ) ;
00224 exit(1) ;
00225 }
00226 scale = 10000.0 / (fmax-fmin) ;
00227
00228 } else if( nx != tempim->nx || ny != tempim->ny ){
00229 fprintf( stderr ,
00230 "file shape doesn't match first image: %s\n" , argv[narg] ) ;
00231 exit(1) ;
00232 }
00233
00234 if( tempim->kind == MRI_short ){
00235 inim[kk] = tempim ;
00236 } else {
00237 inim[kk] = mri_to_short( scale , tempim ) ;
00238 mri_free( tempim ) ;
00239 }
00240 inar[kk] = mri_data_pointer( inim[kk] ) ;
00241 kk++ ;
00242 continue ;
00243
00244 } while( ++narg < argc && kk < nimax ) ;
00245
00246
00247
00248 if( !lcov || !lref || kk < 9 ){
00249 fprintf( stderr , "enough files not given!\n" ) ;
00250 exit(1) ;
00251 }
00252
00253 nimage = kk ;
00254
00255
00256
00257 if( cthresh > 0.0 ){
00258 for( kk=0 ; kk < nimage ; kk++ ){
00259 cchh = fgets( buf , 100 , reffile ) ;
00260 if( cchh == NULL ){
00261 fprintf( stderr , "premature EOF in ref file at line # %d\n" , kk+1 ) ;
00262 exit(1) ;
00263 }
00264 ref[kk] = strtod( buf , NULL ) ;
00265 }
00266 fclose( reffile ) ;
00267 } else {
00268 for( kk=0 ; kk < nimage ; kk++ ) ref[kk] = kk ;
00269 }
00270
00271 for( kk=0 ; kk < nimage ; kk++ ){
00272 wt[kk] = (ref[kk] < 33333) ? 1.0 : 0.0 ;
00273 }
00274
00275 normalize( nimage , ref , wt ) ;
00276
00277
00278
00279 vsum = (double *) malloc( sizeof(double) * nimage ) ;
00280 if( vsum == NULL ){
00281 fprintf( stderr , "cannot malloc space for vector sum\n" ) ;
00282 exit(1) ;
00283 }
00284
00285 if( cvtype >= 2 ){
00286 vsqsum = (double *) malloc( sizeof(double) * nimage*nimage ) ;
00287 if( vsqsum == NULL ){
00288 fprintf( stderr ,"cannot malloc space for cov matrix\n" ) ;
00289 exit(1) ;
00290 }
00291 }
00292
00293
00294
00295 #define COVERR(nn) if(ii<(nn)){fprintf(stderr,"error reading cov file\n");exit(1);}
00296
00297 if( lcov == 1 ){
00298 ii = fread( &numv , sizeof(int) , 1 , covfile ) ; COVERR(1) ;
00299 ii = fread( &mvlen , sizeof(int) , 1 , covfile ) ; COVERR(1) ;
00300 ii = fread( &cvinp , sizeof(int) , 1 , covfile ) ; COVERR(1) ;
00301
00302 if( mvlen != nimage || numv < 0 ){
00303 fprintf( stderr , "cov file has wrong sizes: nv mv = %d %d\n" ,
00304 numv , mvlen ) ;
00305 exit(1) ;
00306 }
00307 if( cvinp != cvtype ){
00308 fprintf( stderr , "cov file has wrong type: %d\n" , cvinp ) ;
00309 exit(1) ;
00310 }
00311
00312 ii = fread( vsum , sizeof(double), nimage , covfile ) ;
00313 COVERR(nimage);
00314 if( cvtype >= 2 ){
00315 ii = fread( vsqsum, sizeof(double), nimage*nimage, covfile ) ;
00316 COVERR(nimage*nimage) ;
00317 }
00318
00319 fclose(covfile) ;
00320 } else if( lcov == -1 ){
00321 mvlen = nimage ;
00322 numv = 0 ;
00323 for( ii=0 ; ii < nimage ; ii++ ) vsum[ii] = 0 ;
00324 if( cvtype >= 2 ){
00325 for( ii=0 ; ii < nimage*nimage ; ii++ ) vsqsum[ii] = 0 ;
00326 }
00327 } else {
00328 fprintf( stderr , "illegal value of lcov occured!\n" ) ;
00329 exit(1) ;
00330 }
00331
00332
00333
00334 for( jj=0 ; jj < ny ; jj++ ){
00335 joff = jj * nx ;
00336 for( ii=0 ; ii < nx ; ii++ ){
00337 ijoff = ii + joff ;
00338
00339
00340
00341 for( kk=0 ; kk < nimage ; kk++ ) vec[kk] = inar[kk][ijoff] ;
00342
00343
00344
00345 if( lmedfilt ){
00346 org[0] = vec[0] ;
00347 org[nimage-1] = vec[nimage-1] ;
00348 for( kk=0 ; kk < nimage-1 ; kk++ )
00349 org[kk] = MEDIAN( vec[kk-1] , vec[kk] , vec[kk+1] ) ;
00350 } else {
00351 for( kk=0 ; kk < nimage ; kk++ ) org[kk] = vec[kk] ;
00352 }
00353
00354
00355
00356 if( cthresh > 0.0 ){
00357 for( kk=0 ; kk < nimage ; kk++ ) vec[kk] = org[kk] ;
00358 normalize( nimage , vec , wt ) ;
00359 sum = 0.0 ;
00360 for( kk=0 ; kk < nimage ; kk++ ) sum += wt[kk] * ref[kk] * vec[kk] ;
00361 } else {
00362 sum = 0.0 ;
00363 }
00364
00365
00366
00367 if( sum >= cthresh ){
00368
00369 numv++ ;
00370
00371 if( ldetrend ) detrend( nimage , org , wt ) ;
00372
00373 for( my=0 ; my < nimage ; my++ ){
00374 vsum[my] += org[my] ;
00375 if( cvtype >= 2 ){
00376 moff = nimage * my ;
00377 for( lx=0 ; lx < nimage ; lx++ )
00378 vsqsum[lx+moff] += org[my]*org[lx] ;
00379 }
00380 }
00381
00382 if( lsave && cthresh > 0.0 ){
00383 if( savefile == NULL ){
00384 savefile = fopen( sfname , "a" ) ;
00385 if( savefile == NULL ){
00386 fprintf( stderr , "cannot open save file\n" ) ;
00387 exit(1) ;
00388 }
00389 }
00390 fprintf( savefile , "# x=%3d y=%3d\n" , ii,jj) ;
00391 for( lx=0 ; lx < nimage ; lx++ )
00392 fprintf( savefile , "%d %12.4e\n" , lx,org[lx] ) ;
00393 }
00394
00395 }
00396
00397 }
00398 }
00399
00400
00401
00402 covfile = fopen( cfname , "w" ) ;
00403 fwrite( &numv , sizeof(int) , 1 , covfile ) ;
00404 fwrite( &mvlen , sizeof(int) , 1 , covfile ) ;
00405 fwrite( &cvtype, sizeof(int) , 1 , covfile ) ;
00406 fwrite( vsum , sizeof(double) , nimage , covfile ) ;
00407 if( cvtype >= 2 ){
00408 fwrite( vsqsum , sizeof(double) , nimage*nimage , covfile ) ;
00409 }
00410 fclose( covfile ) ;
00411
00412 printf( "# vectors in covariance file now %d\n" , numv ) ;
00413
00414 exit(0) ;
00415 }