Doxygen Source Code Documentation
covmat.c File Reference
#include <mrilib.h>#include <string.h>Go to the source code of this file.
Defines | |
| #define | NMAX 1050 |
| #define | COVERR(nn) if(ii<(nn)){fprintf(stderr,"error reading cov file\n");exit(1);} |
Functions | |
| void | detrend (int n, float vec[], float wt[]) |
| void | normalize (int n, float vec[], float wt[]) |
| void | main (int argc, char *argv[]) |
Define Documentation
|
|
|
|
|
Definition at line 25 of file covmat.c. Referenced by main(). |
Function Documentation
|
||||||||||||||||
|
Definition at line 33 of file covmat.c. References vec.
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 }
|
|
||||||||||||
|
convert DTIStudio fiber format data to SUMA segment data Definition at line 103 of file covmat.c. References argc, detrend(), MRI_IMAGE::kind, malloc, MEDIAN, mri_data_pointer(), mri_free(), mri_max(), mri_min(), mri_read(), mri_to_short(), NMAX, normalize(), MRI_IMAGE::nx, MRI_IMAGE::ny, ref, scale, strtod(), and vec.
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 }
|
|
||||||||||||||||
|
Definition at line 79 of file covmat.c. References detrend(), and vec.
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 }
|
