/*****************************************************************************
Major portions of this software are copyrighted by the Medical College
of Wisconsin, 1994-2001, and are released under the Gnu General Public
License, Version 2. See the file README.Copyright for details.
******************************************************************************/
/*
This program calculates the single-factor analysis of variance (ANOVA)
for 3 dimensional AFNI data sets.
File: 3dANOVA.c
Author: B. D. Ward
Date: 09 December 1996
Mod: Incorporated include file 3dANOVA.h.
Date: 15 January 1997
Mod: Added option to check for required disk space.
Date: 23 January 1997
Mod: Added protection against divide by zero.
Date: 13 November 1997
Mod: Extensive changes required to implement the 'bucket' dataset.
Date: 30 December 1997
Mod: Separated 3dANOVA.h and 3dANOVA.lib files.
Date: 5 January 1998
Mod: Continuation of 'bucket' dataset changes.
Date: 9 January 1998
Mod: Added software for statistical tests of individual cell means.
Date: 27 October 1998
Mod: Added changes for incorporating History notes.
Date: 09 September 1999
Mod: Replaced dataset input code with calls to THD_open_dataset,
to allow operator selection of individual sub-bricks for input.
Date: 03 December 1999
Mod: Added call to AFNI_logger.
Date: 15 August 2001
Mod: Modified routine write_afni_data of 3dANOVA.lib so that all output
subbricks will now have the scaled short integer format.
Date: 14 March 2002
Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
Date: 02 December 2002
Mod: Setup to use .1D dataset filenames on output if these are input.
Date: 14 March 2003 - RWCox
Mod: -help menu modified.
Date: 21 July 2005 - P Christidis
Mod: small -help correction
Date: 25 Nov 2005 [rickr]
Mod: Modified contrast t-stat computations, and added -old_method,
-OK, -assume_sph and -debug options.
Date: 08 Dec 2005 [rickr]
*/
/*---------------------------------------------------------------------------*/
#define PROGRAM_NAME "3dANOVA" /* name of this program */
#define PROGRAM_AUTHOR "B. Douglas Ward" /* program author */
#define PROGRAM_INITIAL "09 Dec 1996" /* date of initial program release */
#define PROGRAM_LATEST "21 Jul 2005" /* date of latest program revision */
/*---------------------------------------------------------------------------*/
#define SUFFIX ".3danova" /* suffix for temporary data files */
#include "3dANOVA.h"
#include "3dANOVA.lib"
/*---------------------------------------------------------------------------*/
/*
Routine to display 3dANOVA help menu.
*/
void display_help_menu()
{
printf
(
"This program performs single factor Analysis of Variance (ANOVA)\n"
"on 3D datasets\n"
"\n"
"---------------------------------------------------------------\n"
"\n"
"Usage:\n"
"-----\n"
"\n"
"3dANOVA\n"
" -levels r : r = number of factor levels\n"
"\n"
" -dset 1 filename : data set for factor level 1\n"
" . . .. . .\n"
" -dset 1 filename data set for factor level 1\n"
" . . .. . .\n"
" -dset r filename data set for factor level r\n"
" . . .. . .\n"
" -dset r filename data set for factor level r\n"
"\n"
" [-voxel num] : screen output for voxel # num\n"
"\n"
" [-diskspace] : print out disk space required for\n"
" program execution\n"
"\n"
" [-debug level] : request extra output\n"
"\n"
"The following commands generate individual AFNI 2-sub-brick datasets:\n"
" (In each case, output is written to the file with the specified\n"
" prefix file name.)\n"
"\n"
" [-ftr prefix] : F-statistic for treatment effect\n"
"\n"
" [-mean i prefix] : estimate of factor level i mean\n"
"\n"
" [-diff i j prefix] : difference between factor levels\n"
"\n"
" [-contr c1...cr prefix] : contrast in factor levels\n"
"\n"
"Modified ANOVA computation options: (December, 2005)\n"
"\n"
" ** For details, see %s\n"
"\n"
"[-old_method] request to perform ANOVA using the previous\n"
" functionality (requires -OK, also)\n"
"\n"
"[-OK] confirm you understand that contrasts that\n"
" do not sum to zero have inflated t-stats, and\n"
" contrasts that do sum to zero assume sphericity\n"
" (to be used with -old_method)\n"
"\n"
"[-assume_sph] assume sphericity (zero-sum contrasts, only)\n"
"\n"
" This allows use of the old_method for\n"
" computing contrasts which sum to zero (this\n"
" includes diffs, for instance). Any contrast\n"
" that does not sum to zero is invalid, and\n"
" cannot be used with this option (such as\n"
" ameans).\n"
"\n"
"The following command generates one AFNI 'bucket' type dataset:\n"
"\n"
" [-bucket prefix] : create one AFNI 'bucket' dataset whose\n"
" sub-bricks are obtained by\n"
" concatenating the above output files;\n"
" the output 'bucket' is written to file\n"
" with prefix file name\n"
"\n", ANOVA_MODS_LINK);
printf
(
"N.B.: For this program, the user must specify 1 and only 1 sub-brick\n"
" with each -dset command. That is, if an input dataset contains\n"
" more than 1 sub-brick, a sub-brick selector must be used,\n"
" e.g., -dset 2 'fred+orig[3]'\n"
);
printf
("\n"
"Example of 3dANOVA:\n"
"------------------\n"
"\n"
" Example is based on a study with one factor (independent variable)\n"
" called 'Pictures', with 3 levels:\n"
" (1) Faces, (2) Houses, and (3) Donuts\n"
"\n"
" The ANOVA is being conducted on the data of subjects Fred and Ethel:\n"
"\n"
" 3dANOVA -levels 3 \\\n"
" -dset 1 fred_Faces+tlrc \\\n"
" -dset 1 ethel_Faces+tlrc \\\n"
" \\\n"
" -dset 2 fred_Houses+tlrc \\\n"
" -dset 2 ethel_Houses+tlrc \\\n"
" \\\n"
" -dset 3 fred_Donuts+tlrc \\\n"
" -dset 3 ethel_Donuts+tlrc \\\n"
" \\\n"
" -ftr Pictures \\\n"
" -mean 1 Faces \\\n"
" -mean 2 Houses \\\n"
" -mean 3 Donuts \\\n"
" -diff 1 2 FvsH \\\n"
" -diff 2 3 HvsD \\\n"
" -diff 1 3 FvsD \\\n"
" -contr 1 1 -1 FHvsD \\\n"
" -contr -1 1 1 FvsHD \\\n"
" -contr 1 -1 1 FDvsH \\\n"
" -bucket fred_n_ethel_ANOVA\n"
);
printf("\n" MASTER_SHORTHELP_STRING );
printf("---------------------------------------------------\n"
"Also see HowTo#5 - Group Analysis on the AFNI website:\n"
"http://afni.nimh.nih.gov/pub/dist/HOWTO/howto/ht05_group/html/index.shtml\n"
"\n" );
exit(0);
}
/*---------------------------------------------------------------------------*/
/*
Routine to get user specified anova options.
*/
void get_options (int argc, char ** argv, anova_options * option_data)
{
int nopt = 1; /* input option argument counter */
int ival; /* integer input */
int i, j, k; /* factor level counter */
int nijk; /* number of data files in cell i */
float fval; /* float input */
THD_3dim_dataset * dset=NULL; /* test whether data set exists */
char message[MAX_NAME_LENGTH]; /* error message */
/*----- does user request help menu? -----*/
if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
/*----- add to program log -----*/
AFNI_logger (PROGRAM_NAME,argc,argv);
/*----- initialize the input options -----*/
initialize_options (option_data);
/*----- main loop over input options -----*/
while (nopt < argc)
{
/*----- allocate memory for storing data file names -----*/
if ((option_data->xname == NULL) && (option_data->a > 0))
{
option_data->xname =
(char *****) malloc (sizeof(char ****) * option_data->a);
for (i = 0; i < option_data->a; i++)
{
option_data->xname[i] =
(char ****) malloc (sizeof(char ***) * 1);
for (j = 0; j < 1; j++)
{
option_data->xname[i][j] =
(char ***) malloc (sizeof(char **) * 1);
for (k = 0; k < 1; k++)
{
option_data->xname[i][j][k] =
(char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
}
}
}
}
/*----- -diskspace -----*/
if( strncmp(argv[nopt],"-diskspace",5) == 0 )
{
option_data->diskspace = 1;
nopt++ ; continue ; /* go to next arg */
}
/*----- -debug level 8 Dec 2005 [rickr] -----*/
if( strncmp(argv[nopt],"-debug",4) == 0 ) {
if( ++nopt >= argc ) ANOVA_error("need a level after -debug!") ;
option_data->debug = atoi(argv[nopt]);
nopt++ ; continue ; /* go to next arg */
}
/*----- -datum type -----*/
if( strncmp(argv[nopt],"-datum",6) == 0 ){
if( ++nopt >= argc ) ANOVA_error("need an argument after -datum!") ;
if( strcmp(argv[nopt],"short") == 0 ){
option_data->datum = MRI_short ;
} else if( strcmp(argv[nopt],"float") == 0 ){
option_data->datum = MRI_float ;
} else {
char buf[256] ;
sprintf(buf,"-datum of type '%s' is not supported in 3dANOVA!",
argv[nopt] ) ;
ANOVA_error(buf) ;
}
nopt++ ; continue ; /* go to next arg */
}
/*------------------------------------------------------------*/
/*----- Using the old_method: 08 Dec 2005 [rickr] -----*/
/* if -old_method
if -OK, all contrasts are okay
else if -assume_sph, contrasts adding to 0 are okay
else complain and fail
bits: -old_method = 001, -OK = 010, -assume_sph = 100
valid bit patterns:
000 - use the new method
011 - use the old method
101 - use the old method (only allows zero-sum contrasts)
------------------------------------------------------------*/
/*----- -old_method 23 Nov 2005 [rickr] -----*/
if (strncmp(argv[nopt], "-old_method", 6) == 0)
{
option_data->old_method |= 1;
nopt++;
continue;
}
/*----- -OK: denote both OK and old_method by old_method = 3 -----*/
if (strncmp(argv[nopt], "-OK", 3) == 0)
{
option_data->old_method |= 2;
nopt++;
continue;
}
/*----- -assume_sph: denote assume_sphericity by old_method = 4 ----*/
if (strncmp(argv[nopt], "-assume_sph", 11) == 0)
{
option_data->old_method |= 5; /* also set -old_method bit */
nopt++;
continue;
}
/*------- end old_method checks ------------------------------*/
/*------------------------------------------------------------*/
/*----- -session dirname -----*/
if( strncmp(argv[nopt],"-session",6) == 0 ){
nopt++ ;
if( nopt >= argc ) ANOVA_error("need argument after -session!") ;
strcpy(option_data->session , argv[nopt++]) ;
continue ;
}
/*----- -voxel num -----*/
if (strncmp(argv[nopt], "-voxel", 6) == 0)
{
nopt++;
if (nopt >= argc) ANOVA_error ("need argument after -voxel ");
sscanf (argv[nopt], "%d", &ival);
if (ival <= 0)
ANOVA_error ("illegal argument after -voxel ");
option_data->nvoxel = ival;
nopt++;
continue;
}
/*----- -levels r -----*/
if (strncmp(argv[nopt], "-levels", 7) == 0)
{
nopt++;
if (nopt >= argc) ANOVA_error ("need argument after -levels ");
sscanf (argv[nopt], "%d", &ival);
if ((ival <= 0) || (ival > MAX_LEVELS))
ANOVA_error ("illegal argument after -levels ");
option_data->a = ival;
nopt++;
continue;
}
/*----- -dset level filename -----*/
if (strncmp(argv[nopt], "-dset", 5) == 0)
{
nopt++;
if (nopt+1 >= argc) ANOVA_error ("need 2 arguments after -dset ");
sscanf (argv[nopt], "%d", &ival);
if ((ival <= 0) || (ival > option_data->a))
ANOVA_error ("illegal argument after -dset ");
option_data->na[ival-1] += 1;
nijk = option_data->na[ival-1];
if (nijk > MAX_OBSERVATIONS)
ANOVA_error ("too many data files");
/*--- check whether input files exist ---*/
nopt++;
dset = THD_open_dataset( argv[nopt] ); CHECK_OPEN_ERROR(dset,argv[nopt]);
/*--- check number of selected sub-bricks ---*/
if (DSET_NVALS(dset) != 1)
{
sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
argv[nopt]);
ANOVA_error (message);
}
THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
option_data->xname[ival-1][0][0][nijk-1]
= malloc (sizeof(char) * MAX_NAME_LENGTH);
strcpy (option_data->xname[ival-1][0][0][nijk-1],
argv[nopt]);
nopt++;
continue;
}
/*----- -ftr filename -----*/
if (strncmp(argv[nopt], "-ftr", 4) == 0)
{
nopt++;
if (nopt >= argc) ANOVA_error ("need argument after -ftr ");
option_data->nftr = 1;
option_data->ftrname = malloc (sizeof(char) * MAX_NAME_LENGTH);
strcpy (option_data->ftrname, argv[nopt]);
nopt++;
continue;
}
/*----- -mean level filename -----*/
if (strncmp(argv[nopt], "-mean", 5) == 0)
{
nopt++;
if (nopt+1 >= argc) ANOVA_error ("need 2 arguments after -mean ");
option_data->num_ameans++;
if (option_data->num_ameans > MAX_MEANS)
ANOVA_error ("too many factor level mean estimates");
sscanf (argv[nopt], "%d", &ival);
if ((ival <= 0) || (ival > option_data->a))
ANOVA_error ("illegal argument after -mean ");
option_data->ameans[option_data->num_ameans-1] = ival - 1;
nopt++;
option_data->amname[option_data->num_ameans-1]
= malloc (sizeof(char) * MAX_NAME_LENGTH);
strcpy (option_data->amname[option_data->num_ameans-1], argv[nopt]);
nopt++;
continue;
}
/*----- -diff level1 level2 filename -----*/
if (strncmp(argv[nopt], "-diff", 5) == 0)
{
nopt++;
if (nopt+2 >= argc) ANOVA_error ("need 3 arguments after -diff ");
option_data->num_adiffs++;
if (option_data->num_adiffs > MAX_DIFFS)
ANOVA_error ("too many factor level differences");
sscanf (argv[nopt], "%d", &ival);
if ((ival <= 0) || (ival > option_data->a))
ANOVA_error ("illegal argument after -diff ");
option_data->adiffs[option_data->num_adiffs-1][0] = ival - 1;
nopt++;
sscanf (argv[nopt], "%d", &ival);
if ((ival <= 0) || (ival > option_data->a))
ANOVA_error ("illegal argument after -diff ");
option_data->adiffs[option_data->num_adiffs-1][1] = ival - 1;
nopt++;
option_data->adname[option_data->num_adiffs-1]
= malloc (sizeof(char) * MAX_NAME_LENGTH);
strcpy (option_data->adname[option_data->num_adiffs-1], argv[nopt]);
nopt++;
continue;
}
/*----- -contr c1 ... cr filename -----*/
if (strncmp(argv[nopt], "-contr", 6) == 0)
{
nopt++;
if (nopt + option_data->a >= argc)
ANOVA_error ("need r+1 arguments after -contr ");
option_data->num_acontr++;
if (option_data->num_acontr > MAX_CONTR)
ANOVA_error ("too many factor level contrasts");
for (i = 0; i < option_data->a; i++)
{
sscanf (argv[nopt], "%f", &fval);
option_data->acontr[option_data->num_acontr - 1][i] = fval ;
nopt++;
}
option_data->acname[option_data->num_acontr-1]
= malloc (sizeof(char) * MAX_NAME_LENGTH);
strcpy (option_data->acname[option_data->num_acontr-1], argv[nopt]);
nopt++;
continue;
}
/*----- -bucket filename -----*/
if (strncmp(argv[nopt], "-bucket", 4) == 0)
{
nopt++;
if (nopt >= argc) ANOVA_error ("need argument after -bucket ");
option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
strcpy (option_data->bucket_filename, argv[nopt]);
nopt++;
continue;
}
/*----- unknown command -----*/
sprintf (message,"Unrecognized command line option: %s\n", argv[nopt]);
ANOVA_error (message);
}
if (option_data->old_method == 1)
ANOVA_error("-old_method is insufficient by itself");
}
/*---------------------------------------------------------------------------*/
/*
Routine to check whether temporary files already exist.
*/
void check_temporary_files (anova_options * option_data)
{
char filename[MAX_NAME_LENGTH]; /* temporary file name */
int i; /* file counter */
for (i = 0; i < option_data->a; i++)
{
/*----- temporary file name -----*/
sprintf (filename, "y.%d", i);
/*----- see if file already exists -----*/
check_one_temporary_file (filename);
}
check_one_temporary_file ("ysum");
check_one_temporary_file ("ss");
check_one_temporary_file ("ssto");
check_one_temporary_file ("sstr");
check_one_temporary_file ("sse");
}
/*---------------------------------------------------------------------------*/
/*
Routine to check whether output files already exist.
*/
void check_output_files (anova_options * option_data)
{
int i; /* index */
if (option_data->nftr > 0)
check_one_output_file (option_data, option_data->ftrname);
if (option_data->num_ameans > 0)
for (i = 0; i < option_data->num_ameans; i++)
check_one_output_file (option_data, option_data->amname[i]);
if (option_data->num_adiffs > 0)
for (i = 0; i < option_data->num_adiffs; i++)
check_one_output_file (option_data, option_data->adname[i]);
if (option_data->num_acontr > 0)
for (i = 0; i < option_data->num_acontr; i++)
check_one_output_file (option_data, option_data->acname[i]);
if (option_data->bucket_filename != NULL)
check_one_output_file (option_data, option_data->bucket_filename);
}
/*---------------------------------------------------------------------------*/
/*
Routine to check for valid inputs.
*/
void check_for_valid_inputs (anova_options * option_data)
{
int i; /* factor level index */
char message[MAX_NAME_LENGTH]; /* error message */
if (option_data->a < 2)
ANOVA_error ("must specify number of factor levels (r>1) ");
if (option_data->nt < option_data->a + 1)
ANOVA_error ("too few data sets for ANOVA");
for (i = 0; i < option_data->a; i++)
if (option_data->na[i] == 0)
{
sprintf (message, "level %d has too few data sets", i+1);
ANOVA_error (message);
}
/* check contrasts (show error, and specify 3dANOVA) */
if ( !contrasts_are_valid(option_data, 1, 1) )
ANOVA_error("invalid contrast(s)\n");
}
/*---------------------------------------------------------------------------*/
/*
Routine to calculate the number of data files that have to be stored.
Modified to account for 'bucket' dataset output.
*/
int required_data_files (anova_options * option_data)
{
int r; /* number of factor levels */
int nftr; /* number of F-treatment output files
(0 or 1) */
int nmeans; /* number of estimated mean output files */
int ndiffs; /* number of difference output files */
int ncontr; /* number of contrast output files */
int nout; /* number of output files */
int nmax; /* maximum number of disk files */
/*----- initialize local variables -----*/
r = option_data->a;
nftr = option_data->nftr;
nmeans = option_data->num_ameans;
ndiffs = option_data->num_adiffs;
ncontr = option_data->num_acontr;
nout = nftr + nmeans + ndiffs + ncontr;
nmax = r + 2 + nout;
if (option_data->bucket_filename != NULL)
nmax = max (nmax, 2*nout);
return (nmax);
}
/*---------------------------------------------------------------------------*/
/*
Routine to perform all ANOVA initialization.
*/
void initialize (int argc, char ** argv, anova_options ** option_data)
{
int i; /* factor level index */
/*----- save command line for history notes -----*/
commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
/*----- allocate memory space for input data -----*/
*option_data = (anova_options *) malloc(sizeof(anova_options));
if (*option_data == NULL)
ANOVA_error ("memory allocation error");
/*----- get command line inputs -----*/
get_options(argc, argv, *option_data);
/*----- use first data set to get data set dimensions -----*/
(*option_data)->first_dataset = (*option_data)->xname[0][0][0][0];
get_dimensions (*option_data);
printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
(*option_data)->nx, (*option_data)->ny,
(*option_data)->nz, (*option_data)->nxyz);
if ((*option_data)->nvoxel > (*option_data)->nxyz)
ANOVA_error ("argument of -voxel is too large");
/*----- count number of observations per treatment level -----*/
(*option_data)->nt = 0;
for (i = 0; i < (*option_data)->a; i++)
(*option_data)->nt += (*option_data)->na[i];
/*----- check for valid inputs -----*/
check_for_valid_inputs (*option_data);
/*----- check whether temporary files already exist -----*/
check_temporary_files (*option_data);
/*----- check whether output files already exist -----*/
if( THD_deathcon() ) check_output_files (*option_data);
/*----- check whether there is sufficient disk space -----*/
if ((*option_data)->diskspace) check_disk_space (*option_data);
}
/*---------------------------------------------------------------------------*/
/* 8 Dec 2005 [rickr]
routine to compute the degrees of freedom from a contrast
df = sum_i_in_contr( n_i - 1 )
*/
int df_for_contr(anova_options * option_data, float * contr)
{
int c, df = 0;
for (c = 0; c < option_data->a; c++)
if (contr[c] != 0.0)
df += (option_data->na[c] - 1);
return df;
}
/*---------------------------------------------------------------------------*/
/* 7 Dec 2005 [rickr,gangc]
routine to compute the mean and tstat of an ANOVA contrast
t = L / sqrt(1/df * SL2 * sum_c2)
L = sum_i(contr[i] * Y_i...)
SL = sum_i[ step(|contr[i]|) * ( sum_j(Y_ij^2) - n_i*Y_i.^2 ) ]
sum_c2 = sum_i[ contr[i]^2/n_i ]
note: Y_i. is considered the mean Y_i here
For accuracy, sums are computed using doubles, then copied to float.
The contrast is passed in to allow for more uses of this function.
*/
void calc_contr(anova_options * option_data, float * contr,
float * mean, /* save memory: use to read files */
float * tmean /* save memory: use to sum over obs. */
)
{
const float EPSILON = 1.0e-15; /* protect against divide by zero */
double * L; /* cumulative contrast mean */
double * S1; /* sum(n_i*Y_i^2), for mean Y_i */
double * S2; /* sum(Y_ij^2) {SL = S2-S1} */
double sum_c2; /* sum_i(c_i^2/n_i) */
double denom; /* for computations */
int i, j; /* indices for levels of factors A and C */
int a; /* number of levels */
int n; /* number of observations per cell */
int df; /* degrees of freedom */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
/*----- initialize local variables -----*/
a = option_data->a;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
L = (double *) malloc(sizeof(double)*nxyz);
S1 = (double *) malloc(sizeof(double)*nxyz);
S2 = (double *) malloc(sizeof(double)*nxyz);
if (!L || !S1 || !S2)
ANOVA_error ("calc_contr: unable to allocate sufficient memory");
for (ixyz = 0; ixyz < nxyz; ixyz++) /* init cumulative sums */
L[ixyz] = S1[ixyz] = S2[ixyz] = 0.0;
sum_c2 = 0.0;
if (option_data->debug > 1){
fprintf(stderr,"-d contr (%d levels) = ",a);
for ( i = 0; i < a; i++ ) fprintf(stderr,"%f ",contr[i]);
fprintf(stderr,"\n-d observations = ");
for ( i = 0; i < a; i++ ) fprintf(stderr,"%d ",option_data->na[i]);
fputc('\n',stderr);
}
/*----- loop over factor levels -----*/
for ( i = 0; i < a; i++ )
{
if (contr[i] == 0.0 ) continue; /* then skip this index */
n = option_data->na[i]; /* num observations */
for (ixyz = 0; ixyz < nxyz; ixyz++) /* clear sum over observations */
tmean[ixyz] = 0.0;
/*----- process observations within treatment level -----*/
for (j = 0; j < n; j++)
{
read_afni_data (option_data, option_data->xname[i][0][0][j], mean);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
tmean[ixyz] += mean[ixyz]; /* sum over obs. */
S2[ixyz] += (double)mean[ixyz]*mean[ixyz]; /* sum of squares */
}
}
/*----- tally sum of contrast and squares for the current k -----*/
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
S1[ixyz] += (double)tmean[ixyz]*tmean[ixyz]/n;/* is n_i*mean_Y_i^2 */
L[ixyz] += (double)contr[i]*tmean[ixyz]/n; /* cum. contrast mean */
}
if (option_data->debug > 1 && nvoxel > 0)
fprintf(stderr,"+d mean[%d] = %f, cmean = %f\n",
i,tmean[nvoxel-1]/n,L[nvoxel-1]);
}
/* get df and sum_c2 (seems more readable to put in separate loop) */
df = 0;
for ( i = 0; i < a; i++ )
{
if (contr[i] == 0.0 ) continue; /* then skip this index */
n = option_data->na[i]; /* num observations */
sum_c2 += (double)contr[i]*contr[i]/n; /* see above */
df += n - 1; /* degrees of freedom */
}
/*----- copy final results to float output -----*/
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
denom = sqrt((S2[ixyz] - S1[ixyz]) * sum_c2 / df); /* t denominator */
mean [ixyz] = L[ixyz];
tmean[ixyz] = (denom < EPSILON) ? 0.0 : L[ixyz] / denom;
}
if (option_data->debug > 1 && nvoxel > 0)
fprintf(stderr,"+d S1, S2, sum_c2, df = %f, %f, %f, %d\n",
S1[nvoxel-1], S2[nvoxel-1], sum_c2, df);
/*----- fly, be free! {thud} free! {thud} -----*/
free(L); free(S1); free (S2);
}
/*---------------------------------------------------------------------------*/
/*
Routine to sum over the observations within one treatment level.
Output is stored (temporarily) in data file y"i".3danova, where
"i" = treatment level.
*/
void calculate_y (anova_options * option_data)
{
float * x = NULL; /* pointer to original input data */
float * y = NULL; /* pointer to sum within treatment data */
int i; /* factor level index */
int j; /* observation number index */
int r; /* number of factor levels (treatments) */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
char filename[MAX_NAME_LENGTH]; /* name of file for sum within treatment */
/*----- initialize local variables -----*/
r = option_data->a;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
x = (float *) malloc(sizeof(float)*nxyz);
y = (float *) malloc(sizeof(float)*nxyz);
if ((x == NULL) || (y == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- loop over treatment levels -----*/
for (i = 0; i < r; i++)
{
/*----- sum observations within a treatment level -----*/
volume_zero (y, nxyz);
for (j = 0; j < option_data->na[i]; j++)
{
read_afni_data (option_data,
option_data->xname[i][0][0][j], x);
if (nvoxel > 0)
printf ("y[%d][%d] = %f \n", i+1, j+1, x[nvoxel-1]);
for (ixyz = 0; ixyz < nxyz; ixyz++)
y[ixyz] += x[ixyz];
}
/*----- save the sum for this treatment level -----*/
if (nvoxel > 0)
printf ("y[%d]. = %f \n", i+1, y[nvoxel-1]);
sprintf (filename, "y.%d", i);
volume_write (filename, y, nxyz);
}
/*----- release memory -----*/
free (y); y = NULL;
free (x); x = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to sum over all observations at all treatment levels.
The sum is stored (temporarily) in disk file ysum.3danova.
*/
void calculate_ysum (anova_options * option_data)
{
float * y = NULL; /* pointer to sum within treatment data */
float * ysum = NULL; /* pointer to overall sum data */
int i; /* factor level index */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
char filename[MAX_NAME_LENGTH]; /* sum within treatment input file name */
/*----- initialize local variables -----*/
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
y = (float *) malloc(sizeof(float)*nxyz);
ysum = (float *) malloc(sizeof(float)*nxyz);
if ((y == NULL) || (ysum == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- sum over all observations -----*/
volume_zero (ysum, nxyz);
for (i = 0; i < option_data->a; i++)
{
sprintf (filename, "y.%d", i);
volume_read (filename, y, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
ysum[ixyz] += y[ixyz];
}
if (nvoxel > 0)
printf ("y.. = %f \n", ysum[nvoxel-1]);
volume_write ("ysum",ysum, nxyz);
/*----- release memory -----*/
free (ysum); ysum = NULL;
free (y); y = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to calculate the sum of squares of all observations (SS).
The output is stored (temporarily) in disk file ss.3danova.
*/
void calculate_ss (anova_options * option_data)
{
float * x = NULL; /* pointer to original input data */
float * ss = NULL; /* pointer to sum of squares data */
int i; /* factor level index */
int j; /* observation data index */
int r; /* number of factor levels */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
float xval; /* float input data */
/*----- initialize local variables -----*/
r = option_data->a;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
x = (float *) malloc(sizeof(float)*nxyz);
ss = (float *) malloc(sizeof(float)*nxyz);
if ((x == NULL) || (ss == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- sum the squares of all observations -----*/
volume_zero (ss, nxyz);
for (i = 0; i < r; i++)
for (j = 0; j < option_data->na[i]; j++)
{
read_afni_data (option_data,
option_data->xname[i][0][0][j], x);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
xval = x[ixyz];
ss[ixyz] += xval * xval;
}
}
if (nvoxel > 0)
printf ("SS = %f \n", ss[nvoxel-1]);
volume_write ("ss", ss, nxyz);
/*----- release memory -----*/
free (ss); ss = NULL;
free (x); x = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to calculate total (corrected for the mean) sum of squares (SSTO).
The output is stored (temporarily) in disk file ssto.3danova.
*/
void calculate_ssto (anova_options * option_data)
{
float * ysum = NULL; /* ptr to sum over all observations */
float * ssto = NULL; /* pointer to ssto data */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
int nt; /* total number of observations */
float yval; /* sum over all observations */
/*----- initialize local variables -----*/
nt = option_data->nt;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
ysum = (float *) malloc(sizeof(float)*nxyz);
ssto = (float *) malloc(sizeof(float)*nxyz);
if ((ysum == NULL) || (ssto == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- SSTO = SS - ((YSUM)^2)/nt -----*/
volume_read ("ss", ssto, nxyz);
volume_read ("ysum", ysum, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
yval = ysum[ixyz];
ssto[ixyz] -= yval * yval / nt;
}
/*----- ss data file is no longer needed -----*/
volume_delete ("ss");
/*----- save total (corrected) sum of squares -----*/
if (nvoxel > 0)
printf ("SSTO = %f \n", ssto[nvoxel-1]);
volume_write ("ssto", ssto, nxyz);
/*----- release memory -----*/
free (ssto); ssto = NULL;
free (ysum); ysum = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to calculate the sum of squares due to the treatment (SSTR).
The output is stored (temporarily) in file sstr.3danova.
*/
void calculate_sstr (anova_options * option_data)
{
float * y = NULL; /* input data pointer */
float * sstr = NULL; /* sstr data pointer */
int i; /* factor level index */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
int ni; /* number of observations at level i */
int nt; /* total number of observations */
float yval; /* temporary float value */
char filename[MAX_NAME_LENGTH]; /* input data file name */
/*----- assign local variables -----*/
nt = option_data->nt;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
y = (float *) malloc(sizeof(float)*nxyz);
sstr = (float *) malloc(sizeof(float)*nxyz);
if ((y == NULL) || (sstr == NULL))
ANOVA_error ("unable to allocate sufficient memory");
volume_zero (sstr, nxyz);
for (i = 0; i < option_data->a; i++)
{
ni = option_data->na[i];
sprintf (filename, "y.%d", i);
volume_read (filename, y, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
yval = y[ixyz];
sstr[ixyz] += yval * yval / ni;
}
}
volume_read ("ysum", y, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
yval = y[ixyz];
sstr[ixyz] -= yval * yval / nt;
/*----- protection against round-off error -----*/
if (sstr[ixyz] < 0.0) sstr[ixyz] = 0.0;
}
/*----- ysum data file is no longer needed -----*/
volume_delete ("ysum");
/*----- save treatment sum of squares -----*/
if (nvoxel > 0)
printf ("SSTR = %f \n", sstr[nvoxel-1]);
volume_write ("sstr", sstr, nxyz);
/*----- release memory -----*/
free (sstr); sstr = NULL;
free (y); y = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to calculate the error sum of squares, SSE = SSTO - SSTR .
The result is stored (temporarily) in disk file sse.3danova.
*/
void calculate_sse (anova_options * option_data)
{
float * sstr = NULL; /* pointer to sstr data */
float * sse = NULL; /* pointer to sse data */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
/*----- assign local variables -----*/
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
sstr = (float *) malloc(sizeof(float)*nxyz);
sse = (float *) malloc(sizeof(float)*nxyz);
if ((sstr == NULL) || (sse == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- read SSTO (total sum of squares) -----*/
volume_read ("ssto", sse, nxyz);
/*----- read SSTR (treatment sum of squares) -----*/
volume_read ("sstr", sstr, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
/*----- SSE = SSTO - SSTR -----*/
sse[ixyz] -= sstr[ixyz];
/*----- protection against round-off error -----*/
if (sse[ixyz] < 0.0) sse[ixyz] = 0.0;
}
/*----- ssto data file is no longer needed -----*/
volume_delete ("ssto");
/*----- save error sum of squares -----*/
if (nvoxel > 0)
printf ("SSE = %f \n", sse[nvoxel-1]);
volume_write ("sse", sse, nxyz);
/*----- release memory -----*/
free (sse); sse = NULL;
free (sstr); sstr = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to calculate the F-statistic for treatment, F = MSTR / MSE,
where MSTR = SSTR / (r-1),
and MSE = SSE / (n-r).
The output is stored as a 2 sub-brick AFNI data set. The first sub-brick
contains the square root of MSTR (mean sum of squares due to treatment),
and the second sub-brick contains the corresponding F-statistic.
*/
void calculate_f_statistic (anova_options * option_data)
{
const float EPSILON = 1.0e-10; /* protect against divide by zero */
float * fin = NULL; /* pointer to input float data */
float * mstr = NULL; /* pointer to MSTR data */
float * ftr = NULL; /* pointer to F due-to-treatment */
int r; /* number of factor levels */
int nt; /* total number of observations */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
float fval; /* float value used in calculations */
float mse; /* mean square error */
/*----- initialize local variables -----*/
r = option_data->a;
nt = option_data->nt;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
/*----- Note: the order in which memory allocations take place
is important! -----*/
ftr = (float *) malloc(sizeof(float)*nxyz);
mstr = (float *) malloc(sizeof(float)*nxyz);
fin = (float *) malloc(sizeof(float)*nxyz);
if ((fin == NULL) || (ftr == NULL) || (mstr == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- calculate mean SS due to treatments -----*/
volume_read ("sstr", fin, nxyz);
fval = 1.0 / (r - 1.0);
for (ixyz = 0; ixyz < nxyz; ixyz++)
mstr[ixyz] = fin[ixyz] * fval; /*--- MSTR = SSTR / (r-1) ---*/
if (nvoxel > 0)
printf ("MSTR = %f \n", mstr[nvoxel-1]);
/*----- calculate F-statistic FTR = MSTR / MSE -----*/
volume_read ("sse", fin, nxyz);
fval = 1.0 / (nt - r);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
mse = fin[ixyz] * fval; /*--- MSE = SSE / (nt-r) ---*/
if (mse < EPSILON)
ftr[ixyz] = 0.0;
else
ftr[ixyz] = mstr[ixyz] / mse; /*--- FTR = MSTR / MSE ---*/
}
if (nvoxel > 0)
printf ("FTR = %f \n", ftr[nvoxel-1]);
/*----- release memory -----*/
free (fin); fin = NULL;
/*----- write out afni data file -----*/
for (ixyz = 0; ixyz < nxyz; ixyz++)
mstr[ixyz] = sqrt(mstr[ixyz]); /*-- mstr now holds square root --*/
write_afni_data (option_data, option_data->ftrname, mstr, ftr, r-1, nt-r);
/*----- this data file is no longer needed -----*/
volume_delete ("sstr");
/*----- release memory -----*/
free (mstr); mstr = NULL;
free (ftr); ftr = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to calculate the mean treatment effect at the user specified
treatment level. The output is stored as a 2 sub-brick AFNI data set.
The first sub-brick contains the estimated mean at this treatment level,
and the second sub-brick contains the corresponding t-statistic.
*/
void calculate_means (anova_options * option_data)
{
const float EPSILON = 1.0e-10; /* protect against divide by zero */
float * fin = NULL; /* pointer to float input data */
float * mean = NULL; /* pointer to treatment mean data */
float * tmean = NULL; /* pointer to t-statistic data */
float * contr = NULL; /* for new_method contrast */
int imean; /* output mean option index */
int level; /* factor level index */
int ni; /* number obs. at factor level i */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
int r; /* number of factor levels */
int nt; /* total number of observations */
int df; /* degrees of freedom */
int num_means; /* number of user requested means */
float fval; /* for calculating std. dev. */
float stddev; /* est. std. dev. of factor mean */
char filename[MAX_NAME_LENGTH]; /* input data file name */
/*----- initialize local variables -----*/
r = option_data->a;
nt = option_data->nt;
num_means = option_data->num_ameans;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
mean = (float *) malloc(sizeof(float)*nxyz);
tmean = (float *) malloc(sizeof(float)*nxyz);
contr = (float *) malloc(sizeof(float)*r);
if ((mean == NULL) || (tmean == NULL) || (contr == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- loop over user specified treatment means -----*/
for (imean = 0; imean < num_means; imean++)
{
level = option_data->ameans[imean];
ni = option_data->na[level];
if (option_data->old_method) /* old method 8 Dec 2005 [rickr] */
{
df = nt-r;
/*----- allocate temporary memory space -----*/
fin = (float *) malloc(sizeof(float)*nxyz);
if (fin == NULL)
ANOVA_error ("unable to allocate sufficient memory");
/*----- estimate factor mean for this treatment level -----*/
sprintf (filename, "y.%d", level);
volume_read (filename, fin, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
mean[ixyz] = fin[ixyz] / ni;
/*----- divide by estimated standard deviation of factor mean -----*/
volume_read ("sse", fin, nxyz);
fval = (1.0 / (nt-r)) * (1.0 / ni);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
stddev = sqrt(fin[ixyz] * fval);
if (stddev < EPSILON)
tmean[ixyz] = 0.0;
else
tmean[ixyz] = mean[ixyz] / stddev;
}
/*----- release memory -----*/
free (fin); fin = NULL;
} else { /* new method using contrast (for all cases) */
for(ixyz = 0; ixyz < r; ixyz ++) contr[ixyz] = 0.0; /* clear contr */
contr[level] = 1.0;
calc_contr(option_data, contr, mean, tmean);
df = df_for_contr(option_data, contr);
}
if (nvoxel > 0){
printf ("Mean [%d] = %f \n", level+1, mean[nvoxel-1]);
printf ("t-Mean [%d] = %f, df = %d\n", level+1, tmean[nvoxel-1], df);
}
/*----- write out afni data file -----*/
write_afni_data (option_data, option_data->amname[imean],
mean, tmean, df, 0);
}
/*----- release memory -----*/
free (tmean); free (mean); free (contr);
}
/*---------------------------------------------------------------------------*/
/*
Routine to estimate the difference in the means between two user specified
treatment levels. The output is a 2 sub-brick AFNI data set. The first
sub-brick contains the estimated difference in the means. The second
sub-brick contains the corresponding t-statistic.
*/
void calculate_differences (anova_options * option_data)
{
const float EPSILON = 1.0e-10; /* protect against divide by zero */
float * fin = NULL; /* pointer to float input data */
float * diff = NULL; /* pointer to est. diff. in means */
float * tdiff = NULL; /* pointer to t-statistic data */
float * contr; /* for contrast in new method */
int r; /* number of factor levels */
int nt; /* total number of observations */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
int num_diffs; /* number of user requested diffs. */
int idiff; /* index for requested differences */
int i, j; /* factor level indices */
int ni, nj; /* number of obs. at levels i and j */
int df; /* degrees of freedom */
float fval; /* for calculating std. dev. */
float stddev; /* est. std. dev. of difference */
char filename[MAX_NAME_LENGTH]; /* input file name */
/*----- initialize local variables -----*/
r = option_data->a;
nt = option_data->nt;
num_diffs = option_data->num_adiffs;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
diff = (float *) malloc(sizeof(float)*nxyz);
tdiff = (float *) malloc(sizeof(float)*nxyz);
contr = (float *) malloc(sizeof(float)*r);
if ((diff == NULL) || (tdiff == NULL) || (contr == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- loop over user specified treatment differences -----*/
for (idiff = 0; idiff < num_diffs; idiff++)
{
if (option_data->old_method) /* old method 8 Dec 2005 [rickr] */
{
df = nt - r;
/*----- allocate temporary memory space -----*/
fin = (float *) malloc(sizeof(float)*nxyz);
if (fin == NULL) ANOVA_error ("unable to allocate sufficient memory");
/*----- read first treatment level mean -----*/
i = option_data->adiffs[idiff][0];
ni = option_data->na[i];
sprintf (filename, "y.%d", i);
volume_read (filename, fin, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
diff[ixyz] = fin[ixyz] / ni;
/*----- subtract second treatment level mean -----*/
j = option_data->adiffs[idiff][1];
nj = option_data->na[j];
sprintf (filename, "y.%d", j);
volume_read (filename, fin, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
diff[ixyz] -= fin[ixyz] / nj;
/*----- divide by estimated standard deviation of difference -----*/
volume_read ("sse", fin, nxyz);
fval = (1.0 / (nt-r)) * ((1.0/ni) + (1.0/nj));
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
stddev = sqrt (fin[ixyz] * fval);
if (stddev < EPSILON)
tdiff[ixyz] = 0.0;
else
tdiff[ixyz] = diff[ixyz] / stddev;
}
} else { /* new method */
for (ixyz = 0; ixyz < r; ixyz++ ) contr[ixyz] = 0.0; /* clear old */
contr[option_data->adiffs[idiff][0]] = 1.0;
contr[option_data->adiffs[idiff][1]] = -1.0; /* set difference */
calc_contr(option_data, contr, diff, tdiff);
df = df_for_contr(option_data, contr);
}
if (nvoxel > 0) {
printf ("Diff [%d] - [%d] = %f \n", i+1, j+1, diff[nvoxel-1]);
printf ("t-Diff [%d] - [%d] = %f, df = %d \n",
i+1, j+1, tdiff[nvoxel-1], df);
}
/*----- release memory -----*/
free (fin); fin = NULL;
/*----- write out afni data file -----*/
write_afni_data (option_data, option_data->adname[idiff],
diff, tdiff, df, 0);
}
/*----- release memory -----*/
free (tdiff); free (diff); free (contr);
}
/*---------------------------------------------------------------------------*/
/*
Routine to estimate a user specified contrast in treatment levels.
The output is stored as a 2 sub-brick AFNI data set. The first
sub-brick contains the estimated contrast. The second sub-brick contains
the corresponding t-statistic.
*/
void calculate_contrasts (anova_options * option_data)
{
const float EPSILON = 1.0e-10; /* protect against divide by zero */
float * fin = NULL; /* pointer to float input data */
float * contr = NULL; /* pointer to contrast estimate */
float * tcontr = NULL; /* pointer to t-statistic data */
int r; /* number of factor levels */
int nt; /* total number of observations */
int ixyz, nxyz; /* voxel counters */
int nvoxel; /* output voxel # */
int num_contr; /* number of user requested contrasts */
int icontr; /* index of user requested contrast */
int level; /* factor level index */
int ni; /* number of obs. at factor level i */
int df; /* degrees of freedom */
float fval; /* for calculating std. dev. */
float c; /* contrast coefficient */
float stddev; /* est. std. dev. of contrast */
char filename[MAX_NAME_LENGTH]; /* input data file name */
/*----- initialize local variables -----*/
r = option_data->a;
nt = option_data->nt;
num_contr = option_data->num_acontr;
nxyz = option_data->nxyz;
nvoxel = option_data->nvoxel;
/*----- allocate memory space for calculations -----*/
contr = (float *) malloc(sizeof(float)*nxyz);
tcontr = (float *) malloc(sizeof(float)*nxyz);
fin = (float *) malloc(sizeof(float)*nxyz);
if ((contr == NULL) || (tcontr == NULL) || (fin == NULL))
ANOVA_error ("unable to allocate sufficient memory");
/*----- loop over user specified constrasts -----*/
for (icontr = 0; icontr < num_contr; icontr++)
{
if (option_data->old_method) /* old method 08 Dec 2005 [rickr] */
{
df = nt - r;
volume_zero (contr, nxyz);
fval = 0.0;
for (level = 0; level < r; level++)
{
c = option_data->acontr[icontr][level];
if (c == 0.0) continue;
/*----- add c * treatment level mean to contrast -----*/
sprintf (filename, "y.%d", level);
volume_read (filename, fin, nxyz);
ni = option_data->na[level];
fval += c * c / ni;
for (ixyz = 0; ixyz < nxyz; ixyz++)
contr[ixyz] += c * fin[ixyz] / ni;
}
/*----- divide by estimated standard deviation of the contrast -----*/
volume_read ("sse", fin, nxyz);
for (ixyz = 0; ixyz < nxyz; ixyz++)
{
stddev = sqrt ((fin[ixyz] / (nt-r)) * fval);
if (stddev < EPSILON)
tcontr[ixyz] = 0.0;
else
tcontr[ixyz] = contr[ixyz] / stddev;
}
} else { /* new method */
calc_contr(option_data, option_data->acontr[icontr], contr, tcontr);
df = df_for_contr(option_data, option_data->acontr[icontr]);
}
if (nvoxel > 0){
printf ("Contr no.%d = %f \n", icontr+1, contr[nvoxel-1]);
printf ("t-Contr no.%d = %f, df = %d\n",icontr+1,tcontr[nvoxel-1],df);
}
/*----- write out afni data file -----*/
write_afni_data (option_data, option_data->acname[icontr],
contr, tcontr, df, 0);
}
/*----- release memory -----*/
free (tcontr); tcontr = NULL;
free (contr); contr = NULL;
free (fin); fin = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to perform single factor ANOVA.
*/
void calculate_anova (anova_options * option_data)
{
/*----- sum observations within treatment -----*/
calculate_y (option_data);
/*----- sum all observations -----*/
calculate_ysum (option_data);
/*----- calculate sum of squares -----*/
calculate_ss (option_data);
/*----- calculate total (corrected for the mean) sum of squares -----*/
calculate_ssto (option_data);
/*----- calculate treatment sum of squares -----*/
calculate_sstr (option_data);
/*----- calculate error sum of squares -----*/
calculate_sse (option_data);
}
/*---------------------------------------------------------------------------*/
/*
Routine to analyze the results from a single factor ANOVA.
*/
void analyze_results (anova_options * option_data)
{
/*----- calculate F-statistic for treatment effect -----*/
if (option_data->nftr > 0) calculate_f_statistic (option_data);
/*----- estimate factor level means -----*/
if (option_data->num_ameans > 0) calculate_means (option_data);
/*----- estimate differences in factor level means -----*/
if (option_data->num_adiffs > 0) calculate_differences (option_data);
/*----- calculate contrasts in factor levels -----*/
if (option_data->num_acontr > 0) calculate_contrasts (option_data);
}
/*---------------------------------------------------------------------------*/
/*
Routine to create an AFNI "bucket" output dataset.
*/
void create_bucket (anova_options * option_data)
{
char bucket_str[10000]; /* command line for program 3dbucket */
char refit_str[10000]; /* command line for program 3drefit */
THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
int i; /* file index */
int ibrick; /* sub-brick index number */
char str[100]; /* temporary character string */
/*----- read first dataset -----*/
dset = THD_open_dataset (option_data->first_dataset) ;
CHECK_OPEN_ERROR(dset,option_data->first_dataset) ;
if( DSET_IS_1D(dset) ) USE_1D_filenames(1) ; /* 14 Mar 2003 */
else if( DSET_IS_3D(dset) ) USE_1D_filenames(3) ; /* 21 Mar 2003 */
/*----- make an empty copy of this dataset -----*/
new_dset = EDIT_empty_copy( dset ) ;
THD_delete_3dim_dataset (dset , False); dset = NULL;
EDIT_dset_items (new_dset,
ADN_directory_name, option_data->session,
ADN_none);
/*----- begin command line for program 3dbucket -----*/
strcpy (bucket_str, "3dbucket ");
strcat (bucket_str, " -prefix ");
strcat (bucket_str, option_data->bucket_filename);
/*----- begin command line for program 3drefit -----*/
strcpy (refit_str, "3drefit ");
ibrick = -1;
/*----- make F-statistic sub-bricks -----*/
if (option_data->nftr != 0)
{
add_file_name (new_dset, option_data->ftrname, bucket_str);
ibrick++;
sprintf (str, " -sublabel %d %s:Inten ",
ibrick, option_data->ftrname);
strcat (refit_str, str);
ibrick++;
sprintf (str, " -sublabel %d %s:F-stat ",
ibrick, option_data->ftrname);
strcat (refit_str, str);
}
/*----- make factor level mean sub-bricks -----*/
if (option_data->num_ameans > 0)
for (i = 0; i < option_data->num_ameans; i++)
{
add_file_name (new_dset, option_data->amname[i], bucket_str);
ibrick++;
sprintf (str, " -sublabel %d %s:Mean ",
ibrick, option_data->amname[i]);
strcat (refit_str, str);
ibrick++;
sprintf (str, " -sublabel %d %s:t-stat ",
ibrick, option_data->amname[i]);
strcat (refit_str, str);
}
/*----- make difference in factor level means sub-bricks -----*/
if (option_data->num_adiffs > 0)
for (i = 0; i < option_data->num_adiffs; i++)
{
add_file_name (new_dset, option_data->adname[i], bucket_str);
ibrick++;
sprintf (str, " -sublabel %d %s:Diff ",
ibrick, option_data->adname[i]);
strcat (refit_str, str);
ibrick++;
sprintf (str, " -sublabel %d %s:t-stat ",
ibrick, option_data->adname[i]);
strcat (refit_str, str);
}
/*----- make contrast in factor level means sub-bricks -----*/
if (option_data->num_acontr > 0)
for (i = 0; i < option_data->num_acontr; i++)
{
add_file_name (new_dset, option_data->acname[i], bucket_str);
ibrick++;
sprintf (str, " -sublabel %d %s:Contr ",
ibrick, option_data->acname[i]);
strcat (refit_str, str);
ibrick++;
sprintf (str, " -sublabel %d %s:t-stat ",
ibrick, option_data->acname[i]);
strcat (refit_str, str);
}
/*----- invoke program 3dbucket to generate bucket type output dataset
by concatenating previous output files -----*/
printf("Writing `bucket' dataset ");
printf("into %s\n", option_data->bucket_filename);
fprintf(stderr,"RUNNING COMMAND: %s\n",bucket_str);
system (bucket_str);
/*----- invoke program 3drefit to label individual sub-bricks -----*/
add_file_name (new_dset, option_data->bucket_filename, refit_str);
fprintf(stderr,"RUNNING COMMAND: %s\n",refit_str);
system (refit_str);
/*----- release memory -----*/
THD_delete_3dim_dataset (new_dset , False); new_dset = NULL;
}
/*---------------------------------------------------------------------------*/
/*
Routine to release memory and remove any remaining temporary data files.
If the '-bucket' option has been used, create the bucket dataset and
dispose of the 2 sub-brick datasets.
*/
void terminate (anova_options * option_data)
{
int i;
THD_3dim_dataset * dset=NULL; /* input afni data set pointer */
THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */
char filename[MAX_NAME_LENGTH];
/*----- remove temporary data files -----*/
volume_delete ("sstr");
volume_delete ("sse");
for (i = 0; i < option_data->a; i++)
{
sprintf (filename, "y.%d", i);
volume_delete (filename);
}
/*----- create bucket dataset -----*/
if (option_data->bucket_filename != NULL)
create_bucket (option_data);
/*----- if 'bucket' datset was created, remove the individual 2-subbrick
data files -----*/
if (option_data->bucket_filename != NULL)
{
/*----- read first dataset -----*/
dset = THD_open_dataset (option_data->first_dataset) ;
CHECK_OPEN_ERROR(dset,option_data->first_dataset) ;
/*----- make an empty copy of this dataset -----*/
new_dset = EDIT_empty_copy (dset);
THD_delete_3dim_dataset (dset , False); dset = NULL;
EDIT_dset_items (new_dset,
ADN_directory_name, option_data->session,
ADN_none);
/*----- remove F-statistic data file -----*/
if (option_data->nftr != 0)
remove_dataset (new_dset, option_data->ftrname);
/*----- remove mean factor level data files -----*/
if (option_data->num_ameans > 0)
for (i = 0; i < option_data->num_ameans; i++)
remove_dataset (new_dset, option_data->amname[i]);
/*----- remove difference in factor levels data files -----*/
if (option_data->num_adiffs > 0)
for (i = 0; i < option_data->num_adiffs; i++)
remove_dataset (new_dset, option_data->adname[i]);
/*----- remove contrast in factor levels data files -----*/
if (option_data->num_acontr > 0)
for (i = 0; i < option_data->num_acontr; i++)
remove_dataset (new_dset, option_data->acname[i]);
THD_delete_3dim_dataset (new_dset , False); new_dset = NULL;
}
/*----- deallocate memory -----*/
destroy_anova_options (option_data);
}
/*---------------------------------------------------------------------------*/
/*
Single factor analysis of variance (ANOVA).
*/
int main (int argc, char ** argv)
{
anova_options * option_data = NULL;
#if 0
/*----- Identify software -----*/
printf ("\n\n");
printf ("Program: %s \n", PROGRAM_NAME);
printf ("Author: %s \n", PROGRAM_AUTHOR);
printf ("Initial Release: %s \n", PROGRAM_INITIAL);
printf ("Latest Revision: %s \n", PROGRAM_LATEST);
printf ("\n");
#endif
/*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
mainENTRY("3dANOVA main"); machdep(); PRINT_VERSION("3dANOVA"); AUTHOR(PROGRAM_AUTHOR);
{ int new_argc ; char ** new_argv ;
addto_args( argc , argv , &new_argc , &new_argv ) ;
if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
}
/*----- program initialization -----*/
initialize (argc, argv, &option_data);
/*----- warn user (after any help) -----*/
if( !option_data->old_method )
fprintf(stderr,"\n"
"** Changes have been made for 3dANOVA computations.\n"
" For details, please see:\n"
" %s\n\n", ANOVA_MODS_LINK);
/*----- calculate sums of squares -----*/
calculate_anova (option_data);
/*----- generate requested output -----*/
analyze_results (option_data);
/*----- terminate program -----*/
terminate (option_data);
free (option_data); option_data = NULL;
exit(0);
}
syntax highlighted by Code2HTML, v. 0.9.1