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 }

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 }

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 }

00465 {
00466    double xg , t4 ;
00467    xg = qginv(0.5*p) ;
00468    t4 = stas4( xg , nu ) ;
00469    return t4 ;
00470 }

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 }

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 }