00001
00002
00003
00004
00005
00006
00007 #include <stdio.h>
00008 #include <stdlib.h>
00009 #include <string.h>
00010 #include <math.h>
00011
00012 #include "mrilib.h"
00013
00014
00015
00016 #define MAX_NAME 64
00017 #define SUFFIX_DIFF "diff"
00018 #define SUFFIX_TSPM "tspm"
00019 #define SUFFIX_CORR "corr"
00020
00021 static char TF_pname[MAX_NAME] = "tfim." ;
00022 static MRI_IMARR * TF_set1 = NULL ;
00023 static MRI_IMARR * TF_set2 = NULL ;
00024 static float TF_bval = 0.0 ;
00025 static int TF_use_bval = 0 ;
00026 static float TF_pthresh = 0.0 ;
00027 static int TF_eqcorr = 0 ;
00028 static float TF_eqval = 0.0 ;
00029 static int TF_paired = 0 ;
00030
00031
00032
00033 double qginv( double ) ;
00034 double stas4( double , double ) ;
00035 double stinv( double , double ) ;
00036
00037 void TFIM_syntax( char * ) ;
00038 void TFIM_getopts( int , char * argv[] ) ;
00039
00040
00041
00042 int main( int argc , char * argv[] )
00043 {
00044 int kk , ii,nx,ny,npix , num1,num2 , zerout ;
00045 MRI_IMAGE ** stat_ret ;
00046 MRI_IMAGE * avim1 , * sdim1 , * avim2 , * sdim2 ;
00047 float * avar1 , * sdar1 , * avar2 , * sdar2 ;
00048 MRI_IMAGE * difim , * tspim , * corim , * dofim ;
00049 float * difar , * tspar , * corar , * dofar ;
00050 char name[256] ;
00051 float sdmax , dofbar ;
00052
00053
00054
00055 printf(
00056 "MCW TFIM: t-tests on sets of functional images, by RW Cox\n") ;
00057
00058 machdep() ;
00059
00060 if( argc < 2 ) TFIM_syntax("try tfim -help for usage") ;
00061 else if( strcmp(argv[1],"-help") == 0 ) TFIM_syntax(NULL) ;
00062
00063 TFIM_getopts( argc , argv ) ;
00064
00065 #ifdef TFIM_DEBUG
00066 printf("prefix = %s\n",TF_pname) ;
00067 if( TF_use_bval == 1 ){
00068 printf("bval = %f\n",TF_bval) ;
00069 } else {
00070 printf("set1 has %d images\n",TF_set1->num) ;
00071 }
00072 printf("set2 has %d images\n",TF_set2->num) ;
00073 #endif
00074
00075
00076
00077 nx = TF_set2->imarr[0]->nx ;
00078 ny = TF_set2->imarr[0]->ny ; npix = nx * ny ;
00079
00080 if( TF_set1 != NULL ){
00081 num1 = TF_set1->num ;
00082 for( kk=0 ; kk < num1 ; kk++ ){
00083 (void) mri_stat_seq( TF_set1->imarr[kk] ) ;
00084 mri_free( TF_set1->imarr[kk] ) ;
00085 }
00086 stat_ret = mri_stat_seq( NULL ) ;
00087 avim1 = stat_ret[0] ; avar1 = mri_data_pointer( avim1 ) ;
00088 sdim1 = stat_ret[1] ; sdar1 = mri_data_pointer( sdim1 ) ;
00089 }
00090
00091 num2 = TF_set2->num ;
00092 for( kk=0 ; kk < num2 ; kk++ ){
00093 (void) mri_stat_seq( TF_set2->imarr[kk] ) ;
00094 mri_free( TF_set2->imarr[kk] ) ;
00095 }
00096 stat_ret = mri_stat_seq( NULL ) ;
00097 avim2 = stat_ret[0] ; avar2 = mri_data_pointer( avim2 ) ;
00098 sdim2 = stat_ret[1] ; sdar2 = mri_data_pointer( sdim2 ) ;
00099
00100
00101
00102 difim = mri_new( nx,ny , MRI_float ) ; difar = mri_data_pointer(difim) ;
00103 tspim = mri_new( nx,ny , MRI_float ) ; tspar = mri_data_pointer(tspim) ;
00104 dofim = mri_new( nx,ny , MRI_float ) ; dofar = mri_data_pointer(dofim) ;
00105
00106 zerout = 0 ;
00107
00108 if( TF_use_bval == 1 ){
00109 float scl = 1.0 / sqrt((double) num2) ;
00110
00111 for( ii=0 ; ii < npix ; ii++ ){
00112 difar[ii] = avar2[ii] - TF_bval ;
00113 dofar[ii] = num2 - 1 ;
00114 if( sdar2[ii] > 0 ){
00115 tspar[ii] = difar[ii] / (scl * sdar2[ii]) ;
00116 } else {
00117 tspar[ii] = 0.0 ; zerout++ ;
00118 }
00119 }
00120 } else {
00121 float q1 , q2 ;
00122 float n1i =1.0/num1 , n2i =1.0/num2 ,
00123 n11i=1.0/(num1-1.0), n21i=1.0/(num2-1.0) ;
00124
00125 for( ii=0 ; ii < npix ; ii++ ){
00126 difar[ii] = avar2[ii] - avar1[ii] ;
00127 q1 = SQR(sdar1[ii]) * n1i ;
00128 q2 = SQR(sdar2[ii]) * n2i ;
00129 if( q1 > 0 && q2 > 0 ){
00130 tspar[ii] = difar[ii] / sqrt(q1+q2) ;
00131 dofar[ii] = SQR(q1+q2) / ( q1*q1*n11i + q2*q2*n21i ) ;
00132 } else {
00133 tspar[ii] = 0.0 ; dofar[ii] = num1+num2-2 ; zerout++ ;
00134 }
00135 }
00136 }
00137
00138 if( zerout > 0 ){
00139 printf("** set %d pixels to zero due to 0 variance!\n",zerout) ;
00140 }
00141
00142
00143
00144 if( TF_pthresh > 0.0 ){
00145 float thr ;
00146 if( TF_use_bval == 1 ){
00147 dofbar = dofar[0] ;
00148 thr = stinv( TF_pthresh , dofbar ) ;
00149 printf("-- fixed t-threshold = %g\n",thr) ;
00150 } else {
00151 dofbar = 0.0 ;
00152 for( ii=0 ; ii < npix ; ii++ ) dofbar += dofar[ii] ;
00153 dofbar /= npix ;
00154 thr = stinv( TF_pthresh , dofbar ) ;
00155 printf("-- variable t-threshold about %g (dof mean = %g)\n",
00156 thr,dofbar) ;
00157 }
00158
00159 for( ii=0 ; ii < npix ; ii++ ){
00160 if( TF_use_bval == -1 ) thr = stinv( TF_pthresh , dofar[ii] ) ;
00161 if( fabs(tspar[ii]) < thr ) difar[ii] = 0.0 ;
00162 }
00163 }
00164
00165
00166
00167 sprintf( name , "%s%s" , TF_pname , SUFFIX_DIFF ) ;
00168 printf("-- writing difference file %s\n",name) ;
00169 mri_write( name , difim ) ;
00170
00171 sprintf( name , "%s%s" , TF_pname , SUFFIX_TSPM ) ;
00172 printf("-- writing t-statistic file %s\n",name) ;
00173 mri_write( name , tspim ) ;
00174
00175 if( TF_eqcorr ){
00176 float pth , thr , cth , doff ;
00177
00178 corim = mri_new( nx , ny , MRI_float ) ;
00179 corar = mri_data_pointer( corim ) ;
00180 for( ii=0 ; ii < npix ; ii++ ){
00181 doff = (TF_eqval > 0) ? TF_eqval : dofar[ii] ;
00182 corar[ii] = tspar[ii] / sqrt( doff + SQR(tspar[ii]) ) ;
00183 }
00184
00185 sprintf( name , "%s%s" , TF_pname , SUFFIX_CORR ) ;
00186 printf("-- writing 'correlation' file %s\n",name) ;
00187 mri_write( name , corim ) ;
00188
00189 pth = (TF_pthresh > 0) ? TF_pthresh : 0.001 ;
00190 doff = (TF_eqval > 0) ? TF_eqval : dofbar ;
00191 thr = stinv( pth , doff ) ;
00192 cth = thr / sqrt( doff + SQR(thr) ) ;
00193 printf("-- 'correlation' threshold for p=%g is %g\n", pth,cth ) ;
00194 }
00195
00196 exit(0) ;
00197 }
00198
00199
00200
00201 void TFIM_getopts( int argc , char * argv[] )
00202 {
00203 int nopt = 1 , kk , nx,ny ;
00204
00205
00206
00207 while( nopt < argc ){
00208
00209
00210
00211 if( strncmp(argv[nopt],"-paired",5) == 0 ){
00212 TF_paired = 1 ;
00213 nopt++ ; continue ;
00214 }
00215
00216
00217
00218 if( strncmp(argv[nopt],"-prefix",5) == 0 ){
00219 if( ++nopt >= argc ) TFIM_syntax("-prefix needs a name!") ;
00220 strcpy( TF_pname , argv[nopt] ) ;
00221 kk = strlen(TF_pname) ;
00222 if( TF_pname[kk-1] != '.' ){
00223 TF_pname[kk] = '.' ;
00224 TF_pname[kk+1] = '\0' ;
00225 }
00226 nopt++ ; continue ;
00227 }
00228
00229
00230
00231 if( strncmp(argv[nopt],"-pthresh",5) == 0 ){
00232 char * ch ;
00233
00234 if( ++nopt >= argc ) TFIM_syntax("-pthresh needs a value!");
00235
00236 TF_pthresh = strtod( argv[nopt] , &ch ) ;
00237 if( *ch != '\0' || TF_pthresh <= 0.0 || TF_pthresh > 0.99999 )
00238 TFIM_syntax("value after -pthresh is illegal!") ;
00239
00240 nopt++ ; continue ;
00241 }
00242
00243
00244
00245 if( strncmp(argv[nopt],"-eqcorr",5) == 0 ){
00246 char * ch ;
00247
00248 if( ++nopt >= argc ) TFIM_syntax("-eqcorr needs a value!");
00249
00250 TF_eqval = strtod( argv[nopt] , &ch ) ;
00251 if( *ch != '\0' || TF_eqval < 0.0 )
00252 TFIM_syntax("value after -eqcorr is illegal!") ;
00253
00254 TF_eqcorr = 1 ;
00255 nopt++ ; continue ;
00256 }
00257
00258
00259
00260
00261
00262 if( strncmp(argv[nopt],"-base1",5) == 0 ){
00263 char * ch ;
00264
00265 if( ++nopt >= argc ) TFIM_syntax("-base1 needs a value!");
00266 if( TF_use_bval == -1 ) TFIM_syntax("-base1 with -set1 illegal!");
00267
00268 TF_bval = strtod( argv[nopt] , &ch ) ;
00269 if( *ch != '\0' ) TFIM_syntax("value after -base1 is illegal!") ;
00270
00271 TF_use_bval = 1 ;
00272 nopt++ ; continue ;
00273 }
00274
00275
00276
00277 if( strncmp(argv[nopt],"-set1",5) == 0 ){
00278 if( TF_use_bval == 1 ) TFIM_syntax("-set1 with -base1 illegal!");
00279 TF_use_bval = -1 ;
00280
00281 for( kk=nopt+1 ; kk < argc ; kk++ )
00282 if( strncmp(argv[kk],"-set2",5) == 0 ) break ;
00283
00284 if( kk >= argc ) TFIM_syntax("-set1 not followed by -set2") ;
00285 if( kk-1-nopt <= 0 ) TFIM_syntax("-set1 has no files after it") ;
00286
00287 TF_set1 = mri_read_many_files( kk-1-nopt , argv+(nopt+1) ) ;
00288 if( TF_set1 == NULL )
00289 TFIM_syntax("cannot continue without -set1 images") ;
00290
00291 nopt = kk ; continue ;
00292 }
00293
00294
00295
00296 if( strncmp(argv[nopt],"-set2",5) == 0 ){
00297 if( ++nopt >= argc ) TFIM_syntax("-set2 has not files after it") ;
00298
00299 TF_set2 = mri_read_many_files( argc-nopt , argv+nopt ) ;
00300 if( TF_set2 == NULL )
00301 TFIM_syntax("cannot continue without -set2 images") ;
00302
00303 break ;
00304 }
00305
00306
00307
00308 fprintf(stderr,"*** can't interpret this option: %s\n",argv[nopt]) ;
00309 TFIM_syntax("try tfim -help for usage details") ;
00310 }
00311
00312
00313
00314 if( TF_use_bval == -1 &&
00315 ( TF_set1 == NULL || TF_set1->num < 2 ) )
00316 TFIM_syntax("-set1 has too few files in it!") ;
00317
00318 if( TF_set2 == NULL || TF_set2->num < 2 )
00319 TFIM_syntax("-set2 has too few files in it!") ;
00320
00321 if( TF_use_bval == 1 && TF_paired == 1 )
00322 TFIM_syntax("-paired and -base1 are mutually exclusive!") ;
00323
00324 if( TF_paired == 1 && TF_set1 == NULL )
00325 TFIM_syntax("-paired requires presence of -set1!") ;
00326
00327 if( TF_paired == 1 && TF_set1->num != TF_set2->num ){
00328 char str[256] ;
00329 sprintf(str,"-paired requires equal size images sets,\n"
00330 "but -set1 has %d images and -set2 has %d images" ,
00331 TF_set1->num , TF_set2->num ) ;
00332 TFIM_syntax(str) ;
00333 }
00334
00335
00336
00337 nx = TF_set2->imarr[0]->nx ;
00338 ny = TF_set2->imarr[0]->ny ;
00339
00340 for( kk=1 ; kk < TF_set2->num ; kk++ ){
00341 if( nx != TF_set2->imarr[kk]->nx || ny != TF_set2->imarr[kk]->ny ){
00342 fprintf(stderr,
00343 "*** image %d in -set2 not conformant to image 0\n",
00344 kk) ;
00345 TFIM_syntax("cannot continue with images whose sizes differ!") ;
00346 }
00347 }
00348
00349 if( TF_set1 != NULL ){
00350 for( kk=0 ; kk < TF_set1->num ; kk++ ){
00351 if( nx != TF_set1->imarr[kk]->nx || ny != TF_set1->imarr[kk]->ny ){
00352 fprintf(stderr,
00353 "*** image %d in -set1 not conformant to image 0 in -set2\n",
00354 kk) ;
00355 TFIM_syntax("cannot continue with images whose sizes differ!") ;
00356 }
00357 }
00358 }
00359
00360
00361
00362
00363 if( TF_paired == 1 ){
00364 MRI_IMAGE * im1 , * im2 ;
00365 float * ar1 , * ar2 ;
00366 int ii , npix = nx * ny ;
00367
00368 for( kk=0 ; kk < TF_set1->num ; kk++ ){
00369 im1 = mri_to_float(TF_set1->imarr[kk]) ; mri_free(TF_set1->imarr[kk]) ;
00370 im2 = mri_to_float(TF_set2->imarr[kk]) ; mri_free(TF_set2->imarr[kk]) ;
00371 ar1 = mri_data_pointer(im1) ; ar2 = mri_data_pointer(im2) ;
00372 for( ii=0 ; ii < npix ; ii++ ) ar2[ii] -= ar1[ii] ;
00373 mri_free(im1) ; TF_set2->imarr[kk] = im2 ;
00374 }
00375
00376 FREE_IMARR( TF_set1 ) ; TF_set1 = NULL ;
00377 TF_use_bval = 1 ;
00378 TF_bval = 0.0 ;
00379 }
00380
00381 return ;
00382 }
00383
00384
00385 void TFIM_syntax( char * str )
00386 {
00387 if( str != NULL ){
00388 fprintf(stderr,"*** %s\n",str) ;
00389 exit(-1) ;
00390 }
00391
00392 printf(
00393 "\n"
00394 "Usage 1: tfim [options] -set1 image_files ... -set2 image_files ...\n"
00395 "Usage 2: tfim [options] -base1 bval -set2 image_files ...\n"
00396 "\n"
00397 "In usage 1, the collection of images files after '-set1' and the\n"
00398 "collection after '-set2' are averaged and differenced, and the\n"
00399 "difference is tested for significance with a 2 sample Student t-test.\n"
00400 "\n"
00401 "In usage 2, the collection of image files after '-set2' is averaged\n"
00402 "and then has the constant numerical value 'bval' subtracted, and the\n"
00403 "difference is tested for significance with a 1 sample Student t-test.\n"
00404 "\n"
00405 "N.B.: The input images can be in the usual 'short' or 'byte'\n"
00406 " formats, or in the floating point 'flim' format.\n"
00407 "N.B.: If in either set of images, a given pixel has zero variance\n"
00408 " (i.e., is constant), then the t-test is not performed.\n"
00409 " In that pixel, the .tspm file will be zero.\n"
00410 "\n"
00411 "Options are:\n"
00412 "\n"
00413 " -prefix pname: 'pname' is used as the prefix for the output\n"
00414 " filenames. The output image files are\n"
00415 " + pname.diff = average of set2 minus average of set1\n"
00416 " (or minus 'bval')\n"
00417 " + pname.tspm = t-statistic of difference\n"
00418 " Output images are in the 'flim' (floating pt. image)\n"
00419 " format, and may be converted to 16 bit shorts using\n"
00420 " the program 'ftosh'.\n"
00421 " *** The default 'pname' is 'tfim', if -prefix isn't used.\n"
00422 " -pthresh pval: 'pval' is a numeric value between 0 and 1, giving\n"
00423 " the significance level (per voxel) to threshold the\n"
00424 " output with; voxels with (2-sided) t-statistic\n"
00425 " less significant than 'pval' will have their diff\n"
00426 " output zeroed.\n"
00427 " *** The default is no threshold, if -pthresh isn't used.\n"
00428 " -eqcorr dval: If present, this option means to write out the file\n"
00429 " pname.corr = equivalent correlation statistic\n"
00430 " = t/sqrt(dof+t^2)\n"
00431 " The number 'dval' is the value to use for 'dof' if\n"
00432 " dval is positive. This would typically be the total\n"
00433 " number of data images used in forming the image sets,\n"
00434 " if the image sets are from sfim or fim.\n"
00435 " If dval is zero, then dof is computed from the number\n"
00436 " of images in -set1 and -set2; if these are averages\n"
00437 " from program sfim, then dof will be smallish, which in\n"
00438 " turn means that significant corr values will be higher\n"
00439 " than you may be used to from using program fim.\n"
00440 " *** The default is not to write, if -eqcorr isn't used.\n"
00441 " -paired: If present, this means that -set1 and -set2 should be\n"
00442 " compared using a paired sample t-test. This option is\n"
00443 " illegal with the -base1 option. The number of samples\n"
00444 " in the two sets of images must be equal.\n"
00445 " [This test is implemented by subtracting -set1 images\n"
00446 " from the -set2 images, then testing as in '-base1 0'.]\n"
00447 " *** The default is to do an unpaired test, if -paired isn't\n"
00448 " used. In that case, -set1 and -set2 don't need to have\n"
00449 " the same number of images.\n"
00450 ) ;
00451 exit(0) ;
00452 }
00453
00454
00455
00456
00457
00458
00459
00460
00461
00462
00463
00464 double stinv( double p , double nu )
00465 {
00466 double xg , t4 ;
00467 xg = qginv(0.5*p) ;
00468 t4 = stas4( xg , nu ) ;
00469 return t4 ;
00470 }
00471
00472 double stas4( double x , double nu)
00473 {
00474 double t1,t2,t8,t9,t14,t17,t26,t34,t37 ;
00475 t1 = x*x;
00476 t2 = t1*x;
00477 t8 = t1*t1;
00478 t9 = t8*x;
00479 t14 = nu*nu;
00480 t17 = t8*t2;
00481 t26 = t8*t8;
00482 t34 = t14*t14;
00483 t37 = x+(t2/4+x/4)/nu
00484 +(5.0/96.0*t9+t2/6+x/32)/t14
00485 +(t17/128+19.0/384.0*t9
00486 +17.0/384.0*t2-5.0/128.0*x)/t14/nu
00487 +(79.0/92160.0*t26*x+97.0/11520.0*t17+247.0/15360.0*t9
00488 -t2/48-21.0/2048.0*x)/t34;
00489 return t37 ;
00490 }
00491
00492
00493
00494 double qginv( double p )
00495 {
00496 double dp , dx , dt , ddq , dq ;
00497 int newt ;
00498
00499 dp = (p <= 0.5) ? (p) : (1.0-p) ;
00500
00501 if( dp <= 0.0 ){
00502 dx = 13.0 ;
00503 return ( (p <= 0.5) ? (dx) : (-dx) ) ;
00504 }
00505
00506
00507
00508 dt = sqrt( -2.0 * log(dp) ) ;
00509 dx = dt
00510 - ((.010328e+0*dt + .802853e+0)*dt + 2.525517e+0)
00511 /(((.001308e+0*dt + .189269e+0)*dt + 1.432788e+0)*dt + 1.e+0) ;
00512
00513 #if 0
00514
00515
00516 for( newt=0 ; newt < 3 ; newt++ ){
00517 dq = 0.5e+0 * erfc( dx / 1.414213562373095e+0 ) - dp ;
00518 ddq = exp( -0.5e+0 * dx * dx ) / 2.506628274631000e+0 ;
00519 dx = dx + dq / ddq ;
00520 }
00521 #endif
00522
00523 return ( (p <= 0.5) ? (dx) : (-dx) ) ;
00524 }