00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2001, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008    This program calculates the single-factor analysis of variance (ANOVA)
00009    for 3 dimensional AFNI data sets. 
00010 
00011    File:    3dANOVA.c
00012    Author:  B. D. Ward
00013    Date:    09 December 1996
00014 
00015    Mod:     Incorporated include file 3dANOVA.h.
00016    Date:    15 January 1997
00017 
00018    Mod:     Added option to check for required disk space.
00019    Date:    23 January 1997
00020 
00021    Mod:     Added protection against divide by zero.
00022    Date:    13 November 1997
00023 
00024    Mod:     Extensive changes required to implement the 'bucket' dataset.
00025    Date:    30 December 1997
00026 
00027    Mod:     Separated 3dANOVA.h and 3dANOVA.lib files.
00028    Date:    5 January 1998
00029 
00030    Mod:     Continuation of 'bucket' dataset changes.
00031    Date:    9 January 1998
00032 
00033    Mod:     Added software for statistical tests of individual cell means.
00034    Date:    27 October 1998
00035 
00036    Mod:     Added changes for incorporating History notes.
00037    Date:    09 September 1999
00038 
00039    Mod:     Replaced dataset input code with calls to THD_open_dataset,
00040             to allow operator selection of individual sub-bricks for input.
00041    Date:    03 December 1999
00042 
00043    Mod:     Added call to AFNI_logger.
00044    Date:    15 August 2001
00045 
00046    Mod:     Modified routine write_afni_data of 3dANOVA.lib so that all output
00047             subbricks will now have the scaled short integer format.
00048    Date:    14 March 2002
00049 
00050    Mod:     Set MAX_NAME_LENGTH equal to THD_MAX_NAME.
00051    Date:    02 December 2002
00052 
00053    Mod:     Setup to use .1D dataset filenames on output if these are input.
00054    Date:    14 March 2003 - RWCox
00055    
00056    Mod:     -help menu modified.
00057    Date:    21 July 2005 - P Christidis
00058 */
00059 
00060 /*---------------------------------------------------------------------------*/
00061 
00062 #define PROGRAM_NAME    "3dANOVA"                    /* name of this program */
00063 #define PROGRAM_AUTHOR  "B. Douglas Ward"                  /* program author */
00064 #define PROGRAM_INITIAL "09 Dec 1996"     /* date of initial program release */
00065 #define PROGRAM_LATEST  "21 Jul 2005"     /* date of latest program revision */
00066 
00067 /*---------------------------------------------------------------------------*/
00068 
00069 #define SUFFIX ".3danova"                 /* suffix for temporary data files */
00070 
00071 #include "3dANOVA.h"
00072 #include "3dANOVA.lib"
00073 
00074 
00075 /*---------------------------------------------------------------------------*/
00076 /*
00077    Routine to display 3dANOVA help menu.
00078 */
00079 
00080 void display_help_menu()
00081 {
00082   printf 
00083     (
00084     "This program performs single factor Analysis of Variance (ANOVA)      \n"
00085     "on 3D datasets                                                        \n"
00086     "                                                                      \n"
00087     "---------------------------------------------------------------       \n"
00088     "                                                                      \n"
00089     "Usage:                                                                \n"
00090     "-----                                                                 \n"
00091     "                                                                      \n"     
00092     "3dANOVA                                                               \n"
00093     "   -levels r                   : r = number of factor levels          \n"
00094     "                                                                      \n"
00095     "   -dset 1 filename            : data set for factor level 1          \n"
00096     "         . . .                            . . .                       \n"
00097     "   -dset 1 filename              data set for factor level 1          \n"
00098     "         . . .                            . . .                       \n"
00099     "   -dset r filename              data set for factor level r          \n"
00100     "         . . .                             . . .                      \n"
00101     "   -dset r filename              data set for factor level r          \n"
00102     "                                                                      \n"
00103     "  [-voxel num]                 : screen output for voxel # num        \n"
00104     "                                                                      \n"
00105     "  [-diskspace]                 : print out disk space required for    \n"
00106     "                                 program execution                    \n"
00107     "                                                                      \n"
00108     "The following commands generate individual AFNI 2-sub-brick datasets: \n"
00109     "  (In each case, output is written to the file with the specified     \n"
00110     "   prefix file name.)                                                 \n"
00111     "                                                                      \n"
00112     "  [-ftr prefix]                : F-statistic for treatment effect     \n"
00113     "                                                                      \n"  
00114     "  [-mean i prefix]             : estimate of factor level i mean      \n"
00115     "                                                                      \n"  
00116     "  [-diff i j prefix]           : difference between factor levels     \n"
00117     "                                                                      \n" 
00118     "  [-contr c1...cr prefix]      : contrast in factor levels            \n"
00119     "                                                                      \n"
00120     "The following command generates one AFNI 'bucket' type dataset:       \n"
00121     "                                                                      \n"
00122     "  [-bucket prefix]             : create one AFNI 'bucket' dataset whose \n"
00123     "                                 sub-bricks are obtained by             \n"
00124     "                                 concatenating the above output files;  \n"
00125     "                                 the output 'bucket' is written to file \n"
00126     "                                 with prefix file name                  \n"
00127     "\n");
00128 
00129   printf
00130     (
00131     "N.B.: For this program, the user must specify 1 and only 1 sub-brick  \n"
00132     "      with each -dset command. That is, if an input dataset contains  \n"
00133     "      more than 1 sub-brick, a sub-brick selector must be used,       \n"
00134     "      e.g., -dset 2 'fred+orig[3]'                                    \n"
00135      );
00136    
00137   printf
00138    ("\n"
00139     "Example of 3dANOVA:                                                   \n"
00140     "------------------                                                    \n"
00141     "                                                                      \n"
00142     " Example is based on a study with one factor (independent variable)   \n"
00143     " called 'Pictures', with 3 levels:                                    \n"
00144     "        (1) Faces, (2) Houses, and (3) Donuts                         \n"
00145     "                                                                      \n"
00146     " The ANOVA is being conducted on subject Fred's data:                 \n"
00147     "                                                                      \n"
00148     " 3dANOVA -levels 3                     \\                             \n"
00149     "         -dset 1 fred_Faces+tlrc       \\                             \n"
00150     "         -dset 2 fred_Houses+tlrc      \\                             \n"
00151     "         -dset 3 fred_Donuts+tlrc      \\                             \n"
00152     "         -ftr Pictures                 \\                             \n"
00153     "         -mean 1 Faces                 \\                             \n"
00154     "         -mean 2 Houses                \\                             \n"
00155     "         -mean 3 Donuts                \\                             \n"
00156     "         -diff 1 2 FvsH                \\                             \n"
00157     "         -diff 2 3 HvsD                \\                             \n"
00158     "         -diff 1 3 FvsD                \\                             \n"
00159     "         -contr  1  1 -1 FHvsD         \\                             \n"
00160     "         -contr -1  1  1 FvsHD         \\                             \n"
00161     "         -contr  1 -1  1 FDvsH         \\                             \n"
00162     "         -bucket fred_ANOVA                                           \n"
00163     ); 
00164  
00165   printf("\n" MASTER_SHORTHELP_STRING );
00166 
00167   printf("---------------------------------------------------\n"
00168    "Also see HowTo#5 - Group Analysis on the AFNI website:                \n"
00169    "http://afni.nimh.nih.gov/pub/dist/HOWTO/howto/ht05_group/html/index.shtml\n"
00170      "\n" );
00171      
00172   exit(0);
00173 }
00174 
00175 /*---------------------------------------------------------------------------*/
00176 /*
00177    Routine to get user specified anova options.
00178 */
00179 
00180 void get_options (int argc, char ** argv, anova_options * option_data)
00181 {
00182   int nopt = 1;                  /* input option argument counter */
00183   int ival;                      /* integer input */
00184   int i, j, k;                   /* factor level counter */                   
00185   int nijk;                      /* number of data files in cell i */     
00186   float fval;                    /* float input */
00187   THD_3dim_dataset * dset=NULL;             /* test whether data set exists */
00188   char message[MAX_NAME_LENGTH];            /* error message */
00189 
00190 
00191   /*----- does user request help menu? -----*/
00192   if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();
00193 
00194   /*----- add to program log -----*/
00195   AFNI_logger (PROGRAM_NAME,argc,argv); 
00196 
00197   /*----- initialize the input options -----*/
00198   initialize_options (option_data);
00199 
00200   
00201   /*----- main loop over input options -----*/
00202   while (nopt < argc)
00203     {
00204       
00205       /*----- allocate memory for storing data file names -----*/
00206       if ((option_data->xname == NULL) && (option_data->a > 0))
00207         {
00208           option_data->xname = 
00209             (char *****) malloc (sizeof(char ****) * option_data->a);
00210           for (i = 0;  i < option_data->a;  i++)
00211             {
00212               option_data->xname[i] =
00213                 (char ****) malloc (sizeof(char ***) * 1);
00214               for (j = 0;  j < 1;  j++)
00215                 {
00216                   option_data->xname[i][j] =
00217                     (char ***) malloc (sizeof(char **) * 1);
00218                   for (k = 0;  k < 1;  k++)
00219                     {
00220                       option_data->xname[i][j][k] =
00221                         (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00222                     }
00223                 }
00224             }
00225         }
00226 
00227 
00228       /*-----   -diskspace   -----*/
00229       if( strncmp(argv[nopt],"-diskspace",5) == 0 )
00230         {
00231           option_data->diskspace = 1;
00232           nopt++ ; continue ;  /* go to next arg */
00233         }
00234 
00235       
00236       /*-----   -datum type   -----*/
00237       if( strncmp(argv[nopt],"-datum",6) == 0 ){
00238         if( ++nopt >= argc ) ANOVA_error("need an argument after -datum!") ;
00239         
00240         if( strcmp(argv[nopt],"short") == 0 ){
00241           option_data->datum = MRI_short ;
00242         } else if( strcmp(argv[nopt],"float") == 0 ){
00243           option_data->datum = MRI_float ;
00244         } else {
00245           char buf[256] ;
00246           sprintf(buf,"-datum of type '%s' is not supported in 3dANOVA!",
00247                   argv[nopt] ) ;
00248           ANOVA_error(buf) ;
00249         }
00250         nopt++ ; continue ;  /* go to next arg */
00251       }
00252       
00253       
00254       /*-----   -session dirname    -----*/
00255       if( strncmp(argv[nopt],"-session",6) == 0 ){
00256         nopt++ ;
00257         if( nopt >= argc ) ANOVA_error("need argument after -session!") ;
00258         strcpy(option_data->session , argv[nopt++]) ;
00259         continue ;
00260       }
00261       
00262       
00263       /*-----   -voxel num  -----*/
00264       if (strncmp(argv[nopt], "-voxel", 6) == 0)
00265         {
00266           nopt++;
00267           if (nopt >= argc)  ANOVA_error ("need argument after -voxel ");
00268           sscanf (argv[nopt], "%d", &ival);
00269           if (ival <= 0)
00270             ANOVA_error ("illegal argument after -voxel ");
00271           option_data->nvoxel = ival;
00272           nopt++;
00273           continue;
00274         }
00275       
00276       
00277       /*-----   -levels r  -----*/
00278       if (strncmp(argv[nopt], "-levels", 7) == 0)
00279         {
00280           nopt++;
00281           if (nopt >= argc)  ANOVA_error ("need argument after -levels ");
00282           sscanf (argv[nopt], "%d", &ival);
00283           if ((ival <= 0) || (ival > MAX_LEVELS))
00284             ANOVA_error ("illegal argument after -levels ");
00285           option_data->a = ival;
00286           nopt++;
00287           continue;
00288         }
00289       
00290       
00291       /*-----   -dset level filename   -----*/
00292       if (strncmp(argv[nopt], "-dset", 5) == 0)
00293         {
00294           nopt++;
00295           if (nopt+1 >= argc)  ANOVA_error ("need 2 arguments after -dset ");
00296           sscanf (argv[nopt], "%d", &ival);
00297           if ((ival <= 0) || (ival > option_data->a))
00298             ANOVA_error ("illegal argument after -dset ");
00299           
00300           option_data->na[ival-1] += 1;
00301           nijk = option_data->na[ival-1];
00302           if (nijk > MAX_OBSERVATIONS)
00303             ANOVA_error ("too many data files");
00304           
00305           /*--- check whether input files exist ---*/
00306           nopt++;
00307           dset = THD_open_dataset( argv[nopt] ) ;
00308           if( ! ISVALID_3DIM_DATASET(dset) )
00309             {
00310              sprintf(message,"Unable to open dataset file %s\n", 
00311                      argv[nopt]);
00312              ANOVA_error (message);
00313             }
00314 
00315           /*--- check number of selected sub-bricks ---*/
00316           if (DSET_NVALS(dset) != 1)
00317             {
00318              sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
00319                      argv[nopt]);
00320              ANOVA_error (message);
00321             }
00322 
00323           THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00324           
00325           option_data->xname[ival-1][0][0][nijk-1] 
00326             =  malloc (sizeof(char) * MAX_NAME_LENGTH);
00327           strcpy (option_data->xname[ival-1][0][0][nijk-1], 
00328                   argv[nopt]);
00329 
00330           nopt++;
00331           continue;
00332         }
00333       
00334       
00335       /*-----   -ftr filename   -----*/
00336       if (strncmp(argv[nopt], "-ftr", 4) == 0)
00337         {
00338           nopt++;
00339           if (nopt >= argc)  ANOVA_error ("need argument after -ftr ");
00340           option_data->nftr = 1;
00341           option_data->ftrname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00342           strcpy (option_data->ftrname, argv[nopt]);
00343           nopt++;
00344           continue;
00345         }
00346       
00347       
00348       /*-----   -mean level filename   -----*/
00349       if (strncmp(argv[nopt], "-mean", 5) == 0)
00350         {
00351           nopt++;
00352           if (nopt+1 >= argc)  ANOVA_error ("need 2 arguments after -mean ");
00353           
00354           option_data->num_ameans++;
00355           if (option_data->num_ameans > MAX_MEANS)
00356             ANOVA_error ("too many factor level mean estimates");
00357           
00358           sscanf (argv[nopt], "%d", &ival);
00359           if ((ival <= 0) || (ival > option_data->a))
00360             ANOVA_error ("illegal argument after -mean ");
00361           option_data->ameans[option_data->num_ameans-1] = ival - 1;
00362           nopt++;
00363           
00364           option_data->amname[option_data->num_ameans-1] 
00365             =  malloc (sizeof(char) * MAX_NAME_LENGTH);
00366           strcpy (option_data->amname[option_data->num_ameans-1], argv[nopt]);
00367           nopt++;
00368           continue;
00369         }
00370       
00371       
00372       /*-----   -diff level1 level2 filename   -----*/
00373       if (strncmp(argv[nopt], "-diff", 5) == 0)
00374         {
00375           nopt++;
00376           if (nopt+2 >= argc)  ANOVA_error ("need 3 arguments after -diff ");
00377           
00378           option_data->num_adiffs++;
00379           if (option_data->num_adiffs > MAX_DIFFS)
00380             ANOVA_error ("too many factor level differences");
00381           
00382           sscanf (argv[nopt], "%d", &ival);
00383           if ((ival <= 0) || (ival > option_data->a))
00384             ANOVA_error ("illegal argument after -diff ");
00385           option_data->adiffs[option_data->num_adiffs-1][0] = ival - 1;
00386           nopt++;
00387           
00388           sscanf (argv[nopt], "%d", &ival);
00389           if ((ival <= 0) || (ival > option_data->a))
00390             ANOVA_error ("illegal argument after -diff ");
00391           option_data->adiffs[option_data->num_adiffs-1][1] = ival - 1;
00392           nopt++;
00393           
00394           option_data->adname[option_data->num_adiffs-1] 
00395             =  malloc (sizeof(char) * MAX_NAME_LENGTH);
00396           strcpy (option_data->adname[option_data->num_adiffs-1], argv[nopt]);
00397           nopt++;
00398           continue;
00399         }
00400       
00401       
00402       /*-----   -contr c1 ... cr filename   -----*/
00403       if (strncmp(argv[nopt], "-contr", 6) == 0)
00404         {
00405           nopt++;
00406           if (nopt + option_data->a >= argc)  
00407             ANOVA_error ("need r+1 arguments after -contr ");
00408           
00409           option_data->num_acontr++;
00410           if (option_data->num_acontr > MAX_CONTR)
00411             ANOVA_error ("too many factor level contrasts");
00412           
00413           for (i = 0;  i < option_data->a;  i++)
00414             {
00415               sscanf (argv[nopt], "%f", &fval); 
00416               option_data->acontr[option_data->num_acontr - 1][i] = fval ;
00417               nopt++;
00418             }
00419           
00420           option_data->acname[option_data->num_acontr-1] 
00421             =  malloc (sizeof(char) * MAX_NAME_LENGTH);
00422           strcpy (option_data->acname[option_data->num_acontr-1], argv[nopt]);
00423           nopt++;
00424           continue;
00425         }
00426       
00427       
00428       /*-----   -bucket filename   -----*/
00429       if (strncmp(argv[nopt], "-bucket", 4) == 0)
00430         {
00431           nopt++;
00432           if (nopt >= argc)  ANOVA_error ("need argument after -bucket ");
00433           option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00434           strcpy (option_data->bucket_filename, argv[nopt]);
00435           nopt++;
00436           continue;
00437         }
00438       
00439       
00440       /*----- unknown command -----*/
00441       sprintf (message,"Unrecognized command line option: %s\n", argv[nopt]);
00442       ANOVA_error (message);
00443     }
00444 
00445 }
00446 
00447 /*---------------------------------------------------------------------------*/
00448 /*
00449    Routine to check whether temporary files already exist.
00450 */
00451 
00452 void check_temporary_files (anova_options * option_data)
00453 {
00454    char filename[MAX_NAME_LENGTH];           /* temporary file name */ 
00455 
00456    int i;                                    /* file counter */
00457 
00458    for (i = 0;  i < option_data->a;  i++)
00459      {
00460        /*----- temporary file name -----*/
00461        sprintf (filename, "y.%d", i);
00462 
00463        /*-----  see if file already exists -----*/
00464        check_one_temporary_file (filename);
00465      }
00466 
00467 
00468    check_one_temporary_file ("ysum");
00469    check_one_temporary_file ("ss");
00470    check_one_temporary_file ("ssto");
00471    check_one_temporary_file ("sstr");
00472    check_one_temporary_file ("sse");
00473  
00474 }
00475 
00476 
00477 /*---------------------------------------------------------------------------*/
00478 /*
00479    Routine to check whether output files already exist.
00480 */
00481 
00482 void check_output_files (anova_options * option_data)
00483 {
00484   int i;       /* index */
00485 
00486   if (option_data->nftr > 0)   
00487     check_one_output_file (option_data, option_data->ftrname);
00488 
00489   if (option_data->num_ameans > 0)
00490     for (i = 0;  i < option_data->num_ameans;  i++)
00491       check_one_output_file (option_data, option_data->amname[i]);
00492 
00493   if (option_data->num_adiffs > 0)
00494     for (i = 0;  i < option_data->num_adiffs;  i++)
00495       check_one_output_file (option_data, option_data->adname[i]);
00496 
00497   if (option_data->num_acontr > 0)
00498     for (i = 0;  i < option_data->num_acontr;  i++)
00499       check_one_output_file (option_data, option_data->acname[i]);
00500 
00501   if (option_data->bucket_filename != NULL)
00502     check_one_output_file (option_data, option_data->bucket_filename);
00503 }
00504 
00505 
00506 /*---------------------------------------------------------------------------*/
00507 /*
00508   Routine to check for valid inputs.
00509 */
00510 
00511 void check_for_valid_inputs (anova_options * option_data)
00512 {
00513   int i;                               /* factor level index */
00514   char message[MAX_NAME_LENGTH];       /* error message */
00515   
00516   if (option_data->a < 2)
00517     ANOVA_error ("must specify number of factor levels (r>1) ");
00518   
00519   if (option_data->nt < option_data->a + 1) 
00520     ANOVA_error ("too few data sets for ANOVA");
00521 
00522   for (i = 0;  i < option_data->a;  i++)
00523     if (option_data->na[i] == 0)
00524       {
00525         sprintf (message, "level %d has too few data sets", i+1);
00526         ANOVA_error (message);
00527       }
00528 }
00529 /*---------------------------------------------------------------------------*/
00530 /*
00531   Routine to calculate the number of data files that have to be stored.
00532   Modified to account for 'bucket' dataset output.
00533 */
00534 
00535 int required_data_files (anova_options * option_data)
00536 {
00537   int r;                           /* number of factor levels */
00538   int nftr;                        /* number of F-treatment output files 
00539                                       (0 or 1) */
00540   int nmeans;                      /* number of estimated mean output files */
00541   int ndiffs;                      /* number of difference output files */
00542   int ncontr;                      /* number of contrast output files */
00543   int nout;                        /* number of output files */
00544   int nmax;                        /* maximum number of disk files */
00545 
00546 
00547   /*----- initialize local variables -----*/
00548   r = option_data->a;
00549   nftr = option_data->nftr;
00550   nmeans = option_data->num_ameans;
00551   ndiffs = option_data->num_adiffs;
00552   ncontr = option_data->num_acontr;
00553   nout = nftr + nmeans + ndiffs + ncontr;
00554 
00555   nmax = r + 2 + nout;
00556   if (option_data->bucket_filename != NULL)
00557     nmax = max (nmax, 2*nout);
00558   
00559   return (nmax);
00560 }
00561 
00562 
00563 /*---------------------------------------------------------------------------*/
00564 /*
00565   Routine to perform all ANOVA initialization.
00566 */
00567 
00568 void initialize (int argc,  char ** argv,  anova_options ** option_data)
00569 {
00570   int i;                               /* factor level index */
00571   
00572   
00573 
00574   /*----- save command line for history notes -----*/
00575   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00576 
00577 
00578   /*----- allocate memory space for input data -----*/   
00579   *option_data = (anova_options *) malloc(sizeof(anova_options));
00580   if (*option_data == NULL)
00581     ANOVA_error ("memory allocation error");
00582   
00583   /*----- get command line inputs -----*/
00584   get_options(argc, argv, *option_data);
00585   
00586   /*----- use first data set to get data set dimensions -----*/
00587   (*option_data)->first_dataset = (*option_data)->xname[0][0][0][0];
00588   get_dimensions (*option_data);
00589   printf ("Data set dimensions:  nx = %d  ny = %d  nz = %d  nxyz = %d \n",
00590           (*option_data)->nx, (*option_data)->ny,
00591           (*option_data)->nz, (*option_data)->nxyz);
00592   if ((*option_data)->nvoxel > (*option_data)->nxyz)
00593     ANOVA_error ("argument of -voxel is too large");
00594   
00595   /*----- count number of observations per treatment level -----*/
00596   (*option_data)->nt = 0;
00597   for (i = 0;  i < (*option_data)->a;  i++)
00598     (*option_data)->nt += (*option_data)->na[i];
00599 
00600   /*----- check for valid inputs -----*/
00601   check_for_valid_inputs (*option_data);
00602     
00603   /*----- check whether temporary files already exist -----*/
00604   check_temporary_files (*option_data);
00605   
00606   /*----- check whether output files already exist -----*/
00607   check_output_files (*option_data);
00608   
00609   /*----- check whether there is sufficient disk space -----*/
00610   if ((*option_data)->diskspace)  check_disk_space (*option_data);
00611 }
00612 
00613 
00614 /*---------------------------------------------------------------------------*/
00615 /*
00616    Routine to sum over the observations within one treatment level.
00617    Output is stored (temporarily) in data file y"i".3danova, where
00618    "i" = treatment level. 
00619 */
00620 
00621 void calculate_y (anova_options * option_data)
00622 {
00623    float * x = NULL;                /* pointer to original input data */
00624    float * y = NULL;                /* pointer to sum within treatment data */
00625    int i;                           /* factor level index */
00626    int j;                           /* observation number index */
00627    int r;                           /* number of factor levels (treatments) */
00628    int ixyz, nxyz;                  /* voxel counters */
00629    int nvoxel;                      /* output voxel # */
00630    char filename[MAX_NAME_LENGTH];  /* name of file for sum within treatment */
00631 
00632 
00633    /*----- initialize local variables -----*/
00634    r = option_data->a;
00635    nxyz = option_data->nxyz;
00636    nvoxel = option_data->nvoxel;
00637  
00638    /*----- allocate memory space for calculations -----*/
00639    x = (float *) malloc(sizeof(float)*nxyz);
00640    y = (float *) malloc(sizeof(float)*nxyz);
00641    if ((x == NULL) || (y == NULL))
00642       ANOVA_error ("unable to allocate sufficient memory");
00643   
00644 
00645    /*----- loop over treatment levels -----*/
00646    for (i = 0;  i < r;  i++)
00647    {
00648      /*----- sum observations within a treatment level -----*/
00649       volume_zero (y, nxyz);
00650       for (j = 0;  j < option_data->na[i];  j++)
00651       {
00652          read_afni_data (option_data, 
00653                          option_data->xname[i][0][0][j], x);
00654          if (nvoxel > 0)
00655             printf ("y[%d][%d] = %f \n", i+1, j+1, x[nvoxel-1]);
00656          for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00657             y[ixyz] += x[ixyz];
00658       }
00659  
00660       /*----- save the sum for this treatment level -----*/
00661       if (nvoxel > 0)
00662          printf ("y[%d]. = %f \n", i+1, y[nvoxel-1]); 
00663       sprintf (filename, "y.%d", i);
00664       volume_write (filename, y, nxyz);
00665    }
00666 
00667    /*----- release memory -----*/
00668    free (y);   y = NULL;
00669    free (x);   x = NULL;
00670 
00671 }
00672 
00673 
00674 /*---------------------------------------------------------------------------*/
00675 /*
00676    Routine to  sum over all observations at all treatment levels.
00677    The sum is stored (temporarily) in disk file ysum.3danova.
00678 */
00679 
00680 void calculate_ysum (anova_options * option_data)
00681 {
00682    float * y = NULL;                 /* pointer to sum within treatment data */
00683    float * ysum = NULL;              /* pointer to overall sum data */
00684    int i;                            /* factor level index */
00685    int ixyz, nxyz;                   /* voxel counters */
00686    int nvoxel;                       /* output voxel # */
00687    char filename[MAX_NAME_LENGTH];   /* sum within treatment input file name */
00688 
00689 
00690    /*----- initialize local variables -----*/
00691    nxyz = option_data->nxyz;
00692    nvoxel = option_data->nvoxel;
00693   
00694    /*----- allocate memory space for calculations -----*/
00695    y = (float *) malloc(sizeof(float)*nxyz);
00696    ysum = (float *) malloc(sizeof(float)*nxyz);
00697    if ((y == NULL) || (ysum == NULL))
00698       ANOVA_error ("unable to allocate sufficient memory");
00699   
00700 
00701    /*----- sum over all observations -----*/
00702    volume_zero (ysum, nxyz);
00703    for (i = 0;  i < option_data->a;  i++)
00704    { 
00705       sprintf (filename, "y.%d", i);
00706       volume_read (filename, y, nxyz);
00707       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00708          ysum[ixyz] += y[ixyz];
00709    }
00710 
00711    if (nvoxel > 0)
00712       printf ("y.. = %f \n", ysum[nvoxel-1]);
00713    volume_write ("ysum",ysum, nxyz);
00714 
00715    /*----- release memory -----*/
00716    free (ysum);  ysum = NULL;
00717    free (y);     y = NULL;
00718 
00719 }
00720 
00721 
00722 /*---------------------------------------------------------------------------*/
00723 /*
00724    Routine to calculate the sum of squares of all observations (SS).
00725    The output is stored (temporarily) in disk file ss.3danova.
00726 */
00727 
00728 void calculate_ss (anova_options * option_data)
00729 {
00730    float * x = NULL;                   /* pointer to original input data */
00731    float * ss = NULL;                  /* pointer to sum of squares data */
00732    int i;                              /* factor level index */
00733    int j;                              /* observation data index */
00734    int r;                              /* number of factor levels */
00735    int ixyz, nxyz;                     /* voxel counters */
00736    int nvoxel;                         /* output voxel # */
00737    float xval;                         /* float input data */
00738 
00739  
00740    /*----- initialize local variables -----*/
00741    r = option_data->a;
00742    nxyz = option_data->nxyz;
00743    nvoxel = option_data->nvoxel;
00744   
00745    /*----- allocate memory space for calculations -----*/
00746    x = (float *) malloc(sizeof(float)*nxyz);
00747    ss = (float *) malloc(sizeof(float)*nxyz);
00748    if ((x == NULL) || (ss == NULL))
00749       ANOVA_error ("unable to allocate sufficient memory");
00750 
00751    /*----- sum the squares of all observations -----*/
00752    volume_zero (ss, nxyz);
00753    for (i = 0;  i < r;  i++)
00754       for (j = 0;  j < option_data->na[i];  j++)
00755       {  
00756          read_afni_data (option_data, 
00757                          option_data->xname[i][0][0][j], x);
00758          for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00759          {
00760             xval = x[ixyz];
00761             ss[ixyz] += xval * xval;
00762          }
00763       }
00764 
00765    if (nvoxel > 0)
00766       printf ("SS = %f \n", ss[nvoxel-1]);
00767    volume_write ("ss", ss, nxyz);
00768 
00769    /*----- release memory -----*/
00770    free (ss);  ss = NULL;
00771    free (x);   x = NULL;
00772 
00773 }
00774 
00775 /*---------------------------------------------------------------------------*/
00776 /*
00777    Routine to calculate total (corrected for the mean) sum of squares (SSTO).
00778    The output is stored (temporarily) in disk file ssto.3danova.
00779 */
00780 
00781 void calculate_ssto (anova_options * option_data)
00782 {
00783    float * ysum = NULL;                /* ptr to sum over all observations */
00784    float * ssto = NULL;                /* pointer to ssto data */
00785    int ixyz, nxyz;                     /* voxel counters */
00786    int nvoxel;                         /* output voxel # */
00787    int nt;                             /* total number of observations */
00788    float yval;                         /* sum over all observations */
00789 
00790 
00791    /*----- initialize local variables -----*/
00792    nt = option_data->nt;
00793    nxyz = option_data->nxyz;
00794    nvoxel = option_data->nvoxel;
00795  
00796    /*----- allocate memory space for calculations -----*/
00797    ysum = (float *) malloc(sizeof(float)*nxyz);
00798    ssto = (float *) malloc(sizeof(float)*nxyz);
00799    if ((ysum == NULL) || (ssto == NULL))
00800       ANOVA_error ("unable to allocate sufficient memory");
00801  
00802    /*-----  SSTO = SS - ((YSUM)^2)/nt  -----*/
00803    volume_read ("ss", ssto, nxyz);
00804    volume_read ("ysum", ysum, nxyz); 
00805    for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00806    {
00807       yval = ysum[ixyz];
00808       ssto[ixyz] -= yval * yval / nt;
00809    }
00810 
00811    /*----- ss data file is no longer needed -----*/
00812    volume_delete ("ss");
00813 
00814    /*----- save total (corrected) sum of squares -----*/ 
00815    if (nvoxel > 0)
00816       printf ("SSTO = %f \n", ssto[nvoxel-1]);
00817    volume_write ("ssto", ssto, nxyz);
00818 
00819    /*----- release memory -----*/
00820    free (ssto);   ssto = NULL;
00821    free (ysum);   ysum = NULL;
00822 
00823 }
00824 
00825 /*---------------------------------------------------------------------------*/
00826 /*
00827    Routine to calculate the sum of squares due to the treatment (SSTR).
00828    The output is stored (temporarily) in file sstr.3danova.
00829 */
00830 
00831 void calculate_sstr (anova_options * option_data)
00832 {
00833    float * y = NULL;                   /* input data pointer */
00834    float * sstr = NULL;                /* sstr data pointer */
00835    int i;                              /* factor level index */
00836    int ixyz, nxyz;                     /* voxel counters */
00837    int nvoxel;                         /* output voxel # */
00838    int ni;                             /* number of observations at level i */
00839    int nt;                             /* total number of observations */
00840    float yval;                         /* temporary float value */
00841    char filename[MAX_NAME_LENGTH];     /* input data file name */
00842 
00843 
00844    /*----- assign local variables -----*/
00845    nt = option_data->nt;
00846    nxyz = option_data->nxyz;
00847    nvoxel = option_data->nvoxel;
00848  
00849    /*----- allocate memory space for calculations -----*/
00850    y = (float *) malloc(sizeof(float)*nxyz);
00851    sstr = (float *) malloc(sizeof(float)*nxyz);
00852    if ((y == NULL) || (sstr == NULL))
00853       ANOVA_error ("unable to allocate sufficient memory");
00854  
00855    volume_zero (sstr, nxyz);
00856    for (i = 0;  i < option_data->a;  i++)
00857    {
00858       ni = option_data->na[i];
00859       sprintf (filename, "y.%d", i);
00860       volume_read (filename, y, nxyz);
00861       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00862       {
00863          yval = y[ixyz];
00864          sstr[ixyz] += yval * yval / ni; 
00865       }
00866    }
00867 
00868    volume_read ("ysum", y, nxyz);
00869    for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00870    {
00871       yval = y[ixyz];
00872       sstr[ixyz] -= yval * yval / nt;
00873 
00874       /*----- protection against round-off error -----*/
00875       if (sstr[ixyz] < 0.0)  sstr[ixyz] = 0.0;
00876    }
00877 
00878    /*----- ysum data file is no longer needed -----*/
00879    volume_delete ("ysum");
00880 
00881    /*----- save treatment sum of squares -----*/
00882    if (nvoxel > 0)
00883       printf ("SSTR = %f \n", sstr[nvoxel-1]);
00884    volume_write ("sstr", sstr, nxyz);
00885 
00886    /*----- release memory -----*/
00887    free (sstr);   sstr = NULL;
00888    free (y);      y = NULL;
00889 
00890 }
00891 
00892 
00893 /*---------------------------------------------------------------------------*/
00894 /*
00895    Routine to calculate the error sum of squares,  SSE = SSTO - SSTR .
00896    The result is stored (temporarily) in disk file sse.3danova.
00897 */
00898 
00899 void calculate_sse (anova_options * option_data)
00900 {
00901    float * sstr = NULL;                /* pointer to sstr data */
00902    float * sse = NULL;                 /* pointer to sse data */
00903    int ixyz, nxyz;                     /* voxel counters */
00904    int nvoxel;                         /* output voxel # */
00905 
00906 
00907    /*----- assign local variables -----*/
00908    nxyz = option_data->nxyz;
00909    nvoxel = option_data->nvoxel;
00910    
00911    /*----- allocate memory space for calculations -----*/
00912    sstr = (float *) malloc(sizeof(float)*nxyz);
00913    sse = (float *) malloc(sizeof(float)*nxyz);
00914    if ((sstr == NULL) || (sse == NULL))
00915       ANOVA_error ("unable to allocate sufficient memory");
00916    
00917    /*----- read SSTO (total sum of squares) -----*/
00918    volume_read ("ssto", sse, nxyz);
00919 
00920    /*----- read SSTR (treatment sum of squares) -----*/
00921    volume_read ("sstr", sstr, nxyz);
00922    for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00923    {
00924       /*-----  SSE = SSTO - SSTR  -----*/
00925       sse[ixyz] -= sstr[ixyz];
00926 
00927       /*----- protection against round-off error -----*/
00928       if (sse[ixyz] < 0.0)  sse[ixyz] = 0.0;
00929    }
00930 
00931    /*----- ssto data file is no longer needed -----*/
00932    volume_delete ("ssto");
00933 
00934    /*----- save error sum of squares -----*/
00935    if (nvoxel > 0)
00936       printf ("SSE = %f \n", sse[nvoxel-1]);
00937    volume_write ("sse", sse, nxyz);
00938 
00939    /*----- release memory -----*/
00940    free (sse);    sse  = NULL;
00941    free (sstr);   sstr = NULL;
00942 
00943 }
00944 
00945 
00946 /*---------------------------------------------------------------------------*/
00947 /*
00948    Routine to calculate the F-statistic for treatment, F = MSTR / MSE,
00949       where MSTR = SSTR / (r-1),
00950       and   MSE  = SSE  / (n-r).
00951  
00952    The output is stored as a 2 sub-brick AFNI data set.  The first sub-brick 
00953    contains the square root of MSTR (mean sum of squares due to treatment), 
00954    and the second sub-brick contains the corresponding F-statistic. 
00955 */
00956 
00957 void calculate_f_statistic (anova_options * option_data)
00958 {
00959    const float  EPSILON = 1.0e-10;      /* protect against divide by zero */
00960    float * fin = NULL;                  /* pointer to input float data */
00961    float * mstr = NULL;                 /* pointer to MSTR data */
00962    float * ftr = NULL;                  /* pointer to F due-to-treatment */
00963    int r;                               /* number of factor levels */
00964    int nt;                              /* total number of observations */
00965    int ixyz, nxyz;                      /* voxel counters */
00966    int nvoxel;                          /* output voxel # */
00967    float fval;                          /* float value used in calculations */
00968    float mse;                           /* mean square error */
00969  
00970 
00971    /*----- initialize local variables -----*/
00972    r = option_data->a;
00973    nt = option_data->nt;
00974    nxyz = option_data->nxyz;
00975    nvoxel = option_data->nvoxel;
00976    
00977    /*----- allocate memory space for calculations -----*/
00978    /*----- Note:  the order in which memory allocations take place  
00979                   is important! -----*/
00980    ftr = (float *) malloc(sizeof(float)*nxyz);
00981    mstr = (float *) malloc(sizeof(float)*nxyz);
00982    fin = (float *) malloc(sizeof(float)*nxyz);
00983    if ((fin == NULL) || (ftr == NULL) || (mstr == NULL))
00984       ANOVA_error ("unable to allocate sufficient memory");
00985 
00986 
00987    /*----- calculate mean SS due to treatments -----*/
00988    volume_read ("sstr", fin, nxyz);   
00989    fval = 1.0 / (r - 1.0);
00990    for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00991       mstr[ixyz] = fin[ixyz] * fval;     /*---  MSTR = SSTR / (r-1)  ---*/
00992    if (nvoxel > 0)
00993       printf ("MSTR = %f \n", mstr[nvoxel-1]);
00994 
00995    /*----- calculate F-statistic  FTR = MSTR / MSE  -----*/
00996    volume_read ("sse", fin, nxyz);     
00997    fval = 1.0 / (nt - r);
00998    for (ixyz = 0;  ixyz < nxyz;  ixyz++)
00999      {
01000        mse = fin[ixyz] * fval;            /*---  MSE = SSE / (nt-r)  ---*/
01001        if (mse < EPSILON)
01002          ftr[ixyz] = 0.0;
01003        else
01004          ftr[ixyz] = mstr[ixyz] / mse;      /*---  FTR = MSTR / MSE  ---*/
01005      }
01006    if (nvoxel > 0)
01007       printf ("FTR = %f \n", ftr[nvoxel-1]);
01008 
01009    /*----- release memory -----*/
01010    free (fin);   fin = NULL;
01011 
01012    /*----- write out afni data file -----*/
01013    for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01014       mstr[ixyz] = sqrt(mstr[ixyz]);      /*-- mstr now holds square root --*/
01015    write_afni_data (option_data, option_data->ftrname, mstr, ftr, r-1, nt-r);
01016 
01017    /*----- this data file is no longer needed -----*/
01018    volume_delete ("sstr");
01019 
01020    /*----- release memory -----*/
01021    free (mstr);   mstr = NULL;
01022    free (ftr);    ftr = NULL;
01023 
01024 }
01025 
01026 
01027 /*---------------------------------------------------------------------------*/
01028 /*
01029    Routine to calculate the mean treatment effect at the user specified
01030    treatment level.  The output is stored as a 2 sub-brick AFNI data set.
01031    The first sub-brick contains the estimated mean at this treatment level, 
01032    and the second sub-brick contains the corresponding t-statistic.
01033 */
01034 
01035 void calculate_means (anova_options * option_data)
01036 {
01037    const float  EPSILON = 1.0e-10;    /* protect against divide by zero */
01038    float * fin = NULL;                /* pointer to float input data */
01039    float * mean = NULL;               /* pointer to treatment mean data */
01040    float * tmean = NULL;              /* pointer to t-statistic data */
01041    int imean;                         /* output mean option index */
01042    int level;                         /* factor level index */
01043    int ni;                            /* number obs. at factor level i */
01044    int ixyz, nxyz;                    /* voxel counters */
01045    int nvoxel;                        /* output voxel # */
01046    int r;                             /* number of factor levels */
01047    int nt;                            /* total number of observations */
01048    int num_means;                     /* number of user requested means */
01049    float fval;                        /* for calculating std. dev. */
01050    float stddev;                      /* est. std. dev. of factor mean */
01051    char filename[MAX_NAME_LENGTH];    /* input data file name */
01052  
01053 
01054    /*----- initialize local variables -----*/
01055    r = option_data->a;
01056    nt = option_data->nt;
01057    num_means = option_data->num_ameans;
01058    nxyz = option_data->nxyz;
01059    nvoxel = option_data->nvoxel;
01060 
01061    /*----- allocate memory space for calculations -----*/
01062    mean = (float *) malloc(sizeof(float)*nxyz);
01063    tmean = (float *) malloc(sizeof(float)*nxyz);
01064    if ((mean == NULL) || (tmean == NULL))  
01065       ANOVA_error ("unable to allocate sufficient memory");
01066    
01067    /*----- loop over user specified treatment means -----*/ 
01068    for (imean = 0;  imean < num_means;  imean++)
01069    {
01070       level = option_data->ameans[imean];
01071       ni = option_data->na[level];
01072 
01073       /*----- allocate temporary memory space -----*/
01074       fin = (float *) malloc(sizeof(float)*nxyz);
01075       if (fin == NULL)  ANOVA_error ("unable to allocate sufficient memory");
01076  
01077       /*----- estimate factor mean for this treatment level -----*/
01078       sprintf (filename, "y.%d", level); 
01079       volume_read (filename, fin, nxyz);
01080       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01081          mean[ixyz] =  fin[ixyz] / ni;
01082       if (nvoxel > 0)
01083          printf ("Mean [%d] = %f \n", level+1, mean[nvoxel-1]);
01084 
01085       /*----- divide by estimated standard deviation of factor mean -----*/
01086       volume_read ("sse", fin, nxyz);      
01087       fval = (1.0 / (nt-r)) * (1.0 / ni);
01088       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01089       {
01090          stddev =  sqrt(fin[ixyz] * fval);
01091          if (stddev < EPSILON)
01092            tmean[ixyz] = 0.0;
01093          else
01094            tmean[ixyz] = mean[ixyz] / stddev;
01095       }
01096       if (nvoxel > 0)
01097          printf ("t-Mean [%d] = %f \n", level+1, tmean[nvoxel-1]);
01098  
01099       /*----- release memory -----*/
01100       free (fin);   fin = NULL;
01101 
01102       /*----- write out afni data file -----*/
01103       write_afni_data (option_data, option_data->amname[imean], 
01104                        mean, tmean, nt-r, 0);
01105 
01106    }
01107 
01108    /*----- release memory -----*/
01109    free (tmean);  tmean = NULL;
01110    free (mean);   mean = NULL;
01111 
01112 }
01113 
01114 
01115 /*---------------------------------------------------------------------------*/
01116 /*
01117    Routine to estimate the difference in the means between two user specified
01118    treatment levels.  The output is a 2 sub-brick AFNI data set.  The first
01119    sub-brick contains the estimated difference in the means.  The second
01120    sub-brick contains the corresponding t-statistic.
01121 */
01122 
01123 void calculate_differences (anova_options * option_data)
01124 {
01125    const float  EPSILON = 1.0e-10;     /* protect against divide by zero */
01126    float * fin = NULL;                 /* pointer to float input data */
01127    float * diff = NULL;                /* pointer to est. diff. in means */
01128    float * tdiff = NULL;               /* pointer to t-statistic data */
01129    int r;                              /* number of factor levels */
01130    int nt;                             /* total number of observations */
01131    int ixyz, nxyz;                     /* voxel counters */
01132    int nvoxel;                         /* output voxel # */
01133    int num_diffs;                      /* number of user requested diffs. */
01134    int idiff;                          /* index for requested differences */
01135    int i, j;                           /* factor level indices */
01136    int ni, nj;                         /* number of obs. at levels i and j */
01137    float fval;                         /* for calculating std. dev. */
01138    float stddev;                       /* est. std. dev. of difference */
01139    char filename[MAX_NAME_LENGTH];     /* input file name */
01140 
01141 
01142    /*----- initialize local variables -----*/
01143    r = option_data->a;
01144    nt = option_data->nt;
01145    num_diffs = option_data->num_adiffs;
01146    nxyz = option_data->nxyz;
01147    nvoxel = option_data->nvoxel;
01148    
01149    /*----- allocate memory space for calculations -----*/
01150    diff = (float *) malloc(sizeof(float)*nxyz);
01151    tdiff = (float *) malloc(sizeof(float)*nxyz);
01152    if ((diff == NULL) || (tdiff == NULL))
01153       ANOVA_error ("unable to allocate sufficient memory");
01154 
01155    /*----- loop over user specified treatment differences -----*/
01156    for (idiff = 0;  idiff < num_diffs;  idiff++)
01157    {
01158  
01159       /*----- allocate temporary memory space -----*/
01160       fin = (float *) malloc(sizeof(float)*nxyz);
01161       if (fin == NULL)  ANOVA_error ("unable to allocate sufficient memory");
01162  
01163      /*----- read first treatment level mean -----*/
01164       i = option_data->adiffs[idiff][0];
01165       ni = option_data->na[i];
01166       sprintf (filename, "y.%d", i);
01167       volume_read (filename, fin, nxyz); 
01168       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01169          diff[ixyz] = fin[ixyz] / ni;
01170 
01171       /*----- subtract second treatment level mean -----*/
01172       j = option_data->adiffs[idiff][1];
01173       nj = option_data->na[j];
01174       sprintf (filename, "y.%d", j);
01175       volume_read (filename, fin, nxyz);
01176       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01177          diff[ixyz] -= fin[ixyz] / nj;
01178       if (nvoxel > 0)
01179          printf ("Diff [%d] - [%d] = %f \n", i+1, j+1, diff[nvoxel-1]);
01180 
01181       /*----- divide by estimated standard deviation of difference -----*/
01182       volume_read ("sse", fin, nxyz);
01183       fval = (1.0 / (nt-r)) * ((1.0/ni) + (1.0/nj));
01184       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01185       {
01186          stddev = sqrt (fin[ixyz] * fval);
01187          if (stddev < EPSILON)
01188            tdiff[ixyz] = 0.0;
01189          else
01190            tdiff[ixyz] = diff[ixyz] / stddev;
01191       }          
01192       if (nvoxel > 0)
01193          printf ("t-Diff [%d] - [%d] = %f \n", i+1, j+1, tdiff[nvoxel-1]);
01194 
01195       /*----- release memory -----*/
01196       free (fin);   fin = NULL;
01197 
01198       /*----- write out afni data file -----*/
01199       write_afni_data (option_data, option_data->adname[idiff], 
01200                        diff, tdiff, nt-r, 0);
01201 
01202    }
01203 
01204    /*----- release memory -----*/
01205    free (tdiff);  tdiff = NULL;
01206    free (diff);   diff = NULL;
01207 
01208 }
01209 
01210 
01211 /*---------------------------------------------------------------------------*/
01212 /*
01213    Routine to estimate a user specified contrast in treatment levels.
01214    The output is stored as a 2 sub-brick AFNI data set.  The first
01215    sub-brick contains the estimated contrast.  The second sub-brick contains 
01216    the corresponding t-statistic.
01217 */
01218 
01219 void calculate_contrasts (anova_options * option_data)
01220 {
01221    const float  EPSILON = 1.0e-10;     /* protect against divide by zero */
01222    float * fin = NULL;                 /* pointer to float input data */
01223    float * contr = NULL;               /* pointer to contrast estimate */
01224    float * tcontr = NULL;              /* pointer to t-statistic data */
01225    int r;                              /* number of factor levels */
01226    int nt;                             /* total number of observations */
01227    int ixyz, nxyz;                     /* voxel counters */
01228    int nvoxel;                         /* output voxel # */
01229    int num_contr;                      /* number of user requested contrasts */
01230    int icontr;                         /* index of user requested contrast */
01231    int level;                          /* factor level index */
01232    int ni;                             /* number of obs. at factor level i */
01233    float fval;                         /* for calculating std. dev. */
01234    float c;                            /* contrast coefficient */
01235    float stddev;                       /* est. std. dev. of contrast */
01236    char filename[MAX_NAME_LENGTH];     /* input data file name */
01237 
01238 
01239    /*----- initialize local variables -----*/
01240    r = option_data->a;
01241    nt = option_data->nt;
01242    num_contr = option_data->num_acontr;
01243    nxyz = option_data->nxyz;
01244    nvoxel = option_data->nvoxel;
01245    
01246    /*----- allocate memory space for calculations -----*/
01247    contr  = (float *) malloc(sizeof(float)*nxyz);
01248    tcontr = (float *) malloc(sizeof(float)*nxyz);
01249    if ((contr == NULL) || (tcontr == NULL))
01250       ANOVA_error ("unable to allocate sufficient memory");
01251 
01252 
01253    /*----- loop over user specified constrasts -----*/
01254    for (icontr = 0;  icontr < num_contr;  icontr++)
01255    {
01256       volume_zero (contr, nxyz);
01257       fval = 0.0;
01258 
01259       /*----- allocate temporary memory space -----*/
01260       fin = (float *) malloc(sizeof(float)*nxyz);
01261       if (fin == NULL)  ANOVA_error ("unable to allocate sufficient memory");
01262   
01263       for (level = 0;  level < r;  level++)
01264       {
01265          c = option_data->acontr[icontr][level]; 
01266          if (c == 0.0) continue; 
01267     
01268          /*----- add c * treatment level mean to contrast -----*/
01269          sprintf (filename, "y.%d", level);
01270          volume_read (filename, fin, nxyz);
01271          ni = option_data->na[level];
01272          fval += c * c / ni;
01273          for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01274             contr[ixyz] += c * fin[ixyz] / ni;
01275       }
01276       if (nvoxel > 0)
01277          printf ("Contr no.%d = %f \n", icontr+1, contr[nvoxel-1]);
01278 
01279       /*----- divide by estimated standard deviation of the contrast -----*/
01280       volume_read ("sse", fin, nxyz);
01281       for (ixyz = 0;  ixyz < nxyz;  ixyz++)
01282       {
01283          stddev = sqrt ((fin[ixyz] / (nt-r)) * fval);
01284          if (stddev < EPSILON)
01285            tcontr[ixyz] = 0.0;
01286          else
01287            tcontr[ixyz] = contr[ixyz] / stddev;
01288       }   
01289       if (nvoxel > 0)
01290          printf ("t-Contr no.%d = %f \n", icontr+1, tcontr[nvoxel-1]);
01291 
01292       /*----- release memory -----*/
01293       free (fin);   fin = NULL;
01294 
01295       /*----- write out afni data file -----*/
01296       write_afni_data (option_data, option_data->acname[icontr], 
01297                        contr, tcontr, nt-r, 0);
01298    }
01299 
01300    /*----- release memory -----*/
01301    free (tcontr);   tcontr = NULL;
01302    free (contr);    contr = NULL;
01303 
01304 }
01305 
01306 
01307 /*---------------------------------------------------------------------------*/
01308 /*
01309    Routine to perform single factor ANOVA.
01310 */
01311 
01312 void calculate_anova (anova_options * option_data)
01313 {
01314 
01315    /*-----  sum observations within treatment -----*/
01316    calculate_y (option_data);
01317  
01318    /*-----  sum all observations  -----*/
01319    calculate_ysum (option_data);
01320 
01321    /*-----  calculate sum of squares  -----*/
01322    calculate_ss (option_data);
01323 
01324    /*-----  calculate total (corrected for the mean) sum of squares -----*/
01325    calculate_ssto (option_data);
01326 
01327    /*-----  calculate treatment sum of squares  -----*/
01328    calculate_sstr (option_data);
01329 
01330    /*-----  calculate error sum of squares  -----*/
01331    calculate_sse (option_data);
01332 
01333 }
01334 
01335 
01336 /*---------------------------------------------------------------------------*/
01337 /*
01338    Routine to analyze the results from a single factor ANOVA.
01339 */
01340 
01341 void analyze_results (anova_options * option_data)
01342 {
01343  
01344    /*-----  calculate F-statistic for treatment effect  -----*/
01345    if (option_data->nftr > 0)  calculate_f_statistic (option_data);
01346 
01347    /*-----  estimate factor level means  -----*/
01348    if (option_data->num_ameans > 0)  calculate_means (option_data);
01349 
01350    /*-----  estimate differences in factor level means -----*/
01351    if (option_data->num_adiffs > 0)  calculate_differences (option_data);
01352 
01353    /*-----  calculate contrasts in factor levels -----*/
01354    if (option_data->num_acontr > 0)  calculate_contrasts (option_data);
01355 
01356 }
01357 
01358 
01359 /*---------------------------------------------------------------------------*/
01360 /*
01361    Routine to create an AFNI "bucket" output dataset.
01362 */
01363 
01364 void create_bucket (anova_options * option_data)
01365 
01366 {
01367   char bucket_str[10000];             /* command line for program 3dbucket */
01368   char refit_str[10000];              /* command line for program 3drefit */
01369   THD_3dim_dataset * dset=NULL;       /* input afni data set pointer */
01370   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
01371   int i;                              /* file index */
01372   int ibrick;                         /* sub-brick index number */
01373   char str[100];                      /* temporary character string */
01374 
01375 
01376   /*----- read first dataset -----*/
01377   dset = THD_open_dataset (option_data->first_dataset) ;
01378   if( ! ISVALID_3DIM_DATASET(dset) ){
01379     fprintf(stderr,"*** Unable to open dataset file %s\n",
01380             option_data->first_dataset);
01381     exit(1) ;
01382   }
01383 
01384        if( DSET_IS_1D(dset) ) USE_1D_filenames(1) ; /* 14 Mar 2003 */
01385   else if( DSET_IS_3D(dset) ) USE_1D_filenames(3) ; /* 21 Mar 2003 */
01386   
01387 
01388   /*----- make an empty copy of this dataset -----*/
01389   new_dset = EDIT_empty_copy( dset ) ;
01390   THD_delete_3dim_dataset (dset , False);   dset = NULL;
01391   EDIT_dset_items (new_dset, 
01392                    ADN_directory_name, option_data->session,
01393                    ADN_none);
01394 
01395 
01396   /*----- begin command line for program 3dbucket -----*/
01397   strcpy (bucket_str, "3dbucket ");
01398   strcat (bucket_str, " -prefix ");
01399   strcat (bucket_str, option_data->bucket_filename);
01400 
01401 
01402   /*----- begin command line for program 3drefit -----*/
01403   strcpy (refit_str, "3drefit ");
01404   ibrick = -1;
01405 
01406  
01407   /*----- make F-statistic sub-bricks -----*/
01408   if (option_data->nftr != 0)
01409     {
01410       add_file_name (new_dset, option_data->ftrname, bucket_str);
01411 
01412       ibrick++;
01413       sprintf (str, " -sublabel %d %s:Inten ",
01414                ibrick, option_data->ftrname);
01415       strcat (refit_str, str);
01416 
01417       ibrick++;
01418       sprintf (str, " -sublabel %d %s:F-stat ", 
01419                ibrick, option_data->ftrname);
01420       strcat (refit_str, str);
01421     }
01422       
01423   
01424   /*----- make factor level mean sub-bricks -----*/
01425   if (option_data->num_ameans > 0)
01426     for (i = 0; i < option_data->num_ameans; i++)
01427       {
01428         add_file_name (new_dset, option_data->amname[i], bucket_str);
01429 
01430         ibrick++;
01431         sprintf (str, " -sublabel %d %s:Mean ", 
01432                  ibrick, option_data->amname[i]);
01433         strcat (refit_str, str);
01434 
01435         ibrick++;
01436         sprintf (str, " -sublabel %d %s:t-stat ", 
01437                  ibrick, option_data->amname[i]);
01438         strcat (refit_str, str);
01439       }
01440   
01441 
01442   /*----- make difference in factor level means sub-bricks -----*/
01443   if (option_data->num_adiffs > 0)
01444     for (i = 0; i < option_data->num_adiffs; i++)
01445       {
01446         add_file_name (new_dset, option_data->adname[i], bucket_str);
01447 
01448         ibrick++;
01449         sprintf (str, " -sublabel %d %s:Diff ", 
01450                  ibrick, option_data->adname[i]);
01451         strcat (refit_str, str);
01452 
01453         ibrick++;
01454         sprintf (str, " -sublabel %d %s:t-stat ", 
01455                  ibrick, option_data->adname[i]);
01456         strcat (refit_str, str);
01457       }
01458   
01459 
01460   /*----- make contrast in factor level means sub-bricks -----*/
01461   if (option_data->num_acontr > 0)
01462     for (i = 0; i < option_data->num_acontr; i++)
01463       {
01464         add_file_name (new_dset, option_data->acname[i], bucket_str);
01465 
01466         ibrick++;
01467         sprintf (str, " -sublabel %d %s:Contr ", 
01468                  ibrick, option_data->acname[i]);
01469         strcat (refit_str, str);
01470 
01471         ibrick++;
01472         sprintf (str, " -sublabel %d %s:t-stat ", 
01473                  ibrick, option_data->acname[i]);
01474         strcat (refit_str, str);
01475       }
01476 
01477 
01478   /*----- invoke program 3dbucket to generate bucket type output dataset
01479           by concatenating previous output files -----*/
01480   printf("Writing `bucket' dataset ");
01481   printf("into %s\n", option_data->bucket_filename);
01482   fprintf(stderr,"RUNNING COMMAND: %s\n",bucket_str);
01483   system (bucket_str);
01484 
01485 
01486   /*----- invoke program 3drefit to label individual sub-bricks -----*/
01487   add_file_name (new_dset, option_data->bucket_filename, refit_str);
01488   fprintf(stderr,"RUNNING COMMAND: %s\n",refit_str);
01489   system (refit_str);
01490 
01491 
01492   /*----- release memory -----*/
01493   THD_delete_3dim_dataset (new_dset , False);   new_dset = NULL;
01494   
01495 }
01496 
01497 
01498 /*---------------------------------------------------------------------------*/
01499 /*
01500    Routine to release memory and remove any remaining temporary data files.
01501    If the '-bucket' option has been used, create the bucket dataset and
01502    dispose of the 2 sub-brick datasets.
01503 */
01504 
01505 void terminate (anova_options * option_data)
01506 {
01507   int i;
01508   THD_3dim_dataset * dset=NULL;       /* input afni data set pointer */
01509   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
01510   char filename[MAX_NAME_LENGTH];
01511 
01512 
01513   /*----- remove temporary data files -----*/
01514   volume_delete ("sstr");
01515   volume_delete ("sse");
01516   for (i = 0; i < option_data->a; i++)
01517     {
01518       sprintf (filename, "y.%d", i);
01519       volume_delete (filename);
01520     }
01521   
01522   
01523   /*----- create bucket dataset -----*/
01524   if (option_data->bucket_filename != NULL)
01525     create_bucket (option_data);
01526 
01527   
01528   /*----- if 'bucket' datset was created, remove the individual 2-subbrick
01529           data files -----*/
01530   if (option_data->bucket_filename != NULL)
01531     {
01532 
01533       /*----- read first dataset -----*/
01534       dset = THD_open_dataset (option_data->first_dataset) ;
01535       if( ! ISVALID_3DIM_DATASET(dset) ){
01536         fprintf(stderr,"*** Unable to open dataset file %s\n",
01537                 option_data->first_dataset);
01538         exit(1) ;
01539       }
01540       
01541       /*----- make an empty copy of this dataset -----*/
01542       new_dset = EDIT_empty_copy (dset);
01543       THD_delete_3dim_dataset (dset , False);   dset = NULL;
01544       EDIT_dset_items (new_dset, 
01545                        ADN_directory_name, option_data->session,
01546                        ADN_none);
01547       
01548       /*----- remove F-statistic data file -----*/
01549       if (option_data->nftr != 0)
01550         remove_dataset (new_dset, option_data->ftrname);
01551       
01552       /*----- remove mean factor level data files -----*/
01553       if (option_data->num_ameans > 0)
01554         for (i = 0; i < option_data->num_ameans; i++)
01555           remove_dataset (new_dset, option_data->amname[i]);
01556       
01557       /*----- remove difference in factor levels data files -----*/
01558       if (option_data->num_adiffs > 0)
01559         for (i = 0; i < option_data->num_adiffs; i++)
01560           remove_dataset (new_dset, option_data->adname[i]);
01561       
01562       /*----- remove contrast in factor levels data files -----*/
01563       if (option_data->num_acontr > 0)
01564         for (i = 0; i < option_data->num_acontr; i++)
01565           remove_dataset (new_dset, option_data->acname[i]);
01566       
01567       THD_delete_3dim_dataset (new_dset , False);   new_dset = NULL;
01568     }  
01569 
01570 
01571   /*----- deallocate memory -----*/
01572   destroy_anova_options (option_data);
01573 
01574 }
01575 
01576 
01577 /*---------------------------------------------------------------------------*/
01578 /*
01579    Single factor analysis of variance (ANOVA).
01580 */
01581  
01582 int main (int argc, char ** argv)
01583 {
01584    anova_options * option_data = NULL;
01585 
01586 #if 0
01587    /*----- Identify software -----*/
01588    printf ("\n\n");
01589    printf ("Program:          %s \n", PROGRAM_NAME);
01590    printf ("Author:           %s \n", PROGRAM_AUTHOR); 
01591    printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
01592    printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
01593    printf ("\n");
01594 #endif
01595 
01596 
01597    /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
01598 
01599    mainENTRY("3dANOVA main"); machdep(); PRINT_VERSION("3dANOVA");
01600    { int new_argc ; char ** new_argv ;
01601      addto_args( argc , argv , &new_argc , &new_argv ) ;
01602      if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
01603    }
01604 
01605    /*----- program initialization -----*/
01606    initialize (argc, argv, &option_data);
01607 
01608    /*----- calculate sums of squares -----*/
01609    calculate_anova (option_data);
01610 
01611    /*----- generate requested output -----*/
01612    analyze_results (option_data);
01613 
01614    /*----- terminate program -----*/
01615    terminate (option_data);
01616    free (option_data);   option_data = NULL;
01617 
01618    exit(0);
01619 }