/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2000, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include "mrilib.h"
/*** global data (gasp, shock, horror) ***/
#define MAX_NAME 64
#define SUFFIX_DIFF "diff"
#define SUFFIX_TSPM "tspm"
#define SUFFIX_CORR "corr"
static char TF_pname[MAX_NAME] = "tfim." ;
static MRI_IMARR * TF_set1 = NULL ;
static MRI_IMARR * TF_set2 = NULL ;
static float TF_bval = 0.0 ;
static int TF_use_bval = 0 ;
static float TF_pthresh = 0.0 ;
static int TF_eqcorr = 0 ;
static float TF_eqval = 0.0 ;
static int TF_paired = 0 ;
/*** prototypes ***/
double qginv( double ) ; /* stat functions */
double stas4( double , double ) ;
double stinv( double , double ) ;
void TFIM_syntax( char * ) ;
void TFIM_getopts( int , char * argv[] ) ;
/*** actual program code ***/
int main( int argc , char * argv[] )
{
int kk , ii,nx,ny,npix , num1,num2 , zerout ;
MRI_IMAGE ** stat_ret ;
MRI_IMAGE * avim1 , * sdim1 , * avim2 , * sdim2 ;
float * avar1 , * sdar1 , * avar2 , * sdar2 ;
MRI_IMAGE * difim , * tspim , * corim , * dofim ;
float * difar , * tspar , * corar , * dofar ;
char name[256] ;
float sdmax , dofbar ;
/*----- read inputs -----*/
printf(
"MCW TFIM: t-tests on sets of functional images, by RW Cox\n") ;
machdep() ;
if( argc < 2 ) TFIM_syntax("try tfim -help for usage") ;
else if( strcmp(argv[1],"-help") == 0 ) TFIM_syntax(NULL) ;
TFIM_getopts( argc , argv ) ;
#ifdef TFIM_DEBUG
printf("prefix = %s\n",TF_pname) ;
if( TF_use_bval == 1 ){
printf("bval = %f\n",TF_bval) ;
} else {
printf("set1 has %d images\n",TF_set1->num) ;
}
printf("set2 has %d images\n",TF_set2->num) ;
#endif
/*----- form averages of images -----*/
nx = TF_set2->imarr[0]->nx ;
ny = TF_set2->imarr[0]->ny ; npix = nx * ny ;
if( TF_set1 != NULL ){
num1 = TF_set1->num ;
for( kk=0 ; kk < num1 ; kk++ ){
(void) mri_stat_seq( TF_set1->imarr[kk] ) ;
mri_free( TF_set1->imarr[kk] ) ;
}
stat_ret = mri_stat_seq( NULL ) ;
avim1 = stat_ret[0] ; avar1 = mri_data_pointer( avim1 ) ;
sdim1 = stat_ret[1] ; sdar1 = mri_data_pointer( sdim1 ) ;
}
num2 = TF_set2->num ;
for( kk=0 ; kk < num2 ; kk++ ){
(void) mri_stat_seq( TF_set2->imarr[kk] ) ;
mri_free( TF_set2->imarr[kk] ) ;
}
stat_ret = mri_stat_seq( NULL ) ;
avim2 = stat_ret[0] ; avar2 = mri_data_pointer( avim2 ) ;
sdim2 = stat_ret[1] ; sdar2 = mri_data_pointer( sdim2 ) ;
/*----- process set averages into statistics -----*/
difim = mri_new( nx,ny , MRI_float ) ; difar = mri_data_pointer(difim) ;
tspim = mri_new( nx,ny , MRI_float ) ; tspar = mri_data_pointer(tspim) ;
dofim = mri_new( nx,ny , MRI_float ) ; dofar = mri_data_pointer(dofim) ;
zerout = 0 ; /* count of pixels zero-ed out due to no variance */
if( TF_use_bval == 1 ){
float scl = 1.0 / sqrt((double) num2) ;
for( ii=0 ; ii < npix ; ii++ ){
difar[ii] = avar2[ii] - TF_bval ;
dofar[ii] = num2 - 1 ;
if( sdar2[ii] > 0 ){
tspar[ii] = difar[ii] / (scl * sdar2[ii]) ;
} else {
tspar[ii] = 0.0 ; zerout++ ;
}
}
} else {
float q1 , q2 ;
float n1i =1.0/num1 , n2i =1.0/num2 ,
n11i=1.0/(num1-1.0), n21i=1.0/(num2-1.0) ;
for( ii=0 ; ii < npix ; ii++ ){
difar[ii] = avar2[ii] - avar1[ii] ;
q1 = SQR(sdar1[ii]) * n1i ; /* qi = sigi^2 / numi */
q2 = SQR(sdar2[ii]) * n2i ;
if( q1 > 0 && q2 > 0 ){
tspar[ii] = difar[ii] / sqrt(q1+q2) ;
dofar[ii] = SQR(q1+q2) / ( q1*q1*n11i + q2*q2*n21i ) ;
} else {
tspar[ii] = 0.0 ; dofar[ii] = num1+num2-2 ; zerout++ ;
}
}
}
if( zerout > 0 ){
printf("** set %d pixels to zero due to 0 variance!\n",zerout) ;
}
/** threshold, if desired **/
if( TF_pthresh > 0.0 ){
float thr ;
if( TF_use_bval == 1 ){
dofbar = dofar[0] ;
thr = stinv( TF_pthresh , dofbar ) ;
printf("-- fixed t-threshold = %g\n",thr) ;
} else {
dofbar = 0.0 ;
for( ii=0 ; ii < npix ; ii++ ) dofbar += dofar[ii] ;
dofbar /= npix ;
thr = stinv( TF_pthresh , dofbar ) ;
printf("-- variable t-threshold about %g (dof mean = %g)\n",
thr,dofbar) ;
}
for( ii=0 ; ii < npix ; ii++ ){
if( TF_use_bval == -1 ) thr = stinv( TF_pthresh , dofar[ii] ) ;
if( fabs(tspar[ii]) < thr ) difar[ii] = 0.0 ;
}
}
/*----- write output images -----*/
sprintf( name , "%s%s" , TF_pname , SUFFIX_DIFF ) ;
printf("-- writing difference file %s\n",name) ;
mri_write( name , difim ) ;
sprintf( name , "%s%s" , TF_pname , SUFFIX_TSPM ) ;
printf("-- writing t-statistic file %s\n",name) ;
mri_write( name , tspim ) ;
if( TF_eqcorr ){
float pth , thr , cth , doff ;
corim = mri_new( nx , ny , MRI_float ) ;
corar = mri_data_pointer( corim ) ;
for( ii=0 ; ii < npix ; ii++ ){
doff = (TF_eqval > 0) ? TF_eqval : dofar[ii] ;
corar[ii] = tspar[ii] / sqrt( doff + SQR(tspar[ii]) ) ;
}
sprintf( name , "%s%s" , TF_pname , SUFFIX_CORR ) ;
printf("-- writing 'correlation' file %s\n",name) ;
mri_write( name , corim ) ;
pth = (TF_pthresh > 0) ? TF_pthresh : 0.001 ;
doff = (TF_eqval > 0) ? TF_eqval : dofbar ;
thr = stinv( pth , doff ) ;
cth = thr / sqrt( doff + SQR(thr) ) ;
printf("-- 'correlation' threshold for p=%g is %g\n", pth,cth ) ;
}
exit(0) ;
}
/*---------------------------------------------------------------------*/
void TFIM_getopts( int argc , char * argv[] )
{
int nopt = 1 , kk , nx,ny ;
/*--- scan argument list ---*/
while( nopt < argc ){
/** -paired **/
if( strncmp(argv[nopt],"-paired",5) == 0 ){
TF_paired = 1 ;
nopt++ ; continue ; /* skip to next arg */
}
/** -prefix pname **/
if( strncmp(argv[nopt],"-prefix",5) == 0 ){
if( ++nopt >= argc ) TFIM_syntax("-prefix needs a name!") ;
strcpy( TF_pname , argv[nopt] ) ;
kk = strlen(TF_pname) ;
if( TF_pname[kk-1] != '.' ){
TF_pname[kk] = '.' ;
TF_pname[kk+1] = '\0' ;
}
nopt++ ; continue ; /* skip to next arg */
}
/** -pthresh pval **/
if( strncmp(argv[nopt],"-pthresh",5) == 0 ){
char * ch ;
if( ++nopt >= argc ) TFIM_syntax("-pthresh needs a value!");
TF_pthresh = strtod( argv[nopt] , &ch ) ;
if( *ch != '\0' || TF_pthresh <= 0.0 || TF_pthresh > 0.99999 )
TFIM_syntax("value after -pthresh is illegal!") ;
nopt++ ; continue ; /* skip to next arg */
}
/** -eqcorr dval **/
if( strncmp(argv[nopt],"-eqcorr",5) == 0 ){
char * ch ;
if( ++nopt >= argc ) TFIM_syntax("-eqcorr needs a value!");
TF_eqval = strtod( argv[nopt] , &ch ) ;
if( *ch != '\0' || TF_eqval < 0.0 )
TFIM_syntax("value after -eqcorr is illegal!") ;
TF_eqcorr = 1 ;
nopt++ ; continue ; /* skip to next arg */
}
/** after this point, the options are no longer 'free floating' **/
/** -base1 bval **/
if( strncmp(argv[nopt],"-base1",5) == 0 ){
char * ch ;
if( ++nopt >= argc ) TFIM_syntax("-base1 needs a value!");
if( TF_use_bval == -1 ) TFIM_syntax("-base1 with -set1 illegal!");
TF_bval = strtod( argv[nopt] , &ch ) ;
if( *ch != '\0' ) TFIM_syntax("value after -base1 is illegal!") ;
TF_use_bval = 1 ;
nopt++ ; continue ; /* skip to next arg */
}
/** -set1 file file ... **/
if( strncmp(argv[nopt],"-set1",5) == 0 ){
if( TF_use_bval == 1 ) TFIM_syntax("-set1 with -base1 illegal!");
TF_use_bval = -1 ;
for( kk=nopt+1 ; kk < argc ; kk++ )
if( strncmp(argv[kk],"-set2",5) == 0 ) break ;
if( kk >= argc ) TFIM_syntax("-set1 not followed by -set2") ;
if( kk-1-nopt <= 0 ) TFIM_syntax("-set1 has no files after it") ;
TF_set1 = mri_read_many_files( kk-1-nopt , argv+(nopt+1) ) ;
if( TF_set1 == NULL )
TFIM_syntax("cannot continue without -set1 images") ;
nopt = kk ; continue ; /* skip to arg that matched -set2 */
}
/** -set2 file file ... */
if( strncmp(argv[nopt],"-set2",5) == 0 ){
if( ++nopt >= argc ) TFIM_syntax("-set2 has not files after it") ;
TF_set2 = mri_read_many_files( argc-nopt , argv+nopt ) ;
if( TF_set2 == NULL )
TFIM_syntax("cannot continue without -set2 images") ;
break ; /* end of possible inputs */
}
/** get to here is bad news! **/
fprintf(stderr,"*** can't interpret this option: %s\n",argv[nopt]) ;
TFIM_syntax("try tfim -help for usage details") ;
}
/*--- check arguments for OK-ositiness ---*/
if( TF_use_bval == -1 &&
( TF_set1 == NULL || TF_set1->num < 2 ) )
TFIM_syntax("-set1 has too few files in it!") ;
if( TF_set2 == NULL || TF_set2->num < 2 )
TFIM_syntax("-set2 has too few files in it!") ;
if( TF_use_bval == 1 && TF_paired == 1 )
TFIM_syntax("-paired and -base1 are mutually exclusive!") ;
if( TF_paired == 1 && TF_set1 == NULL )
TFIM_syntax("-paired requires presence of -set1!") ;
if( TF_paired == 1 && TF_set1->num != TF_set2->num ){
char str[256] ;
sprintf(str,"-paired requires equal size images sets,\n"
"but -set1 has %d images and -set2 has %d images" ,
TF_set1->num , TF_set2->num ) ;
TFIM_syntax(str) ;
}
/* check images for consistency */
nx = TF_set2->imarr[0]->nx ;
ny = TF_set2->imarr[0]->ny ;
for( kk=1 ; kk < TF_set2->num ; kk++ ){
if( nx != TF_set2->imarr[kk]->nx || ny != TF_set2->imarr[kk]->ny ){
fprintf(stderr,
"*** image %d in -set2 not conformant to image 0\n",
kk) ;
TFIM_syntax("cannot continue with images whose sizes differ!") ;
}
}
if( TF_set1 != NULL ){
for( kk=0 ; kk < TF_set1->num ; kk++ ){
if( nx != TF_set1->imarr[kk]->nx || ny != TF_set1->imarr[kk]->ny ){
fprintf(stderr,
"*** image %d in -set1 not conformant to image 0 in -set2\n",
kk) ;
TFIM_syntax("cannot continue with images whose sizes differ!") ;
}
}
}
/* if paired t-test, subtract set1 from set2
to convert it into the equivalent base level test vs. 0 */
if( TF_paired == 1 ){
MRI_IMAGE * im1 , * im2 ;
float * ar1 , * ar2 ;
int ii , npix = nx * ny ;
for( kk=0 ; kk < TF_set1->num ; kk++ ){
im1 = mri_to_float(TF_set1->imarr[kk]) ; mri_free(TF_set1->imarr[kk]) ;
im2 = mri_to_float(TF_set2->imarr[kk]) ; mri_free(TF_set2->imarr[kk]) ;
ar1 = mri_data_pointer(im1) ; ar2 = mri_data_pointer(im2) ;
for( ii=0 ; ii < npix ; ii++ ) ar2[ii] -= ar1[ii] ;
mri_free(im1) ; TF_set2->imarr[kk] = im2 ;
}
FREE_IMARR( TF_set1 ) ; TF_set1 = NULL ;
TF_use_bval = 1 ;
TF_bval = 0.0 ;
}
return ;
}
/*---------------------------------------------------------------------*/
void TFIM_syntax( char * str )
{
if( str != NULL ){
fprintf(stderr,"*** %s\n",str) ;
exit(-1) ;
}
printf(
"\n"
"Usage 1: tfim [options] -set1 image_files ... -set2 image_files ...\n"
"Usage 2: tfim [options] -base1 bval -set2 image_files ...\n"
"\n"
"In usage 1, the collection of images files after '-set1' and the\n"
"collection after '-set2' are averaged and differenced, and the\n"
"difference is tested for significance with a 2 sample Student t-test.\n"
"\n"
"In usage 2, the collection of image files after '-set2' is averaged\n"
"and then has the constant numerical value 'bval' subtracted, and the\n"
"difference is tested for significance with a 1 sample Student t-test.\n"
"\n"
"N.B.: The input images can be in the usual 'short' or 'byte'\n"
" formats, or in the floating point 'flim' format.\n"
"N.B.: If in either set of images, a given pixel has zero variance\n"
" (i.e., is constant), then the t-test is not performed.\n"
" In that pixel, the .tspm file will be zero.\n"
"\n"
"Options are:\n"
"\n"
" -prefix pname: 'pname' is used as the prefix for the output\n"
" filenames. The output image files are\n"
" + pname.diff = average of set2 minus average of set1\n"
" (or minus 'bval')\n"
" + pname.tspm = t-statistic of difference\n"
" Output images are in the 'flim' (floating pt. image)\n"
" format, and may be converted to 16 bit shorts using\n"
" the program 'ftosh'.\n"
" *** The default 'pname' is 'tfim', if -prefix isn't used.\n"
" -pthresh pval: 'pval' is a numeric value between 0 and 1, giving\n"
" the significance level (per voxel) to threshold the\n"
" output with; voxels with (2-sided) t-statistic\n"
" less significant than 'pval' will have their diff\n"
" output zeroed.\n"
" *** The default is no threshold, if -pthresh isn't used.\n"
" -eqcorr dval: If present, this option means to write out the file\n"
" pname.corr = equivalent correlation statistic\n"
" = t/sqrt(dof+t^2)\n"
" The number 'dval' is the value to use for 'dof' if\n"
" dval is positive. This would typically be the total\n"
" number of data images used in forming the image sets,\n"
" if the image sets are from sfim or fim.\n"
" If dval is zero, then dof is computed from the number\n"
" of images in -set1 and -set2; if these are averages\n"
" from program sfim, then dof will be smallish, which in\n"
" turn means that significant corr values will be higher\n"
" than you may be used to from using program fim.\n"
" *** The default is not to write, if -eqcorr isn't used.\n"
" -paired: If present, this means that -set1 and -set2 should be\n"
" compared using a paired sample t-test. This option is\n"
" illegal with the -base1 option. The number of samples\n"
" in the two sets of images must be equal.\n"
" [This test is implemented by subtracting -set1 images\n"
" from the -set2 images, then testing as in '-base1 0'.]\n"
" *** The default is to do an unpaired test, if -paired isn't\n"
" used. In that case, -set1 and -set2 don't need to have\n"
" the same number of images.\n"
) ;
exit(0) ;
}
/*----------------------------------------------------------------------
code for inverse of central t distribution
Inputs: p = double sided tail probability
nu = degrees of freedom
Output: T such that P( |t| > T ) = p
This version is only good for nu >= 5, since it uses the
approximations in Abramowitz and Stegun, Eq. 26.7.5 (p. 949).
------------------------------------------------------------------------*/
double stinv( double p , double nu )
{
double xg , t4 ;
xg = qginv(0.5*p) ;
t4 = stas4( xg , nu ) ;
return t4 ;
}
double stas4( double x , double nu) /* this code generated by Maple */
{
double t1,t2,t8,t9,t14,t17,t26,t34,t37 ;
t1 = x*x;
t2 = t1*x;
t8 = t1*t1;
t9 = t8*x;
t14 = nu*nu;
t17 = t8*t2;
t26 = t8*t8;
t34 = t14*t14;
t37 = x+(t2/4+x/4)/nu
+(5.0/96.0*t9+t2/6+x/32)/t14
+(t17/128+19.0/384.0*t9
+17.0/384.0*t2-5.0/128.0*x)/t14/nu
+(79.0/92160.0*t26*x+97.0/11520.0*t17+247.0/15360.0*t9
-t2/48-21.0/2048.0*x)/t34;
return t37 ;
}
/*** given p, return x such that Q(x)=p, for 0 < p < 1 ***/
double qginv( double p )
{
double dp , dx , dt , ddq , dq ;
int newt ;
dp = (p <= 0.5) ? (p) : (1.0-p) ; /* make between 0 and 0.5 */
if( dp <= 0.0 ){
dx = 13.0 ;
return ( (p <= 0.5) ? (dx) : (-dx) ) ;
}
/** Step 1: use 26.2.23 from Abramowitz and Stegun **/
dt = sqrt( -2.0 * log(dp) ) ;
dx = dt
- ((.010328e+0*dt + .802853e+0)*dt + 2.525517e+0)
/(((.001308e+0*dt + .189269e+0)*dt + 1.432788e+0)*dt + 1.e+0) ;
#if 0
/** Step 2: do 3 Newton steps to improve this **/
for( newt=0 ; newt < 3 ; newt++ ){
dq = 0.5e+0 * erfc( dx / 1.414213562373095e+0 ) - dp ;
ddq = exp( -0.5e+0 * dx * dx ) / 2.506628274631000e+0 ;
dx = dx + dq / ddq ;
}
#endif
return ( (p <= 0.5) ? (dx) : (-dx) ) ; /* return with correct sign */
}
syntax highlighted by Code2HTML, v. 0.9.1