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 } |