00001 /*** program COVMAT:
00002      inputs:  (optional) covariance matrix file,
00003               reference time series file,
00004               a sequence of image files.
00005      output:  (updated) covariance matrix file, with all vectors
00006                 correlated above threshold "added" into covariance matrix,
00007               (optional) file with all such vectors saved for display
00008 
00009      Uses libmri.a and mrilib.h (MRI file library routines).
00010 
00011      RW Cox, December 1993
00012 ***/
00013 
00014 /*** modified January 1994 to include FIM-like "ignore" locations in
00015      reference time series.  The mechanism is to use a weighted inner
00016      product everywhere -- now, the weights are 0 or 1, but can easily
00017      be modified later.
00018 ***/
00019 
00020 /****************************************************************************/
00021 
00022 #include <mrilib.h>
00023 #include <string.h>
00024 
00025 #define NMAX 1050   /* maximum length of time series of images */
00026 
00027 /***************************************************************************/
00028 
00029 /*** detrend:
00030      routine to remove unwanted components from time series
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 ;             /* initialization flag */
00039    static float cf0,cf1 , lf0,lf1 ;   /* to be initialized */
00040 
00041  /*** initialize coefficients for detrending ***/
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 ;     /* constant factor for sum0 */
00053       cf1 = -sum1 / det ;     /* constant factor for sum1 */
00054       lf0 = cf1 ;             /* linear factor for sum0 */
00055       lf1 =  sum0 / det ;     /* linear factor for sum1 */
00056    }
00057 
00058  /*** remove mean and linear trend ***/
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 /*** normalize:
00076      routine to scale a time series to unit vector
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 /*** the main enchilada ***/
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] ;                           /* Jan 1994 addition */
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 ;               /* Jan 1994 addition */
00121 
00122    int ldetrend = FALSE , lmedfilt = FALSE , lsave = FALSE ;
00123 
00124  /*** inputs ***/
00125 
00126    narg = 1 ;   /* argument counter */
00127    kk   = 0 ;   /* image counter */
00128    do {
00129 
00130     /*** deal with switches ***/
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     /*** deal with files ***/
00180 
00181       if( !lcov ){                               /* 1st file = cov matrix */
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 ){                            /* 2nd file = ref time series */
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] ) ;    /* rest of files = images */
00208 
00209       if( tempim == NULL ) exit(1) ;       /* check for errors on read */
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 ){           /* 1st image file ==> do some initializations */
00216 
00217          nx = tempim->nx ;     /* save dimensions */
00218          ny = tempim->ny ;
00219 
00220          fmin = mri_min( tempim ) ;   /* compute a scale factor */
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 ){  /* check dimensions */
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 ){   /* convert all inputs to shorts */
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] ) ;   /* pointer to data */
00241       kk++ ;
00242       continue ;
00243 
00244    } while( ++narg < argc && kk < nimax ) ;  /* end of loop over arguments */
00245 
00246  /*** check for reasonable inputs at this point */
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  /*** read ref file now ***/
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++ ){            /* set weights */
00272       wt[kk] = (ref[kk] < 33333) ? 1.0 : 0.0 ;
00273    }
00274 
00275    normalize( nimage , ref , wt ) ;   /* make into a unit vector */
00276 
00277  /*** allocate space for cov file stuff ***/
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  /*** read cov file now, if present ***/
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  /*** do all pixels ***/
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        /* load data from (ii,jj) pixel into vec[] */
00340 
00341          for( kk=0 ; kk < nimage ; kk++ ) vec[kk] = inar[kk][ijoff] ;
00342 
00343        /* if desired, median filter into org[] */
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        /* dot product with reference */
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        /* if we still like this pixel, do stuff with it */
00366 
00367          if( sum >= cthresh ){
00368 
00369             numv++ ;  /* another vector has passed */
00370 
00371             if( ldetrend ) detrend( nimage , org , wt ) ;
00372 
00373             for( my=0 ; my < nimage ; my++ ){  /* form sum & sum-of-products */
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 ){      /* save timeseries? */
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          } /* end sum >= cthresh */
00396 
00397       }  /* end ii */
00398    }  /* end jj */
00399 
00400  /*** save covariance ***/
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 }