Doxygen Source Code Documentation
tfim.c File Reference
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "mrilib.h"
Go to the source code of this file.
Defines | |
#define | MAX_NAME 64 |
#define | SUFFIX_DIFF "diff" |
#define | SUFFIX_TSPM "tspm" |
#define | SUFFIX_CORR "corr" |
Functions | |
double | qginv (double) |
double | stas4 (double, double) |
double | stinv (double, double) |
void | TFIM_syntax (char *) |
void | TFIM_getopts (int, char *argv[]) |
int | main (int argc, char *argv[]) |
Variables | |
char | TF_pname [MAX_NAME] = "tfim." |
MRI_IMARR * | TF_set1 = NULL |
MRI_IMARR * | TF_set2 = NULL |
float | TF_bval = 0.0 |
int | TF_use_bval = 0 |
float | TF_pthresh = 0.0 |
int | TF_eqcorr = 0 |
float | TF_eqval = 0.0 |
int | TF_paired = 0 |
Define Documentation
|
|
|
Definition at line 19 of file tfim.c. Referenced by main(). |
|
Definition at line 17 of file tfim.c. Referenced by main(). |
|
Definition at line 18 of file tfim.c. Referenced by main(). |
Function Documentation
|
find asymmetry measures in 1/2 spaces perp to L-R * Definition at line 42 of file tfim.c. References argc, MRI_IMARR::imarr, machdep(), mri_data_pointer(), mri_free(), mri_new(), mri_stat_seq(), mri_write(), name, MRI_IMARR::num, MRI_IMAGE::nx, MRI_IMAGE::ny, SQR, stinv(), SUFFIX_CORR, SUFFIX_DIFF, SUFFIX_TSPM, TF_bval, TF_eqval, TF_pname, TF_pthresh, TF_use_bval, TFIM_getopts(), TFIM_syntax(), and thr.
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 /*----- read inputs -----*/ 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 /*----- form averages of images -----*/ 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 /*----- process set averages into statistics -----*/ 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 ; /* count of pixels zero-ed out due to no variance */ 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 ; /* qi = sigi^2 / numi */ 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 /** threshold, if desired **/ 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 /*----- write output images -----*/ 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 } |
|
solve for x by a modified newton-raphson method * Definition at line 494 of file tfim.c. Referenced by beta_t2z(), binomial_t2z(), chisq_t2z(), correl_t2z(), fstat_t2z(), gamma_t2z(), main(), poisson_t2z(), set_unusuality_tail(), stinv(), studave_p2t(), student_t2z(), THD_outlier_count(), and UC_unusuality().
00495 { 00496 double dp , dx , dt , ddq , dq ; 00497 int newt ; 00498 00499 dp = (p <= 0.5) ? (p) : (1.0-p) ; /* make between 0 and 0.5 */ 00500 00501 if( dp <= 0.0 ){ 00502 dx = 13.0 ; 00503 return ( (p <= 0.5) ? (dx) : (-dx) ) ; 00504 } 00505 00506 /** Step 1: use 26.2.23 from Abramowitz and Stegun **/ 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 /** Step 2: do 3 Newton steps to improve this **/ 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) ) ; /* return with correct sign */ 00524 } |
|
Definition at line 472 of file tfim.c. Referenced by stinv().
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 } |
|
Definition at line 464 of file tfim.c. References p, qginv(), and stas4(). Referenced by main().
|
|
threshold, if desired * Definition at line 201 of file tfim.c. References argc, FREE_IMARR, MRI_IMARR::imarr, mri_data_pointer(), mri_free(), mri_read_many_files(), mri_to_float(), MRI_IMARR::num, MRI_IMAGE::nx, MRI_IMAGE::ny, strtod(), TF_bval, TF_eqcorr, TF_eqval, TF_paired, TF_pname, TF_pthresh, TF_use_bval, and TFIM_syntax(). Referenced by main().
00202 { 00203 int nopt = 1 , kk , nx,ny ; 00204 00205 /*--- scan argument list ---*/ 00206 00207 while( nopt < argc ){ 00208 00209 /** -paired **/ 00210 00211 if( strncmp(argv[nopt],"-paired",5) == 0 ){ 00212 TF_paired = 1 ; 00213 nopt++ ; continue ; /* skip to next arg */ 00214 } 00215 00216 /** -prefix pname **/ 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 ; /* skip to next arg */ 00227 } 00228 00229 /** -pthresh pval **/ 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 ; /* skip to next arg */ 00241 } 00242 00243 /** -eqcorr dval **/ 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 ; /* skip to next arg */ 00256 } 00257 00258 /** after this point, the options are no longer 'free floating' **/ 00259 00260 /** -base1 bval **/ 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 ; /* skip to next arg */ 00273 } 00274 00275 /** -set1 file file ... **/ 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 ; /* skip to arg that matched -set2 */ 00292 } 00293 00294 /** -set2 file file ... */ 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 ; /* end of possible inputs */ 00304 } 00305 00306 /** get to here is bad news! **/ 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 /*--- check arguments for OK-ositiness ---*/ 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 /* check images for consistency */ 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 /* if paired t-test, subtract set1 from set2 00361 to convert it into the equivalent base level test vs. 0 */ 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 } |
|
get to here is bad news! * Definition at line 385 of file tfim.c. Referenced by main(), and TFIM_getopts().
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 } |
Variable Documentation
|
Definition at line 24 of file tfim.c. Referenced by main(), and TFIM_getopts(). |
|
Definition at line 27 of file tfim.c. Referenced by TFIM_getopts(). |
|
Definition at line 28 of file tfim.c. Referenced by main(), and TFIM_getopts(). |
|
Definition at line 29 of file tfim.c. Referenced by TFIM_getopts(). |
|
Definition at line 21 of file tfim.c. Referenced by main(), and TFIM_getopts(). |
|
Definition at line 26 of file tfim.c. Referenced by main(), and TFIM_getopts(). |
|
|
|
|
|
Definition at line 25 of file tfim.c. Referenced by main(), and TFIM_getopts(). |