00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009
00010
00011 #define METH_MEAN 0
00012 #define METH_SLOPE 1
00013 #define METH_SIGMA 2
00014 #define METH_CVAR 3
00015
00016 #define METH_MEDIAN 4
00017 #define METH_MAD 5
00018
00019 #define METH_MAX 6
00020 #define METH_MIN 7
00021
00022 #define METH_DW 8
00023
00024 #define METH_SIGMA_NOD 9
00025 #define METH_CVAR_NOD 10
00026
00027 #define METH_AUTOCORR 11
00028 #define METH_AUTOREGP 12
00029
00030 #define METH_ABSMAX 13
00031
00032 #define METH_ARGMAX 14
00033 #define METH_ARGMIN 15
00034 #define METH_ARGABSMAX 16
00035
00036 #define MAX_NUM_OF_METHS 20
00037 static int meth[MAX_NUM_OF_METHS] = {METH_MEAN};
00038 static int nmeths = 0;
00039 static char prefix[THD_MAX_PREFIX] = "stat" ;
00040 static int datum = MRI_float ;
00041 static char meth_names[][20] = {"Mean","Slope","Std Dev","Coef of Var","Median",
00042 "Med Abs Dev", "Max", "Min", "Durbin-Watson", "Std Dev (NOD)",
00043 "Coef Var(NOD)","AutoCorr","AutoReg","Absolute Max",
00044 "ArgMax","ArgMin","ArgAbsMax"};
00045 static void STATS_tsfunc( double tzero , double tdelta ,
00046 int npts , float ts[] , double ts_mean ,
00047 double ts_slope , void * ud , int nbriks, float * val ) ;
00048
00049 static void autocorr( int npts, float ints[], int numVals, float outcoeff[] ) ;
00050
00051 int main( int argc , char * argv[] )
00052 {
00053 THD_3dim_dataset * old_dset , * new_dset ;
00054 int nopt, nbriks, ii ;
00055 int addBriks = 0;
00056 int numMultBriks,methIndex,brikIndex;
00057
00058
00059 if( argc < 2 || strcmp(argv[1],"-help") == 0 ){
00060 printf("Usage: 3dTstat [options] dataset\n"
00061 "Computes one or more voxel-wise statistics for a 3D+time dataset\n"
00062 "and stores them in a bucket dataset.\n"
00063 "\n"
00064 "Options:\n"
00065 " -mean = compute mean of input voxels [DEFAULT]\n"
00066 " -slope = compute mean slope of input voxels vs. time\n"
00067 " -stdev = compute standard deviation of input voxels\n"
00068 " [N.B.: this is computed after ]\n"
00069 " [ the slope has been removed]\n"
00070 " -cvar = compute coefficient of variation of input\n"
00071 " voxels = stdev/fabs(mean)\n"
00072 " **N.B.: You can add NOD to the end of the above 2\n"
00073 " options to turn off detrending, as in\n"
00074 " -stdevNOD or -cvarNOD\n"
00075 "\n"
00076 " -MAD = compute MAD (median absolute deviation) of\n"
00077 " input voxels = median(|voxel-median(voxel)|)\n"
00078 " [N.B.: the trend is NOT removed for this]\n"
00079 " -DW = compute Durbin-Watson Statistic of\n"
00080 " input voxels\n"
00081 " [N.B.: the trend is removed for this]\n"
00082 " -median = compute median of input voxels [undetrended]\n"
00083 " -min = compute minimum of input voxels [undetrended]\n"
00084 " -max = compute maximum of input voxels [undetrended]\n"
00085 " -absmax = compute absolute maximum of input voxels [undetrended]\n"
00086 " -argmin = index of minimum of input voxels [undetrended]\n"
00087 " -argmax = index of maximum of input voxels [undetrended]\n"
00088 " -argabsmax = index of absolute maximum of input voxels [undetrended]\n"
00089 "\n"
00090 " -prefix p = use string 'p' for the prefix of the\n"
00091 " output dataset [DEFAULT = 'stat']\n"
00092 " -datum d = use data type 'd' for the type of storage\n"
00093 " of the output, where 'd' is one of\n"
00094 " 'byte', 'short', or 'float' [DEFAULT=float]\n"
00095 " -autocorr n = compute autocorrelation function and return\n"
00096 " first n coefficients\n"
00097 " -autoreg n = compute autoregression coefficients and return\n"
00098 " first n coefficients\n"
00099 " [N.B.: -autocorr 0 and/or -autoreg 0 will return coefficients\n"
00100 " equal to the length of the input data\n"
00101 "\n"
00102 "The output is a bucket dataset. The input dataset\n"
00103 "may use a sub-brick selection list, as in program 3dcalc.\n"
00104 ) ;
00105 printf("\n" MASTER_SHORTHELP_STRING ) ;
00106 exit(0) ;
00107 }
00108
00109 mainENTRY("3dTstat main"); machdep(); AFNI_logger("3dTstat",argc,argv);
00110 PRINT_VERSION("3dTstat");
00111
00112 nopt = 1 ;
00113 nbriks = 0 ;
00114 nmeths = 0 ;
00115 while( nopt < argc && argv[nopt][0] == '-' ){
00116
00117
00118
00119 if( strcmp(argv[nopt],"-median") == 0 ){
00120 meth[nmeths++] = METH_MEDIAN ;
00121 nbriks++ ;
00122 nopt++ ; continue ;
00123 }
00124
00125 if( strcmp(argv[nopt],"-DW") == 0 ){
00126 meth[nmeths++] = METH_DW ;
00127 nbriks++ ;
00128 nopt++ ; continue ;
00129 }
00130
00131 if( strcmp(argv[nopt],"-MAD") == 0 ){
00132 meth[nmeths++] = METH_MAD ;
00133 nbriks++ ;
00134 nopt++ ; continue ;
00135 }
00136
00137 if( strcmp(argv[nopt],"-mean") == 0 ){
00138 meth[nmeths++] = METH_MEAN ;
00139 nbriks++ ;
00140 nopt++ ; continue ;
00141 }
00142
00143 if( strcmp(argv[nopt],"-slope") == 0 ){
00144 meth[nmeths++] = METH_SLOPE ;
00145 nbriks++ ;
00146 nopt++ ; continue ;
00147 }
00148
00149 if( strcmp(argv[nopt],"-stdev") == 0 ||
00150 strcmp(argv[nopt],"-sigma") == 0 ){
00151
00152 meth[nmeths++] = METH_SIGMA ;
00153 nbriks++ ;
00154 nopt++ ; continue ;
00155 }
00156
00157 if( strcmp(argv[nopt],"-cvar") == 0 ){
00158 meth[nmeths++] = METH_CVAR ;
00159 nbriks++ ;
00160 nopt++ ; continue ;
00161 }
00162
00163 if( strcmp(argv[nopt],"-stdevNOD") == 0 ||
00164 strcmp(argv[nopt],"-sigmaNOD") == 0 ){
00165
00166 meth[nmeths++] = METH_SIGMA_NOD ;
00167 nbriks++ ;
00168 nopt++ ; continue ;
00169 }
00170
00171 if( strcmp(argv[nopt],"-cvarNOD") == 0 ){
00172 meth[nmeths++] = METH_CVAR_NOD ;
00173 nbriks++ ;
00174 nopt++ ; continue ;
00175 }
00176
00177 if( strcmp(argv[nopt],"-min") == 0 ){
00178 meth[nmeths++] = METH_MIN ;
00179 nbriks++ ;
00180 nopt++ ; continue ;
00181 }
00182
00183 if( strcmp(argv[nopt],"-max") == 0 ){
00184 meth[nmeths++] = METH_MAX ;
00185 nbriks++ ;
00186 nopt++ ; continue ;
00187 }
00188
00189 if( strcmp(argv[nopt],"-absmax") == 0 ){
00190 meth[nmeths++] = METH_ABSMAX ;
00191 nbriks++ ;
00192 nopt++ ; continue ;
00193 }
00194
00195 if( strcmp(argv[nopt],"-argmin") == 0 ){
00196 meth[nmeths++] = METH_ARGMIN ;
00197 nbriks++ ;
00198 nopt++ ; continue ;
00199 }
00200
00201 if( strcmp(argv[nopt],"-argmax") == 0 ){
00202 meth[nmeths++] = METH_ARGMAX ;
00203 nbriks++ ;
00204 nopt++ ; continue ;
00205 }
00206
00207 if( strcmp(argv[nopt],"-argabsmax") == 0 ){
00208 meth[nmeths++] = METH_ARGABSMAX ;
00209 nbriks++ ;
00210 nopt++ ; continue ;
00211 }
00212
00213 if( strcmp(argv[nopt],"-autocorr") == 0 ){
00214 meth[nmeths++] = METH_AUTOCORR ;
00215 if( ++nopt >= argc ){
00216 fprintf(stderr,"*** -autocorr needs an argument!\n"); exit(1);
00217 }
00218 meth[nmeths++] = atoi(argv[nopt++]);
00219 if (meth[nmeths - 1] == 0) {
00220 addBriks++;
00221 } else {
00222 nbriks+=meth[nmeths - 1] ;
00223 }
00224 continue ;
00225 }
00226
00227 if( strcmp(argv[nopt],"-autoreg") == 0 ){
00228 meth[nmeths++] = METH_AUTOREGP ;
00229 if( ++nopt >= argc ){
00230 fprintf(stderr,"*** -autoreg needs an argument!\n"); exit(1);
00231 }
00232 meth[nmeths++] = atoi(argv[nopt++]);
00233 if (meth[nmeths - 1] == 0) {
00234 addBriks++;
00235 } else {
00236 nbriks+=meth[nmeths - 1] ;
00237 }
00238 continue ;
00239 }
00240
00241
00242
00243 if( strcmp(argv[nopt],"-prefix") == 0 ){
00244 if( ++nopt >= argc ){
00245 fprintf(stderr,"*** -prefix needs an argument!\n"); exit(1);
00246 }
00247 MCW_strncpy(prefix,argv[nopt],THD_MAX_PREFIX) ;
00248 if( !THD_filename_ok(prefix) ){
00249 fprintf(stderr,"*** %s is not a valid prefix!\n",prefix); exit(1);
00250 }
00251 nopt++ ; continue ;
00252 }
00253
00254
00255
00256 if( strcmp(argv[nopt],"-datum") == 0 ){
00257 if( ++nopt >= argc ){
00258 fprintf(stderr,"*** -datum needs an argument!\n"); exit(1);
00259 }
00260 if( strcmp(argv[nopt],"short") == 0 ){
00261 datum = MRI_short ;
00262 } else if( strcmp(argv[nopt],"float") == 0 ){
00263 datum = MRI_float ;
00264 } else if( strcmp(argv[nopt],"byte") == 0 ){
00265 datum = MRI_byte ;
00266 } else {
00267 fprintf(stderr,"-datum of type '%s' is not supported!\n",
00268 argv[nopt] ) ;
00269 exit(1) ;
00270 }
00271 nopt++ ; continue ;
00272 }
00273
00274
00275
00276 fprintf(stderr,"*** Unknown option: %s\n",argv[nopt]) ; exit(1) ;
00277 }
00278
00279
00280
00281 if (nmeths == 0) nmeths = 1;
00282 if (nbriks == 0 && addBriks == 0) nbriks = 1;
00283
00284
00285
00286 if( nopt >= argc ){
00287 fprintf(stderr,"*** No input dataset!?\n"); exit(1);
00288 }
00289
00290 old_dset = THD_open_dataset( argv[nopt] ) ;
00291 if( !ISVALID_DSET(old_dset) ){
00292 fprintf(stderr,"*** Can't open dataset %s\n",argv[nopt]); exit(1);
00293 }
00294
00295 if( DSET_NVALS(old_dset) < 2 ){
00296 fprintf(stderr,"*** Can't use dataset with < 2 values per voxel!\n") ;
00297 exit(1) ;
00298 }
00299
00300 if( DSET_NUM_TIMES(old_dset) < 2 ){
00301 fprintf(stderr,"--- Input dataset is not 3D+time!\n"
00302 "--- Adding an artificial time axis with dt=1.0\n" ) ;
00303 EDIT_dset_items( old_dset ,
00304 ADN_ntt , DSET_NVALS(old_dset) ,
00305 ADN_ttorg , 0.0 ,
00306 ADN_ttdel , 1.0 ,
00307 ADN_tunits , UNITS_SEC_TYPE ,
00308 NULL ) ;
00309 }
00310
00311
00312
00313
00314 nbriks += ((DSET_NVALS(old_dset)-1) * addBriks);
00315
00316
00317
00318 new_dset = MAKER_4D_to_typed_fbuc(
00319 old_dset ,
00320 prefix ,
00321 datum ,
00322 0 ,
00323 0 ,
00324 nbriks ,
00325 STATS_tsfunc ,
00326 NULL
00327 ) ;
00328
00329 if( new_dset != NULL ){
00330 tross_Copy_History( old_dset , new_dset ) ;
00331 tross_Make_History( "3dTstat" , argc,argv , new_dset ) ;
00332 for (methIndex = 0,brikIndex = 0; methIndex < nmeths; methIndex++, brikIndex++) {
00333 if ((meth[methIndex] == METH_AUTOCORR) || (meth[methIndex] == METH_AUTOREGP)) {
00334 numMultBriks = meth[methIndex+1];
00335 if (numMultBriks == 0) numMultBriks = DSET_NVALS(old_dset);
00336 for (ii = 1; ii <= numMultBriks; ii++) {
00337 char tmpstr[25];
00338 if (meth[methIndex] == METH_AUTOREGP) {
00339 sprintf(tmpstr,"%s[%d](%d)",meth_names[meth[methIndex]],numMultBriks,ii);
00340 } else {
00341 sprintf(tmpstr,"%s(%d)",meth_names[meth[methIndex]],ii);
00342 }
00343 EDIT_BRICK_LABEL(new_dset, (brikIndex + ii - 1), tmpstr) ;
00344 }
00345 methIndex++;
00346 brikIndex += numMultBriks - 1;
00347 } else {
00348 EDIT_BRICK_LABEL(new_dset, brikIndex, meth_names[meth[methIndex]]) ;
00349 }
00350 }
00351 DSET_write( new_dset ) ;
00352 fprintf(stderr,"--- Output dataset %s\n",DSET_BRIKNAME(new_dset)) ;
00353 } else {
00354 fprintf(stderr,"*** Unable to compute output dataset!\n") ;
00355 exit(1) ;
00356 }
00357
00358 exit(0) ;
00359 }
00360
00361
00362
00363
00364
00365 static void STATS_tsfunc( double tzero, double tdelta ,
00366 int npts, float ts[],
00367 double ts_mean, double ts_slope,
00368 void * ud, int nbriks, float * val )
00369 {
00370 static int nvox , ncall ;
00371 int meth_index, ii , out_index;
00372 float* ts_det;
00373
00374
00375
00376 if( val == NULL ){
00377
00378 if( npts > 0 ){
00379
00380 nvox = npts ;
00381 ncall = 0 ;
00382
00383 } else {
00384
00385
00386
00387 }
00388 return ;
00389 }
00390
00391
00392
00393
00394 ts_det = (float*)calloc(npts, sizeof(float));
00395 memcpy( ts_det, ts, npts * sizeof(float));
00396 for( ii = 0; ii < npts; ii++) ts_det[ii] -=
00397 (ts_mean - (ts_slope * (npts - 1) * tdelta/2) + ts_slope * tdelta * ii) ;
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410 for (meth_index = 0, out_index = 0 ; meth_index < nmeths; meth_index++, out_index++) {
00411 switch( meth[meth_index] ){
00412
00413 default:
00414 case METH_MEAN: val[out_index] = ts_mean ; break ;
00415
00416 case METH_SLOPE: val[out_index] = ts_slope ; break ;
00417
00418 case METH_CVAR_NOD:
00419 case METH_SIGMA_NOD:
00420 case METH_CVAR:
00421 case METH_SIGMA:{
00422 register int ii ;
00423 register double sum ;
00424
00425 sum = 0.0 ;
00426 if((meth[meth_index] == METH_CVAR) || (meth[meth_index] == METH_SIGMA )){
00427 for( ii=0 ; ii < npts ; ii++ ) sum += ts_det[ii] * ts_det[ii] ;
00428 } else {
00429 for( ii=0 ; ii < npts ; ii++ ) sum += (ts[ii]-ts_mean)
00430 *(ts[ii]-ts_mean) ;
00431 }
00432
00433 sum = sqrt( sum/(npts-1) ) ;
00434
00435 if((meth[meth_index] == METH_SIGMA) || (meth[meth_index] == METH_SIGMA_NOD)) val[out_index] = sum ;
00436 else if( ts_mean != 0.0 ) val[out_index] = sum / fabs(ts_mean) ;
00437 else val[out_index] = 0.0 ;
00438 }
00439 break ;
00440
00441
00442
00443
00444 case METH_MEDIAN:{
00445 float* ts_copy;
00446 ts_copy = (float*)calloc(npts, sizeof(float));
00447 memcpy( ts_copy, ts, npts * sizeof(float));
00448 val[out_index] = qmed_float( npts , ts_copy ) ;
00449 free(ts_copy);
00450 }
00451 break ;
00452
00453 case METH_MAD:{
00454 float* ts_copy;
00455 register int ii ;
00456 register float vm ;
00457 ts_copy = (float*)calloc(npts, sizeof(float));
00458 memcpy( ts_copy, ts, npts * sizeof(float));
00459 vm = qmed_float( npts , ts_copy ) ;
00460 for( ii=0 ; ii < npts ; ii++ ) ts_copy[ii] = fabs(ts_copy[ii]-vm) ;
00461 val[out_index] = qmed_float( npts , ts_copy ) ;
00462 free(ts_copy);
00463 }
00464 break ;
00465
00466 case METH_ARGMIN:
00467 case METH_MIN:{
00468 register int ii,outdex=0 ;
00469 register float vm=ts[0] ;
00470 for( ii=1 ; ii < npts ; ii++ ) if( ts[ii] < vm ) {
00471 vm = ts[ii] ;
00472 outdex = ii ;
00473 }
00474 if (meth[meth_index] == METH_MIN) {
00475 val[out_index] = vm ;
00476 } else {
00477 val[out_index] = outdex ;
00478 }
00479 }
00480 break ;
00481
00482 case METH_DW:{
00483 register int ii ;
00484 register float den=ts[0]*ts[0] ;
00485 register float num=0 ;
00486 for( ii=1 ; ii < npts ; ii++ ) {
00487 num = num + (ts_det[ii] - ts_det[ii-1])
00488 *(ts_det[ii] - ts_det[ii-1]);
00489 den = den + ts_det[ii] * ts_det[ii];
00490 }
00491 if (den == 0) {
00492 val[out_index] = 0 ;
00493 } else {
00494 val[out_index] = num/den ;
00495 }
00496 }
00497 break ;
00498
00499 case METH_ABSMAX:
00500 case METH_ARGMAX:
00501 case METH_ARGABSMAX:
00502 case METH_MAX:{
00503 register int ii, outdex=0 ;
00504 register float vm=ts[0] ;
00505 if ((meth[meth_index] == METH_ABSMAX) || (meth[meth_index] == METH_ARGABSMAX)) {
00506 vm = fabs(vm) ;
00507 for( ii=1 ; ii < npts ; ii++ ) {
00508 if( fabs(ts[ii]) > vm ) {
00509 vm = fabs(ts[ii]) ;
00510 outdex = ii ;
00511 }
00512 }
00513 } else {
00514 for( ii=1 ; ii < npts ; ii++ ) {
00515 if( ts[ii] > vm ) {
00516 vm = ts[ii] ;
00517 outdex = ii ;
00518 }
00519 }
00520 }
00521
00522 if ((meth[meth_index] == METH_ABSMAX) || (meth[meth_index] == METH_MAX)) {
00523 val[out_index] = vm ;
00524 } else {
00525 val[out_index] = outdex ;
00526 }
00527 }
00528 break ;
00529
00530 case METH_AUTOCORR:{
00531 int numVals;
00532 float* ts_corr;
00533
00534
00535
00536
00537
00538
00539 numVals = meth[meth_index + 1];
00540 if (numVals == 0) numVals = npts - 1;
00541 ts_corr = (float*)calloc(numVals,sizeof(float));
00542
00543 autocorr(npts,ts,numVals,ts_corr);
00544
00545
00546
00547 for( ii = 0; ii < numVals; ii++) {
00548 val[out_index + ii] = ts_corr[ii];
00549 }
00550
00551
00552
00553
00554 meth_index++;
00555
00556
00557
00558
00559 out_index+=(numVals - 1);
00560 free(ts_corr);
00561 }
00562 break ;
00563
00564 case METH_AUTOREGP:{
00565 int numVals,kk,mm;
00566 float *ts_corr, *y, *z;
00567 float alpha, beta, tmp;
00568
00569
00570
00571
00572
00573
00574
00575 numVals = meth[meth_index + 1];
00576 if (numVals == 0) numVals = npts - 1;
00577
00578
00579
00580
00581
00582
00583 ts_corr = (float*)malloc((numVals) * sizeof(float));
00584 y = (float*)malloc(numVals * sizeof(float));
00585 z = (float*)malloc(numVals * sizeof(float));
00586
00587
00588
00589 autocorr(npts,ts,numVals,ts_corr);
00590
00591
00592
00593
00594
00595
00596
00597
00598
00599
00600 y[0] = -ts_corr[0];
00601 alpha = -ts_corr[0];
00602 beta = 1;
00603 for (kk = 0 ; kk < (numVals - 1) ; kk++ ) {
00604 beta = (1 - alpha * alpha ) * beta ;
00605 tmp = 0;
00606 for ( mm = 0 ; mm <= kk ; mm++) {
00607 tmp = tmp + ts_corr[kk - mm] * y[mm];
00608 }
00609 alpha = - ( ts_corr[kk+1] + tmp ) / beta ;
00610 for ( mm = 0 ; mm <= kk ; mm++ ) {
00611 z[mm] = y[mm] + alpha * y[kk - mm];
00612 }
00613 for ( mm = 0 ; mm <= kk ; mm++ ) {
00614 y[mm] = z[mm];
00615 }
00616 y[kk+1] = alpha ;
00617 }
00618
00619
00620
00621
00622 for( ii = 0; ii < numVals; ii++) {
00623 val[out_index + ii] = y[ii];
00624 if (!finite(y[ii])) {
00625 fprintf(stderr,"BAD NUMBER y[%d] = %f, Call# %d\n",ii,y[ii],ncall);
00626 }
00627 }
00628
00629
00630
00631
00632 meth_index++;
00633
00634
00635
00636
00637 out_index+=(numVals - 1);
00638 free(ts_corr);
00639 free(y);
00640 free(z);
00641 }
00642 break ;
00643 }
00644 }
00645
00646 free(ts_det);
00647 ncall++ ; return ;
00648 }
00649
00650
00651 static void autocorr( int npts, float in_ts[], int numVals, float outcoeff[] )
00652 {
00653
00654
00655 int ii,nfft;
00656 double scaler;
00657 complex * cxar = NULL;
00658
00659
00660
00661 nfft = csfft_nextup_one35(npts * 2 - 1);
00662
00663
00664 cxar = (complex *) calloc( sizeof(complex) , nfft ) ;
00665
00666
00667 for( ii=0 ; ii < npts ; ii++ ){
00668 cxar[ii].r = in_ts[ii]; cxar[ii].i = 0.0;
00669 }
00670
00671 for( ii=npts ; ii < nfft ; ii++ ){ cxar[ii].r = cxar[ii].i = 0.0; }
00672
00673
00674
00675
00676
00677
00678 csfft_cox( -1 , nfft , cxar ) ;
00679
00680
00681 for( ii=0 ; ii < nfft ; ii++ ){
00682 cxar[ii].r = CSQR(cxar[ii]) ; cxar[ii].i = 0 ;
00683 }
00684
00685
00686
00687
00688
00689 csfft_scale_inverse( 1 ) ;
00690 csfft_cox( 1 , nfft, cxar ) ;
00691
00692
00693
00694
00695
00696
00697
00698
00699 scaler = cxar[0].r/npts;
00700 for (ii = 0 ; ii < numVals ; ii++ ) {
00701 outcoeff[ii] = cxar[ii+1].r/((npts - (ii+1)) * scaler);
00702 }
00703 free(cxar);
00704 cxar = NULL;
00705 }