00001
00002
00003
00004
00005
00006
00007 #include "mrilib.h"
00008
00009 #undef TTDEBUG
00010 #undef USE_PTHRESH
00011
00012
00013
00014
00015
00016 static EDIT_options TT_edopt ;
00017
00018 static int TT_paired = 0 ;
00019 static int TT_pooled = 1 ;
00020 static float TT_pthresh = 0.0 ;
00021 static int TT_use_bval = 0 ;
00022 static float TT_bval = 0.0 ;
00023 static int TT_use_editor = 0 ;
00024 static int TT_be_quiet = 0 ;
00025 static int TT_workmem = 12 ;
00026
00027 #define MEGA 1048576
00028
00029 static THD_string_array * TT_set1 = NULL ;
00030 static THD_string_array * TT_set2 = NULL ;
00031
00032 static char TT_session[THD_MAX_NAME] = "./" ;
00033 static char TT_prefix[THD_MAX_PREFIX] = "tdif" ;
00034 static char TT_label[THD_MAX_LABEL] = "\0" ;
00035
00036 static int TT_datum = ILLEGAL_TYPE ;
00037
00038 static char TT_dof_prefix[THD_MAX_PREFIX] = "\0" ;
00039
00040
00041 void TT_read_opts( int , char ** ) ;
00042 void TT_syntax(char *) ;
00043
00044
00045
00046
00047
00048 #ifdef TTDEBUG
00049 # define DUMP1 fprintf(stderr,"ARG: %s\n",argv[nopt])
00050 # define DUMP2 fprintf(stderr,"ARG: %s %s\n",argv[nopt],argv[nopt+1])
00051 # define DUMP3 fprintf(stderr,"ARG: %s %s %s\n",argv[nopt],argv[nopt+1],argv[nopt+2])
00052 #else
00053 # define DUMP1
00054 # define DUMP2
00055 # define DUMP3
00056 #endif
00057
00058 void TT_read_opts( int argc , char * argv[] )
00059 {
00060 int nopt = 1 ;
00061 float val ;
00062 int ival , kk ;
00063
00064 INIT_EDOPT( &TT_edopt ) ;
00065
00066 while( nopt < argc && argv[nopt][0] == '-' ){
00067
00068
00069
00070 ival = EDIT_check_argv( argc , argv , nopt , &TT_edopt ) ;
00071 if( ival > 0 ){
00072 TT_use_editor = 1 ;
00073 nopt += ival ; continue ;
00074 }
00075
00076
00077
00078 if( strncmp(argv[nopt],"-quiet",6) == 0 ){
00079 DUMP1 ;
00080 TT_be_quiet = 1 ;
00081 nopt++ ; continue ;
00082 }
00083
00084
00085
00086 if( strncmp(argv[nopt],"-workmem",6) == 0 ){
00087 nopt++ ;
00088 if( nopt >= argc ) TT_syntax("need argument after -workmem!") ;
00089 ival = strtol( argv[nopt] , NULL , 10 ) ;
00090 if( ival <= 0 ) TT_syntax("illegal argument after -workmem!") ;
00091 TT_workmem = ival ;
00092 nopt++ ; continue ;
00093 }
00094
00095
00096
00097 if( strncmp(argv[nopt],"-datum",6) == 0 ){
00098 if( ++nopt >= argc ) TT_syntax("need an argument after -datum!") ;
00099
00100 if( strcmp(argv[nopt],"short") == 0 ){
00101 TT_datum = MRI_short ;
00102 } else if( strcmp(argv[nopt],"float") == 0 ){
00103 TT_datum = MRI_float ;
00104 } else {
00105 char buf[256] ;
00106 sprintf(buf,"-datum of type '%s' is not supported in 3dttest!",
00107 argv[nopt] ) ;
00108 TT_syntax(buf) ;
00109 }
00110 nopt++ ; continue ;
00111 }
00112
00113
00114
00115 if( strncmp(argv[nopt],"-session",6) == 0 ){
00116 DUMP2 ;
00117 nopt++ ;
00118 if( nopt >= argc ) TT_syntax("need argument after -session!") ;
00119 MCW_strncpy( TT_session , argv[nopt++] , THD_MAX_NAME ) ;
00120 continue ;
00121 }
00122
00123
00124
00125 if( strncmp(argv[nopt],"-prefix",6) == 0 ){
00126 DUMP2 ;
00127 nopt++ ;
00128 if( nopt >= argc ) TT_syntax("need argument after -prefix!") ;
00129 MCW_strncpy( TT_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00130 continue ;
00131 }
00132
00133 #if 0
00134
00135
00136 if( strncmp(argv[nopt],"-label",6) == 0 ){
00137 DUMP2 ;
00138 nopt++ ;
00139 if( nopt >= argc ) TT_syntax("need argument after -label!") ;
00140 MCW_strncpy( TT_label , argv[nopt++] , THD_MAX_LABEL ) ;
00141 continue ;
00142 }
00143 #endif
00144
00145
00146
00147 if( strncmp(argv[nopt],"-paired",6) == 0 ){
00148 TT_paired = 1 ;
00149 nopt++ ; continue ;
00150 }
00151
00152
00153
00154 if( strncmp(argv[nopt],"-unpooled",6) == 0 ){
00155 TT_pooled = 0 ;
00156 nopt++ ; continue ;
00157 }
00158
00159
00160
00161 if( strncmp(argv[nopt],"-dof_prefix",6) == 0 ){
00162 DUMP2 ;
00163 nopt++ ;
00164 if( nopt >= argc ) TT_syntax("need argument after -dof_prefix!") ;
00165 MCW_strncpy( TT_dof_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00166 continue ;
00167 }
00168
00169 #ifdef USE_PTHRESH
00170
00171
00172 if( strncmp(argv[nopt],"-pthresh",6) == 0 ){
00173 char * ch ;
00174
00175 if( ++nopt >= argc ) TT_syntax("-pthresh needs a value!");
00176 TT_pthresh = strtod( argv[nopt] , &ch ) ;
00177 if( *ch != '\0' || TT_pthresh <= 0.0 || TT_pthresh > 0.99999 )
00178 TT_syntax("value after -pthresh is illegal!") ;
00179
00180 fprintf(stderr,
00181 "*** -pthresh not implemented yet, will ignore!\n" ) ;
00182 TT_pthresh = 0.0 ;
00183
00184 nopt++ ; continue ;
00185 }
00186 #endif
00187
00188
00189
00190
00191
00192 if( strncmp(argv[nopt],"-base1",6) == 0 ){
00193 char * ch ;
00194
00195 if( ++nopt >= argc ) TT_syntax("-base1 needs a value!");
00196 if( TT_use_bval == -1 ) TT_syntax("-base1 with -set1 illegal!");
00197 TT_bval = strtod( argv[nopt] , &ch ) ;
00198 if( *ch != '\0' ) TT_syntax("value after -base1 is illegal!") ;
00199 TT_use_bval = 1 ;
00200 nopt++ ; continue ;
00201 }
00202
00203
00204
00205 if( strncmp(argv[nopt],"-set1",6) == 0 ){
00206 if( TT_use_bval == 1 ) TT_syntax("-set1 with -base1 illegal!");
00207 TT_use_bval = -1 ;
00208 INIT_SARR( TT_set1 ) ;
00209 for( kk=nopt+1 ; kk < argc ; kk++ ){
00210 if( strncmp(argv[kk],"-set2",6) == 0 ) break ;
00211 ADDTO_SARR( TT_set1 , argv[kk] ) ;
00212 }
00213
00214 if( kk >= argc ) TT_syntax("-set1 not followed by -set2") ;
00215 if( kk-1-nopt <= 0 ) TT_syntax("-set1 has no datasets after it") ;
00216 if( TT_set1->num < 2 ) TT_syntax("-set1 doesn't have enough datasets") ;
00217 nopt = kk ; continue ;
00218 }
00219
00220
00221
00222 if( strncmp(argv[nopt],"-set2",6) == 0 ){
00223 INIT_SARR( TT_set2 ) ;
00224 for( kk=nopt+1 ; kk < argc ; kk++ ){
00225 ADDTO_SARR( TT_set2 , argv[kk] ) ;
00226 }
00227 if( TT_set2->num < 2 ) TT_syntax("-set2 doesn't have enough datasets") ;
00228 break ;
00229 }
00230
00231
00232
00233 fprintf(stderr,"*** unrecognized option %s\n",argv[nopt]) ;
00234 exit(1) ;
00235
00236 }
00237
00238 if( strlen(TT_label) == 0 ){
00239 MCW_strncpy(TT_label,TT_prefix,THD_MAX_LABEL) ;
00240 }
00241
00242
00243
00244 if( TT_use_bval == 0 )
00245 TT_syntax("neither -base1 nor -set1 is present!") ;
00246
00247 if( TT_use_bval == -1 &&
00248 ( TT_set1 == NULL || TT_set1->num < 2 ) )
00249 TT_syntax("-set1 has too few datasets in it!") ;
00250
00251 if( TT_set2 == NULL || TT_set2->num < 2 )
00252 TT_syntax("-set2 has too few datasets in it!") ;
00253
00254 if( TT_use_bval == 1 && TT_paired == 1 )
00255 TT_syntax("-paired and -base1 are mutually exclusive!") ;
00256
00257 if( TT_paired == 1 && TT_set1 == NULL )
00258 TT_syntax("-paired requires presence of -set1!") ;
00259
00260 if( TT_paired == 1 && TT_set1->num != TT_set2->num ){
00261 char str[256] ;
00262 sprintf(str,"-paired requires equal size dataset collections,\n"
00263 "but -set1 has %d datasets and -set2 has %d datasets" ,
00264 TT_set1->num , TT_set2->num ) ;
00265 TT_syntax(str) ;
00266 }
00267
00268 if( TT_pooled == 0 && TT_paired == 1 )
00269 TT_syntax("-paired and -unpooled are mutually exclusive!") ;
00270
00271 if( TT_pooled == 0 && TT_use_bval == 1 )
00272 TT_syntax("-base1 and -unpooled are mutually exclusive!") ;
00273
00274 if( TT_pooled == 1 && TT_dof_prefix[0] != '\0' )
00275 fprintf(stderr,"** WARNING: -dof_prefix is used only with -unpooled!\n");
00276
00277 #ifdef TTDEBUG
00278 printf("*** finished with options\n") ;
00279 #endif
00280
00281 return ;
00282 }
00283
00284
00285
00286 void TT_syntax(char * msg)
00287 {
00288 if( msg != NULL ){
00289 fprintf(stderr,"*** %s\n",msg) ;
00290 exit(1) ;
00291 }
00292
00293 printf(
00294 "Gosset (Student) t-test sets of 3D datasets\n"
00295 "Usage 1: 3dttest [options] -set1 datasets ... -set2 datasets ...\n"
00296 " for comparing the means of 2 sets of datasets (voxel by voxel).\n"
00297 "\n"
00298 "Usage 2: 3dttest [options] -base1 bval -set2 datasets ...\n"
00299 " for comparing the mean of 1 set of datasets against a constant.\n"
00300 "\n"
00301
00302 "OUTPUTS:\n"
00303 " A single dataset is created that is the voxel-by-voxel difference\n"
00304 " of the mean of set2 minus the mean of set1 (or minus 'bval').\n"
00305 " The output dataset will be of the intensity+Ttest ('fitt') type.\n"
00306 " The t-statistic at each voxel can be used as an interactive\n"
00307 " thresholding tool in AFNI.\n"
00308 #ifdef USE_PTHRESH
00309 " If the -pthresh option (below) is used to threshold for significance,\n"
00310 " then the t-values will NOT be stored for each voxel; that is, the output\n"
00311 " dataset will be of the intensity-only ('fim') type.\n"
00312 #endif
00313 "\n"
00314
00315 "t-TESTING OPTIONS:\n"
00316 " -set1 datasets ... = Specifies the collection of datasets to put into\n"
00317 " the first set. The mean of set1 will be tested\n"
00318 " with a 2-sample t-test against the mean of set2.\n"
00319 " N.B.: -set1 and -base1 are mutually exclusive!\n"
00320 " -base1 bval = 'bval' is a numerical value that the mean of set2\n"
00321 " will be tested against with a 1-sample t-test.\n"
00322 " -set2 datasets ... = Specifies the collection of datasets to put into\n"
00323 " the second set. There must be at least 2 datasets\n"
00324 " in each of set1 (if used) and set2.\n"
00325 " -paired = Specifies the use of a paired-sample t-test to\n"
00326 " compare set1 and set2. If this option is used,\n"
00327 " set1 and set2 must have the same cardinality.\n"
00328 " N.B.: A paired test is intended for use when the set1 and set2\n"
00329 " dataset function values may be pairwise correlated.\n"
00330 " If they are in fact uncorrelated, this test has less\n"
00331 " statistical 'power' than the unpaired (default) t-test.\n"
00332 " This loss of power is the price that is paid for\n"
00333 " insurance against pairwise correlations.\n"
00334 " -unpooled = Specifies that the variance estimates for set1 and\n"
00335 " set2 be computed separately (not pooled together).\n"
00336 " This only makes sense if -paired is NOT given.\n"
00337 " N.B.: If this option is used, the number of degrees\n"
00338 " of freedom per voxel is a variable, rather\n"
00339 " than a constant.\n"
00340 " -dof_prefix ddd = If '-unpooled' is also used, then a dataset with\n"
00341 " prefix 'ddd' will be created that contains the\n"
00342 " degrees of freedom (DOF) in each voxel.\n"
00343 " You can convert the t-value in the -prefix\n"
00344 " dataset to a z-score using the -dof_prefix dataset\n"
00345 " using commands like so:\n"
00346 " 3dcalc -a 'pname+orig[1]' -b ddd+orig \\\n"
00347 " -datum float -prefix ddd_zz -expr 'fitt_t2z(a,b)'\n"
00348 " 3drefit -substatpar 0 fizt ddd_zz+orig\n"
00349 " At present, AFNI is incapable of directly dealing\n"
00350 " with datasets whose DOF parameter varies between\n"
00351 " voxels. Converting to a z-score (with no parameters)\n"
00352 " is one way of getting around this difficulty.\n"
00353 #ifdef USE_PTHRESH
00354 " -pthresh pval = 'pval' is a probability level (i.e., from 0 to 1)\n"
00355 " at which to threshold the output, per voxel.\n"
00356 " N.B.: NOT IMPLEMENTED YET!\n"
00357 #endif
00358 " -workmem mega = 'mega' specifies the number of megabytes of RAM\n"
00359 " to use for statistical workspace. It defaults\n"
00360 " to 12. The program will run faster if this is\n"
00361 " larger (see the NOTES section below).\n"
00362 "\n"
00363 "The -base1 or -set1 command line switches must follow all other options\n"
00364 "(including those described below) except for the -set2 switch.\n"
00365 "\n"
00366
00367 "INPUT EDITING OPTIONS: The same as are available in 3dmerge.\n"
00368 "\n"
00369
00370 "OUTPUT OPTIONS: these options control the output files.\n"
00371 " -session dirname = Write output into given directory (default=./)\n"
00372 " -prefix pname = Use 'pname' for the output directory prefix\n"
00373 " (default=tdif)\n"
00374 #if 0
00375 " -label string = Use 'string' for the label in the output\n"
00376 " dataset (the label is used for switching\n"
00377 " between datasets in AFNI)\n"
00378 #endif
00379 " -datum type = Use 'type' to store the output difference\n"
00380 " in the means; 'type' may be short or float.\n"
00381 " How the default is determined is described\n"
00382 " in the notes below.\n"
00383 "\n"
00384 "NOTES:\n"
00385 " ** To economize on memory, 3dttest makes multiple passes through\n"
00386 " the input datasets. On each pass, the entire editing process\n"
00387 " will be carried out again. For efficiency's sake, it is\n"
00388 " better to carry out the editing using 3dmerge to produce\n"
00389 " temporary datasets, and then run 3dttest on them. This applies\n"
00390 " with particular force if a 'blurring' option is used.\n"
00391 " Note also that editing a dataset requires that it be read into\n"
00392 " memory in its entirety (so that the disk file is not altered).\n"
00393 " This will increase the memory needs of the program far beyond\n"
00394 " the level set by the -workmem option.\n"
00395 " ** The input datasets are specified by their .HEAD files,\n"
00396 " but their .BRIK files must exist also! This program cannot\n"
00397 " 'warp-on-demand' from other datasets.\n"
00398 " ** This program cannot deal with time-dependent or complex-valued datasets!\n"
00399 " By default, the output dataset function values will be shorts if the\n"
00400 " first input dataset is byte- or short-valued; otherwise they will be\n"
00401 " floats. This behavior may be overridden using the -datum option.\n"
00402 ) ;
00403 printf("\n" MASTER_SHORTHELP_STRING ) ;
00404 exit(0) ;
00405 }
00406
00407
00408
00409 static float ptable[] = { 0.5 , 0.2 , 0.1 , 0.05 , 0.01 , 0.001 , 0.0001 , 0.00001 } ;
00410
00411 int main( int argc , char * argv[] )
00412 {
00413 int nx,ny,nz , nxyz , ii,kk , num1,num2 , num_tt=0 , iv ,
00414 output_type , output_nvals , piece , num_piece , piece_len , fim_offset ;
00415 float dx,dy,dz , dxyz ,
00416 num1_inv , num2_inv , num1m1_inv , num2m1_inv , dof , tthr ,
00417 dd,tt,q1,q2 , f1,f2 , tt_max=0.0 ;
00418 THD_3dim_dataset * dset=NULL , * new_dset=NULL , * qset=NULL ;
00419 float * av1 , * av2 , * sd1 , * sd2 , * ffim , * gfim ;
00420
00421 void * vsp ;
00422 short * tsp , * tsar ;
00423 float * fsp , * fsar ;
00424
00425 void * vdif ;
00426 short * sdif , * sdar ;
00427 float * fdif , * fdar ;
00428
00429 void * vfim ;
00430 char cbuf[THD_MAX_NAME] ;
00431 float fbuf[MAX_STAT_AUX] , fimfac , fimfacinv ;
00432 int output_datum ;
00433 int piece_size ;
00434 float npiece , memuse ;
00435
00436 float *dofbrik=NULL , *dofar=NULL ;
00437 THD_3dim_dataset *dof_dset=NULL ;
00438
00439
00440
00441 if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) TT_syntax(NULL) ;
00442
00443
00444
00445 mainENTRY("3dttest main"); machdep() ; PRINT_VERSION("3dttest") ;
00446
00447 { int new_argc ; char ** new_argv ;
00448 addto_args( argc , argv , &new_argc , &new_argv ) ;
00449 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00450 }
00451
00452 AFNI_logger("3dttest",argc,argv) ;
00453
00454 TT_read_opts( argc , argv ) ;
00455
00456 if( ! TT_be_quiet )
00457 printf("3dttest: t-tests of 3D datasets, by RW Cox\n") ;
00458
00459
00460
00461 dset = THD_open_dataset( TT_set2->ar[0] ) ;
00462 if( ! ISVALID_3DIM_DATASET(dset) ){
00463 fprintf(stderr,"*** Unable to open dataset file %s\n",TT_set2->ar[0]);
00464 exit(1) ;
00465 }
00466
00467 nx = dset->daxes->nxx ;
00468 ny = dset->daxes->nyy ;
00469 nz = dset->daxes->nzz ; nxyz = nx * ny * nz ;
00470 dx = fabs(dset->daxes->xxdel) ;
00471 dy = fabs(dset->daxes->yydel) ;
00472 dz = fabs(dset->daxes->zzdel) ; dxyz = dx * dy * dz ;
00473
00474 #ifdef TTDEBUG
00475 printf("*** nx=%d ny=%d nz=%d\n",nx,ny,nz) ;
00476 #endif
00477
00478
00479
00480 #ifdef TTDEBUG
00481 printf("*** making empty dataset\n") ;
00482 #endif
00483
00484 new_dset = EDIT_empty_copy( dset ) ;
00485
00486 tross_Make_History( "3dttest" , argc,argv , new_dset ) ;
00487
00488 strcpy( cbuf , dset->self_name ) ; strcat( cbuf , "+TT" ) ;
00489
00490 iv = DSET_PRINCIPAL_VALUE(dset) ;
00491
00492 if( TT_datum >= 0 ){
00493 output_datum = TT_datum ;
00494 } else {
00495 output_datum = DSET_BRICK_TYPE(dset,iv) ;
00496 if( output_datum == MRI_byte ) output_datum = MRI_short ;
00497 }
00498
00499 #ifdef TTDEBUG
00500 printf(" ** datum = %s\n",MRI_TYPE_name[output_datum]) ;
00501 #endif
00502
00503 iv = EDIT_dset_items( new_dset ,
00504 ADN_prefix , TT_prefix ,
00505 ADN_label1 , TT_prefix ,
00506 ADN_directory_name , TT_session ,
00507 ADN_self_name , cbuf ,
00508 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE ,
00509 ADN_func_type , FUNC_TT_TYPE ,
00510 ADN_nvals , FUNC_nvals[FUNC_TT_TYPE] ,
00511 ADN_datum_all , output_datum ,
00512 ADN_none ) ;
00513
00514 if( iv > 0 ){
00515 fprintf(stderr,
00516 "*** %d errors in attempting to create output dataset!\n",iv ) ;
00517 exit(1) ;
00518 }
00519
00520 if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
00521 fprintf(stderr,
00522 "*** Output dataset file %s already exists--cannot continue!\a\n",
00523 new_dset->dblk->diskptr->header_name ) ;
00524 exit(1) ;
00525 }
00526
00527 #ifdef TTDEBUG
00528 printf("*** deleting exemplar dataset\n") ;
00529 #endif
00530
00531 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00532
00533
00534
00535 #define MTEST(ptr) \
00536 if((ptr)==NULL) \
00537 ( fprintf(stderr,"*** Cannot allocate memory for statistics!\n" \
00538 "*** Try using the -workmem option to reduce memory needs,\n" \
00539 "*** or create more swap space in the operating system.\n" \
00540 ), exit(0) )
00541
00542
00543
00544 npiece = 3.0 ;
00545 if( TT_paired ) npiece += 1.0 ;
00546 else if( TT_set1 != NULL ) npiece += 2.0 ;
00547
00548 npiece += mri_datum_size(output_datum) / (float) sizeof(float) ;
00549 npiece += mri_datum_size(output_datum) / (float) sizeof(float) ;
00550
00551 piece_size = TT_workmem * MEGA / ( npiece * sizeof(float) ) ;
00552 if( piece_size > nxyz ) piece_size = nxyz ;
00553
00554 #ifdef TTDEBUG
00555 printf("*** malloc-ing space for statistics: %g float arrays of length %d\n",
00556 npiece,piece_size) ;
00557 #endif
00558
00559 av2 = (float *) malloc( sizeof(float) * piece_size ) ; MTEST(av2) ;
00560 sd2 = (float *) malloc( sizeof(float) * piece_size ) ; MTEST(sd2) ;
00561 ffim = (float *) malloc( sizeof(float) * piece_size ) ; MTEST(ffim) ;
00562 num2 = TT_set2->num ;
00563
00564 if( TT_paired ){
00565 av1 = sd1 = NULL ;
00566 gfim = (float *) malloc( sizeof(float) * piece_size ) ; MTEST(gfim) ;
00567 num1 = num2 ;
00568 } else if( TT_set1 != NULL ){
00569 av1 = (float *) malloc( sizeof(float) * piece_size ) ; MTEST(av1) ;
00570 sd1 = (float *) malloc( sizeof(float) * piece_size ) ; MTEST(sd1) ;
00571 gfim = NULL ;
00572 num1 = TT_set1->num ;
00573 } else {
00574 av1 = sd1 = NULL ;
00575 gfim = NULL ;
00576 num1 = 0 ;
00577 }
00578
00579 vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ; MTEST(vdif) ;
00580 vsp = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ; MTEST(vsp) ;
00581
00582
00583
00584 if( TT_pooled == 0 && TT_dof_prefix[0] != '\0' ){
00585 dofbrik = (float *) malloc( sizeof(float) * nxyz ) ; MTEST(dofbrik) ;
00586
00587 dof_dset = EDIT_empty_copy( new_dset ) ;
00588
00589 tross_Make_History( "3dttest" , argc,argv , dof_dset ) ;
00590
00591 EDIT_dset_items( dof_dset ,
00592 ADN_prefix , TT_dof_prefix ,
00593 ADN_directory_name , TT_session ,
00594 ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE,
00595 ADN_func_type , FUNC_BUCK_TYPE ,
00596 ADN_nvals , 1 ,
00597 ADN_datum_all , MRI_float ,
00598 ADN_none ) ;
00599
00600 if( THD_is_file(dof_dset->dblk->diskptr->header_name) ){
00601 fprintf(stderr,
00602 "*** -dof_prefix dataset file %s already exists--cannot continue!\a\n",
00603 dof_dset->dblk->diskptr->header_name ) ;
00604 exit(1) ;
00605 }
00606
00607 EDIT_substitute_brick( dof_dset , 0 , MRI_float , dofbrik ) ;
00608 }
00609
00610
00611
00612 if( ! TT_be_quiet ){
00613 memuse = sizeof(float) * piece_size * npiece
00614 + ( mri_datum_size(output_datum) + sizeof(short) ) * nxyz ;
00615
00616 if( dofbrik != NULL ) memuse += sizeof(float) * nxyz ;
00617
00618 printf("--- allocated %d Megabytes memory for internal use\n",(int)(memuse/MEGA)) ;
00619 }
00620
00621 mri_fix_data_pointer( vdif , DSET_BRICK(new_dset,0) ) ;
00622 mri_fix_data_pointer( vsp , DSET_BRICK(new_dset,1) ) ;
00623
00624 switch( output_datum ){
00625 default: fprintf(stderr,"Illegal input data type %d = %s!\n",
00626 output_datum , MRI_TYPE_name[output_datum] ) ;
00627 exit(1) ;
00628
00629 case MRI_short: sdif = (short *) vdif ; tsp = (short *) vsp ; break ;
00630 case MRI_float: fdif = (float *) vdif ; fsp = (float *) vsp ; break ;
00631 }
00632
00633 num2_inv = 1.0 / num2 ; num2m1_inv = 1.0 / (num2-1) ;
00634 if( num1 > 0 ){
00635 num1_inv = 1.0 / num1 ; num1m1_inv = 1.0 / (num1-1) ;
00636 }
00637
00638
00639
00640
00641
00642 #define DOPEN(ds,name) \
00643 do{ int pv ; (ds) = THD_open_dataset((name)) ; \
00644 if( !ISVALID_3DIM_DATASET((ds)) ){ \
00645 fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
00646 if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || (ds)->daxes->nzz!=nz ){ \
00647 fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; } \
00648 if( DSET_NUM_TIMES((ds)) > 1 ){ \
00649 fprintf(stderr,"*** Can't use time-dependent data: %s\n",(name));exit(1); } \
00650 if( TT_use_editor ) EDIT_one_dataset( (ds), &TT_edopt ) ; \
00651 else THD_load_datablock( (ds)->dblk ) ; \
00652 pv = DSET_PRINCIPAL_VALUE((ds)) ; \
00653 if( DSET_ARRAY((ds),pv) == NULL ){ \
00654 fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \
00655 if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \
00656 fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1); } \
00657 break ; } while (0)
00658
00659
00660
00661 #define SUB_POINTER(ds,vv,ind,ptr) \
00662 do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \
00663 default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \
00664 case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \
00665 (ptr) = (void *)( fim + (ind) ) ; \
00666 } break ; \
00667 case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \
00668 (ptr) = (void *)( fim + (ind) ) ; \
00669 } break ; \
00670 case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \
00671 (ptr) = (void *)( fim + (ind) ) ; \
00672 } break ; } break ; } while(0)
00673
00674
00675
00676 num_piece = (nxyz + piece_size - 1) / piece_size ;
00677
00678 nice(2) ;
00679
00680 for( piece=0 ; piece < num_piece ; piece++ ){
00681
00682 fim_offset = piece * piece_size ;
00683 piece_len = (piece < num_piece-1 ) ? piece_size
00684 : (nxyz - fim_offset) ;
00685
00686 #ifdef TTDEBUG
00687 printf("*** start of piece %d: length=%d offset=%d\n",piece,piece_len,fim_offset) ;
00688 #else
00689 if( ! TT_be_quiet ){
00690 printf("--- starting piece %d/%d (%d voxels) ",piece+1,num_piece,piece_len) ;
00691 fflush(stdout) ;
00692 }
00693 #endif
00694
00695
00696
00697 for( ii=0 ; ii < piece_len ; ii++ ) av2[ii] = 0.0 ;
00698 for( ii=0 ; ii < piece_len ; ii++ ) sd2[ii] = 0.0 ;
00699
00700 for( kk=0 ; kk < num2 ; kk++ ){
00701
00702
00703
00704 DOPEN(dset,TT_set2->ar[kk]) ;
00705 iv = DSET_PRINCIPAL_VALUE(dset) ;
00706
00707 #ifndef TTDEBUG
00708 if( ! TT_be_quiet ){ printf(".") ; fflush(stdout) ; }
00709 #else
00710 printf(" ** opened dataset file %s\n",TT_set2->ar[kk]);
00711 #endif
00712
00713 if( piece == 0 && kk == 0 ){
00714 fimfac = DSET_BRICK_FACTOR(dset,iv) ;
00715 if( fimfac == 0.0 ) fimfac = 1.0 ;
00716 fimfacinv = 1.0 / fimfac ;
00717 #ifdef TTDEBUG
00718 printf(" ** set fimfac = %g\n",fimfac) ;
00719 #endif
00720 }
00721
00722
00723
00724 SUB_POINTER(dset,iv,fim_offset,vfim) ;
00725 EDIT_coerce_scale_type( piece_len , DSET_BRICK_FACTOR(dset,iv) ,
00726 DSET_BRICK_TYPE(dset,iv),vfim ,
00727 MRI_float ,ffim ) ;
00728 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00729
00730
00731
00732 if( TT_paired ){
00733 DOPEN(dset,TT_set1->ar[kk]) ;
00734 iv = DSET_PRINCIPAL_VALUE(dset) ;
00735
00736 #ifndef TTDEBUG
00737 if( ! TT_be_quiet ){ printf(".") ; fflush(stdout) ; }
00738 #else
00739 printf(" ** opened dataset file %s\n",TT_set1->ar[kk]);
00740 #endif
00741
00742 SUB_POINTER(dset,iv,fim_offset,vfim) ;
00743 EDIT_coerce_scale_type( piece_len , DSET_BRICK_FACTOR(dset,iv) ,
00744 DSET_BRICK_TYPE(dset,iv),vfim ,
00745 MRI_float ,gfim ) ;
00746 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00747
00748 for( ii=0 ; ii < piece_len ; ii++ ) ffim[ii] -= gfim[ii] ;
00749 }
00750
00751 #ifdef TTDEBUG
00752 printf(" * adding into av2 and sd2\n") ;
00753 #endif
00754
00755
00756
00757 for( ii=0 ; ii < piece_len ; ii++ ){
00758 dd = ffim[ii] ; av2[ii] += dd ; sd2[ii] += dd * dd ;
00759 }
00760
00761 }
00762
00763
00764
00765 #ifdef TTDEBUG
00766 printf(" ** forming mean and sigma of set2\n") ;
00767 #endif
00768
00769 for( ii=0 ; ii < piece_len ; ii++ ){
00770 av2[ii] *= num2_inv ;
00771 dd = (sd2[ii] - num2*av2[ii]*av2[ii]) ;
00772 sd2[ii] = (dd > 0.0) ? sqrt( num2m1_inv * dd ) : 0.0 ;
00773 }
00774
00775
00776
00777 if( ! TT_paired && TT_set1 != NULL ){
00778
00779 for( ii=0 ; ii < piece_len ; ii++ ) av1[ii] = 0.0 ;
00780 for( ii=0 ; ii < piece_len ; ii++ ) sd1[ii] = 0.0 ;
00781
00782 for( kk=0 ; kk < num1 ; kk++ ){
00783 DOPEN(dset,TT_set1->ar[kk]) ;
00784 iv = DSET_PRINCIPAL_VALUE(dset) ;
00785
00786 #ifndef TTDEBUG
00787 if( ! TT_be_quiet ){ printf(".") ; fflush(stdout) ; }
00788 #else
00789 printf(" ** opened dataset file %s\n",TT_set1->ar[kk]);
00790 #endif
00791
00792 SUB_POINTER(dset,iv,fim_offset,vfim) ;
00793 EDIT_coerce_scale_type( piece_len , DSET_BRICK_FACTOR(dset,iv) ,
00794 DSET_BRICK_TYPE(dset,iv),vfim ,
00795 MRI_float ,ffim ) ;
00796 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00797
00798 #ifdef TTDEBUG
00799 printf(" * adding into av1 and sd1\n") ;
00800 #endif
00801
00802 for( ii=0 ; ii < piece_len ; ii++ ){
00803 dd = ffim[ii] ; av1[ii] += dd ; sd1[ii] += dd * dd ;
00804 }
00805 }
00806
00807
00808
00809 #ifdef TTDEBUG
00810 printf(" ** forming mean and sigma of set1\n") ;
00811 #endif
00812
00813 for( ii=0 ; ii < piece_len ; ii++ ){
00814 av1[ii] *= num1_inv ;
00815 dd = (sd1[ii] - num1*av1[ii]*av1[ii]) ;
00816 sd1[ii] = (dd > 0.0) ? sqrt( num1m1_inv * dd ) : 0.0 ;
00817 }
00818 }
00819
00820
00821
00822 #ifndef TTDEBUG
00823 if( ! TT_be_quiet ){ printf("+") ; fflush(stdout) ; }
00824 #else
00825 printf(" ** computing t-tests next\n") ;
00826 #endif
00827
00828 switch( output_datum ){
00829 case MRI_short:
00830 sdar = sdif + fim_offset ;
00831 tsar = tsp + fim_offset ;
00832 break ;
00833 case MRI_float:
00834 fdar = fdif + fim_offset ;
00835 fsar = fsp + fim_offset ;
00836 break ;
00837 }
00838
00839
00840
00841 #define DIFASS switch( output_datum ){ \
00842 case MRI_short: sdar[ii] = (short) (fimfacinv*dd) ; break ; \
00843 case MRI_float: fdar[ii] = (float) dd ; break ; }
00844
00845 #define TOP_SS 32700
00846 #define TOP_TT (32700.0/FUNC_TT_SCALE_SHORT)
00847
00848 if( TT_paired || TT_use_bval == 1 ){
00849
00850 f2 = 1.0 / sqrt( (double) num2 ) ;
00851 for( ii=0 ; ii < piece_len ; ii++ ){
00852 dd = av2[ii] - TT_bval ; DIFASS ;
00853 if( sd2[ii] > 0.0 ){
00854 num_tt++ ;
00855 tt = dd / (f2 * sd2[ii]) ;
00856 switch( output_datum ){
00857 case MRI_short:
00858 tsar[ii] = (tt>TOP_TT) ? (TOP_SS)
00859 : (tt<-TOP_TT) ? (-TOP_SS)
00860 : ((short)(FUNC_TT_SCALE_SHORT*tt)) ;
00861 break ;
00862 case MRI_float:
00863 fsar[ii] = tt ;
00864 break ;
00865 }
00866 tt = fabs(tt) ; if( tt > tt_max ) tt_max = tt ;
00867 } else {
00868 switch( output_datum ){
00869 case MRI_short: tsar[ii] = 0 ; break ;
00870 case MRI_float: fsar[ii] = 0 ; break ;
00871 }
00872 }
00873 }
00874
00875 #ifdef TTDEBUG
00876 printf(" ** paired or bval test: num_tt = %d\n",num_tt) ;
00877 #endif
00878
00879 } else if( TT_pooled ){
00880
00881 f1 = (num1-1.0) * (1.0/num1 + 1.0/num2) / (num1+num2-2.0) ;
00882 f2 = (num2-1.0) * (1.0/num1 + 1.0/num2) / (num1+num2-2.0) ;
00883 for( ii=0 ; ii < piece_len ; ii++ ){
00884 dd = av2[ii] - av1[ii] ; DIFASS ;
00885 q1 = f1 * sd1[ii]*sd1[ii] + f2 * sd2[ii]*sd2[ii] ;
00886 if( q1 > 0.0 ){
00887 num_tt++ ;
00888 tt = dd / sqrt(q1) ;
00889 switch( output_datum ){
00890 case MRI_short:
00891 tsar[ii] = (tt>TOP_TT) ? (TOP_SS)
00892 : (tt<-TOP_TT) ? (-TOP_SS)
00893 : ((short)(FUNC_TT_SCALE_SHORT*tt)) ;
00894 break ;
00895 case MRI_float:
00896 fsar[ii] = tt ;
00897 break ;
00898 }
00899 tt = fabs(tt) ; if( tt > tt_max ) tt_max = tt ;
00900 } else {
00901 switch( output_datum ){
00902 case MRI_short: tsar[ii] = 0 ; break ;
00903 case MRI_float: fsar[ii] = 0 ; break ;
00904 }
00905 }
00906 }
00907
00908 #ifdef TTDEBUG
00909 printf(" ** pooled test: num_tt = %d\n",num_tt) ;
00910 #endif
00911
00912 } else {
00913
00914
00915 if( dofbrik != NULL ) dofar = dofbrik + fim_offset ;
00916
00917 for( ii=0 ; ii < piece_len ; ii++ ){
00918 dd = av2[ii] - av1[ii] ; DIFASS ;
00919 q1 = num1_inv * sd1[ii]*sd1[ii] ;
00920 q2 = num2_inv * sd2[ii]*sd2[ii] ;
00921 if( q1>0.0 && q2>0.0 ){
00922 num_tt++ ;
00923 tt = dd / sqrt(q1+q2) ;
00924 switch( output_datum ){
00925 case MRI_short:
00926 tsar[ii] = (tt>TOP_TT) ? (TOP_SS)
00927 : (tt<-TOP_TT) ? (-TOP_SS)
00928 : ((short)(FUNC_TT_SCALE_SHORT*tt)) ;
00929 break ;
00930 case MRI_float:
00931 fsar[ii] = tt ;
00932 break ;
00933 }
00934 tt = fabs(tt) ; if( tt > tt_max ) tt_max = tt ;
00935
00936 if( dofar != NULL )
00937 dofar[ii] = (q1+q2)*(q1+q2)
00938 / (num1m1_inv*q1*q1 + num2m1_inv*q2*q2) ;
00939 } else {
00940 switch( output_datum ){
00941 case MRI_short: tsar[ii] = 0 ; break ;
00942 case MRI_float: fsar[ii] = 0 ; break ;
00943 }
00944 if( dofar != NULL ) dofar[ii] = 1.0 ;
00945 }
00946 }
00947
00948 #ifdef TTDEBUG
00949 printf(" ** unpooled test: num_tt = %d\n",num_tt) ;
00950 #endif
00951 }
00952
00953 #ifndef TTDEBUG
00954 if( ! TT_be_quiet ){ printf("\n") ; fflush(stdout) ; }
00955 #endif
00956
00957 }
00958
00959 if( TT_paired ){
00960 printf("--- Number of degrees of freedom = %d (paired test)\n",num2-1) ;
00961 dof = num2 - 1 ;
00962 } else if( TT_use_bval == 1 ){
00963 printf("--- Number of degrees of freedom = %d (1-sample test)\n",num2-1) ;
00964 dof = num2 - 1 ;
00965 } else {
00966 printf("--- Number of degrees of freedom = %d (2-sample test)\n",num1+num2-2) ;
00967 dof = num1+num2-2 ;
00968 if( ! TT_pooled )
00969 printf(" (For unpooled variance estimate, this is only approximate!)\n") ;
00970 }
00971
00972 printf("--- Number of t-tests performed = %d out of %d voxels\n",num_tt,nxyz) ;
00973 printf("--- Largest |t| value found = %g\n",tt_max) ;
00974
00975 kk = sizeof(ptable) / sizeof(float) ;
00976 for( ii=0 ; ii < kk ; ii++ ){
00977 tt = student_p2t( ptable[ii] , dof ) ;
00978 printf("--- Double sided tail p = %8f at t = %8f\n" , ptable[ii] , tt ) ;
00979 }
00980
00981 printf("--- Writing combined dataset into %s\n", DSET_BRIKNAME(new_dset) ) ;
00982
00983 fbuf[0] = dof ;
00984 for( ii=1 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
00985 (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
00986
00987 fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
00988 fbuf[1] = (output_datum == MRI_short ) ? 1.0 / FUNC_TT_SCALE_SHORT : 0.0 ;
00989 (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
00990
00991 THD_load_statistics( new_dset ) ;
00992 THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
00993
00994 if( dof_dset != NULL ){
00995 DSET_write( dof_dset ) ;
00996 fprintf(stderr,"--- Wrote unpooled DOF dataset into %s\n", DSET_BRIKNAME(dof_dset) ) ;
00997 }
00998
00999 exit(0) ;
01000 }
