00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062 #define PROGRAM_NAME "3dANOVA"
00063 #define PROGRAM_AUTHOR "B. Douglas Ward"
00064 #define PROGRAM_INITIAL "09 Dec 1996"
00065 #define PROGRAM_LATEST "21 Jul 2005"
00066
00067
00068
00069 #define SUFFIX ".3danova"
00070
00071 #include "3dANOVA.h"
00072 #include "3dANOVA.lib"
00073
00074
00075
00076
00077
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
00178
00179
00180 void get_options (int argc, char ** argv, anova_options * option_data)
00181 {
00182 int nopt = 1;
00183 int ival;
00184 int i, j, k;
00185 int nijk;
00186 float fval;
00187 THD_3dim_dataset * dset=NULL;
00188 char message[MAX_NAME_LENGTH];
00189
00190
00191
00192 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00193
00194
00195 AFNI_logger (PROGRAM_NAME,argc,argv);
00196
00197
00198 initialize_options (option_data);
00199
00200
00201
00202 while (nopt < argc)
00203 {
00204
00205
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
00229 if( strncmp(argv[nopt],"-diskspace",5) == 0 )
00230 {
00231 option_data->diskspace = 1;
00232 nopt++ ; continue ;
00233 }
00234
00235
00236
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 ;
00251 }
00252
00253
00254
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
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
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
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
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
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
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
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
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
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
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
00441 sprintf (message,"Unrecognized command line option: %s\n", argv[nopt]);
00442 ANOVA_error (message);
00443 }
00444
00445 }
00446
00447
00448
00449
00450
00451
00452 void check_temporary_files (anova_options * option_data)
00453 {
00454 char filename[MAX_NAME_LENGTH];
00455
00456 int i;
00457
00458 for (i = 0; i < option_data->a; i++)
00459 {
00460
00461 sprintf (filename, "y.%d", i);
00462
00463
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
00480
00481
00482 void check_output_files (anova_options * option_data)
00483 {
00484 int i;
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
00509
00510
00511 void check_for_valid_inputs (anova_options * option_data)
00512 {
00513 int i;
00514 char message[MAX_NAME_LENGTH];
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
00532
00533
00534
00535 int required_data_files (anova_options * option_data)
00536 {
00537 int r;
00538 int nftr;
00539
00540 int nmeans;
00541 int ndiffs;
00542 int ncontr;
00543 int nout;
00544 int nmax;
00545
00546
00547
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
00566
00567
00568 void initialize (int argc, char ** argv, anova_options ** option_data)
00569 {
00570 int i;
00571
00572
00573
00574
00575 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00576
00577
00578
00579 *option_data = (anova_options *) malloc(sizeof(anova_options));
00580 if (*option_data == NULL)
00581 ANOVA_error ("memory allocation error");
00582
00583
00584 get_options(argc, argv, *option_data);
00585
00586
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
00596 (*option_data)->nt = 0;
00597 for (i = 0; i < (*option_data)->a; i++)
00598 (*option_data)->nt += (*option_data)->na[i];
00599
00600
00601 check_for_valid_inputs (*option_data);
00602
00603
00604 check_temporary_files (*option_data);
00605
00606
00607 check_output_files (*option_data);
00608
00609
00610 if ((*option_data)->diskspace) check_disk_space (*option_data);
00611 }
00612
00613
00614
00615
00616
00617
00618
00619
00620
00621 void calculate_y (anova_options * option_data)
00622 {
00623 float * x = NULL;
00624 float * y = NULL;
00625 int i;
00626 int j;
00627 int r;
00628 int ixyz, nxyz;
00629 int nvoxel;
00630 char filename[MAX_NAME_LENGTH];
00631
00632
00633
00634 r = option_data->a;
00635 nxyz = option_data->nxyz;
00636 nvoxel = option_data->nvoxel;
00637
00638
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
00646 for (i = 0; i < r; i++)
00647 {
00648
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
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
00668 free (y); y = NULL;
00669 free (x); x = NULL;
00670
00671 }
00672
00673
00674
00675
00676
00677
00678
00679
00680 void calculate_ysum (anova_options * option_data)
00681 {
00682 float * y = NULL;
00683 float * ysum = NULL;
00684 int i;
00685 int ixyz, nxyz;
00686 int nvoxel;
00687 char filename[MAX_NAME_LENGTH];
00688
00689
00690
00691 nxyz = option_data->nxyz;
00692 nvoxel = option_data->nvoxel;
00693
00694
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
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
00716 free (ysum); ysum = NULL;
00717 free (y); y = NULL;
00718
00719 }
00720
00721
00722
00723
00724
00725
00726
00727
00728 void calculate_ss (anova_options * option_data)
00729 {
00730 float * x = NULL;
00731 float * ss = NULL;
00732 int i;
00733 int j;
00734 int r;
00735 int ixyz, nxyz;
00736 int nvoxel;
00737 float xval;
00738
00739
00740
00741 r = option_data->a;
00742 nxyz = option_data->nxyz;
00743 nvoxel = option_data->nvoxel;
00744
00745
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
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
00770 free (ss); ss = NULL;
00771 free (x); x = NULL;
00772
00773 }
00774
00775
00776
00777
00778
00779
00780
00781 void calculate_ssto (anova_options * option_data)
00782 {
00783 float * ysum = NULL;
00784 float * ssto = NULL;
00785 int ixyz, nxyz;
00786 int nvoxel;
00787 int nt;
00788 float yval;
00789
00790
00791
00792 nt = option_data->nt;
00793 nxyz = option_data->nxyz;
00794 nvoxel = option_data->nvoxel;
00795
00796
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
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
00812 volume_delete ("ss");
00813
00814
00815 if (nvoxel > 0)
00816 printf ("SSTO = %f \n", ssto[nvoxel-1]);
00817 volume_write ("ssto", ssto, nxyz);
00818
00819
00820 free (ssto); ssto = NULL;
00821 free (ysum); ysum = NULL;
00822
00823 }
00824
00825
00826
00827
00828
00829
00830
00831 void calculate_sstr (anova_options * option_data)
00832 {
00833 float * y = NULL;
00834 float * sstr = NULL;
00835 int i;
00836 int ixyz, nxyz;
00837 int nvoxel;
00838 int ni;
00839 int nt;
00840 float yval;
00841 char filename[MAX_NAME_LENGTH];
00842
00843
00844
00845 nt = option_data->nt;
00846 nxyz = option_data->nxyz;
00847 nvoxel = option_data->nvoxel;
00848
00849
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
00875 if (sstr[ixyz] < 0.0) sstr[ixyz] = 0.0;
00876 }
00877
00878
00879 volume_delete ("ysum");
00880
00881
00882 if (nvoxel > 0)
00883 printf ("SSTR = %f \n", sstr[nvoxel-1]);
00884 volume_write ("sstr", sstr, nxyz);
00885
00886
00887 free (sstr); sstr = NULL;
00888 free (y); y = NULL;
00889
00890 }
00891
00892
00893
00894
00895
00896
00897
00898
00899 void calculate_sse (anova_options * option_data)
00900 {
00901 float * sstr = NULL;
00902 float * sse = NULL;
00903 int ixyz, nxyz;
00904 int nvoxel;
00905
00906
00907
00908 nxyz = option_data->nxyz;
00909 nvoxel = option_data->nvoxel;
00910
00911
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
00918 volume_read ("ssto", sse, nxyz);
00919
00920
00921 volume_read ("sstr", sstr, nxyz);
00922 for (ixyz = 0; ixyz < nxyz; ixyz++)
00923 {
00924
00925 sse[ixyz] -= sstr[ixyz];
00926
00927
00928 if (sse[ixyz] < 0.0) sse[ixyz] = 0.0;
00929 }
00930
00931
00932 volume_delete ("ssto");
00933
00934
00935 if (nvoxel > 0)
00936 printf ("SSE = %f \n", sse[nvoxel-1]);
00937 volume_write ("sse", sse, nxyz);
00938
00939
00940 free (sse); sse = NULL;
00941 free (sstr); sstr = NULL;
00942
00943 }
00944
00945
00946
00947
00948
00949
00950
00951
00952
00953
00954
00955
00956
00957 void calculate_f_statistic (anova_options * option_data)
00958 {
00959 const float EPSILON = 1.0e-10;
00960 float * fin = NULL;
00961 float * mstr = NULL;
00962 float * ftr = NULL;
00963 int r;
00964 int nt;
00965 int ixyz, nxyz;
00966 int nvoxel;
00967 float fval;
00968 float mse;
00969
00970
00971
00972 r = option_data->a;
00973 nt = option_data->nt;
00974 nxyz = option_data->nxyz;
00975 nvoxel = option_data->nvoxel;
00976
00977
00978
00979
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
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;
00992 if (nvoxel > 0)
00993 printf ("MSTR = %f \n", mstr[nvoxel-1]);
00994
00995
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;
01001 if (mse < EPSILON)
01002 ftr[ixyz] = 0.0;
01003 else
01004 ftr[ixyz] = mstr[ixyz] / mse;
01005 }
01006 if (nvoxel > 0)
01007 printf ("FTR = %f \n", ftr[nvoxel-1]);
01008
01009
01010 free (fin); fin = NULL;
01011
01012
01013 for (ixyz = 0; ixyz < nxyz; ixyz++)
01014 mstr[ixyz] = sqrt(mstr[ixyz]);
01015 write_afni_data (option_data, option_data->ftrname, mstr, ftr, r-1, nt-r);
01016
01017
01018 volume_delete ("sstr");
01019
01020
01021 free (mstr); mstr = NULL;
01022 free (ftr); ftr = NULL;
01023
01024 }
01025
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035 void calculate_means (anova_options * option_data)
01036 {
01037 const float EPSILON = 1.0e-10;
01038 float * fin = NULL;
01039 float * mean = NULL;
01040 float * tmean = NULL;
01041 int imean;
01042 int level;
01043 int ni;
01044 int ixyz, nxyz;
01045 int nvoxel;
01046 int r;
01047 int nt;
01048 int num_means;
01049 float fval;
01050 float stddev;
01051 char filename[MAX_NAME_LENGTH];
01052
01053
01054
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
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
01068 for (imean = 0; imean < num_means; imean++)
01069 {
01070 level = option_data->ameans[imean];
01071 ni = option_data->na[level];
01072
01073
01074 fin = (float *) malloc(sizeof(float)*nxyz);
01075 if (fin == NULL) ANOVA_error ("unable to allocate sufficient memory");
01076
01077
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
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
01100 free (fin); fin = NULL;
01101
01102
01103 write_afni_data (option_data, option_data->amname[imean],
01104 mean, tmean, nt-r, 0);
01105
01106 }
01107
01108
01109 free (tmean); tmean = NULL;
01110 free (mean); mean = NULL;
01111
01112 }
01113
01114
01115
01116
01117
01118
01119
01120
01121
01122
01123 void calculate_differences (anova_options * option_data)
01124 {
01125 const float EPSILON = 1.0e-10;
01126 float * fin = NULL;
01127 float * diff = NULL;
01128 float * tdiff = NULL;
01129 int r;
01130 int nt;
01131 int ixyz, nxyz;
01132 int nvoxel;
01133 int num_diffs;
01134 int idiff;
01135 int i, j;
01136 int ni, nj;
01137 float fval;
01138 float stddev;
01139 char filename[MAX_NAME_LENGTH];
01140
01141
01142
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
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
01156 for (idiff = 0; idiff < num_diffs; idiff++)
01157 {
01158
01159
01160 fin = (float *) malloc(sizeof(float)*nxyz);
01161 if (fin == NULL) ANOVA_error ("unable to allocate sufficient memory");
01162
01163
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
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
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
01196 free (fin); fin = NULL;
01197
01198
01199 write_afni_data (option_data, option_data->adname[idiff],
01200 diff, tdiff, nt-r, 0);
01201
01202 }
01203
01204
01205 free (tdiff); tdiff = NULL;
01206 free (diff); diff = NULL;
01207
01208 }
01209
01210
01211
01212
01213
01214
01215
01216
01217
01218
01219 void calculate_contrasts (anova_options * option_data)
01220 {
01221 const float EPSILON = 1.0e-10;
01222 float * fin = NULL;
01223 float * contr = NULL;
01224 float * tcontr = NULL;
01225 int r;
01226 int nt;
01227 int ixyz, nxyz;
01228 int nvoxel;
01229 int num_contr;
01230 int icontr;
01231 int level;
01232 int ni;
01233 float fval;
01234 float c;
01235 float stddev;
01236 char filename[MAX_NAME_LENGTH];
01237
01238
01239
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
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
01254 for (icontr = 0; icontr < num_contr; icontr++)
01255 {
01256 volume_zero (contr, nxyz);
01257 fval = 0.0;
01258
01259
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
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
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
01293 free (fin); fin = NULL;
01294
01295
01296 write_afni_data (option_data, option_data->acname[icontr],
01297 contr, tcontr, nt-r, 0);
01298 }
01299
01300
01301 free (tcontr); tcontr = NULL;
01302 free (contr); contr = NULL;
01303
01304 }
01305
01306
01307
01308
01309
01310
01311
01312 void calculate_anova (anova_options * option_data)
01313 {
01314
01315
01316 calculate_y (option_data);
01317
01318
01319 calculate_ysum (option_data);
01320
01321
01322 calculate_ss (option_data);
01323
01324
01325 calculate_ssto (option_data);
01326
01327
01328 calculate_sstr (option_data);
01329
01330
01331 calculate_sse (option_data);
01332
01333 }
01334
01335
01336
01337
01338
01339
01340
01341 void analyze_results (anova_options * option_data)
01342 {
01343
01344
01345 if (option_data->nftr > 0) calculate_f_statistic (option_data);
01346
01347
01348 if (option_data->num_ameans > 0) calculate_means (option_data);
01349
01350
01351 if (option_data->num_adiffs > 0) calculate_differences (option_data);
01352
01353
01354 if (option_data->num_acontr > 0) calculate_contrasts (option_data);
01355
01356 }
01357
01358
01359
01360
01361
01362
01363
01364 void create_bucket (anova_options * option_data)
01365
01366 {
01367 char bucket_str[10000];
01368 char refit_str[10000];
01369 THD_3dim_dataset * dset=NULL;
01370 THD_3dim_dataset * new_dset=NULL;
01371 int i;
01372 int ibrick;
01373 char str[100];
01374
01375
01376
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) ;
01385 else if( DSET_IS_3D(dset) ) USE_1D_filenames(3) ;
01386
01387
01388
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
01397 strcpy (bucket_str, "3dbucket ");
01398 strcat (bucket_str, " -prefix ");
01399 strcat (bucket_str, option_data->bucket_filename);
01400
01401
01402
01403 strcpy (refit_str, "3drefit ");
01404 ibrick = -1;
01405
01406
01407
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
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
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
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
01479
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
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
01493 THD_delete_3dim_dataset (new_dset , False); new_dset = NULL;
01494
01495 }
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505 void terminate (anova_options * option_data)
01506 {
01507 int i;
01508 THD_3dim_dataset * dset=NULL;
01509 THD_3dim_dataset * new_dset=NULL;
01510 char filename[MAX_NAME_LENGTH];
01511
01512
01513
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
01524 if (option_data->bucket_filename != NULL)
01525 create_bucket (option_data);
01526
01527
01528
01529
01530 if (option_data->bucket_filename != NULL)
01531 {
01532
01533
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
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
01549 if (option_data->nftr != 0)
01550 remove_dataset (new_dset, option_data->ftrname);
01551
01552
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
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
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
01572 destroy_anova_options (option_data);
01573
01574 }
01575
01576
01577
01578
01579
01580
01581
01582 int main (int argc, char ** argv)
01583 {
01584 anova_options * option_data = NULL;
01585
01586 #if 0
01587
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
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
01606 initialize (argc, argv, &option_data);
01607
01608
01609 calculate_anova (option_data);
01610
01611
01612 analyze_results (option_data);
01613
01614
01615 terminate (option_data);
01616 free (option_data); option_data = NULL;
01617
01618 exit(0);
01619 }