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
00063
00064
00065
00066
00067 #define PROGRAM_NAME "3dANOVA2"
00068 #define PROGRAM_AUTHOR "B. Douglas Ward"
00069 #define PROGRAM_INITIAL "09 Dec 1996"
00070 #define PROGRAM_LATEST "02 Aug 2005"
00071
00072
00073
00074 #define SUFFIX ".3danova2"
00075
00076 #include "3dANOVA.h"
00077 #include "3dANOVA.lib"
00078
00079
00080
00081
00082
00083
00084
00085 void display_help_menu()
00086 {
00087 printf
00088 (
00089 "This program performs two-factor ANOVA on 3D data sets \n\n"
00090 "Usage: \n"
00091 "3dANOVA2 \n"
00092 "-type k type of ANOVA model to be used: \n"
00093 " k=1 fixed effects model (A and B fixed) \n"
00094 " k=2 random effects model (A and B random) \n"
00095 " k=3 mixed effects model (A fixed, B random) \n"
00096 " \n"
00097 "-alevels a a = number of levels of factor A \n"
00098 "-blevels b b = number of levels of factor B \n"
00099 "-dset 1 1 filename data set for level 1 of factor A \n"
00100 " and level 1 of factor B \n"
00101 " . . . . . . \n"
00102 " \n"
00103 "-dset i j filename data set for level i of factor A \n"
00104 " and level j of factor B \n"
00105 " . . . . . . \n"
00106 " \n"
00107 "-dset a b filename data set for level a of factor A \n"
00108 " and level b of factor B \n"
00109 " \n"
00110 "[-voxel num] screen output for voxel # num \n"
00111 "[-diskspace] print out disk space required for \n"
00112 " program execution \n"
00113 " \n"
00114 " \n"
00115 "The following commands generate individual AFNI 2 sub-brick datasets: \n"
00116 " (In each case, output is written to the file with the specified \n"
00117 " prefix file name.) \n"
00118 " \n"
00119 "[-ftr prefix] F-statistic for treatment effect \n"
00120 "[-fa prefix] F-statistic for factor A effect \n"
00121 "[-fb prefix] F-statistic for factor B effect \n"
00122 "[-fab prefix] F-statistic for interaction \n"
00123 "[-amean i prefix] estimate mean of factor A level i \n"
00124 "[-bmean j prefix] estimate mean of factor B level j \n"
00125 "[-xmean i j prefix] estimate mean of cell at level i of \n"
00126 " factor A, level j of factor B \n"
00127 "[-adiff i j prefix] difference between levels i and j of \n"
00128 " factor A \n"
00129 "[-bdiff i j prefix] difference between levels i and j of \n"
00130 " factor B \n"
00131 "[-xdiff i j k l prefix] difference between cell mean at A=i,B=j \n"
00132 " and cell mean at A=k,B=l \n"
00133 "[-acontr c1 ... ca prefix] contrast in factor A levels \n"
00134 "[-bcontr c1 ... cb prefix] contrast in factor B levels \n"
00135 "[-xcontr c11 ... c1b c21 ... c2b ... ca1 ... cab prefix] \n"
00136 " contrast in cell means \n"
00137 " \n"
00138 " \n"
00139 "The following command generates one AFNI 'bucket' type dataset: \n"
00140 " \n"
00141 "[-bucket prefix] create one AFNI 'bucket' dataset whose \n"
00142 " sub-bricks are obtained by concatenating \n"
00143 " the above output files; the output 'bucket'\n"
00144 " is written to file with prefix file name \n"
00145 "\n");
00146
00147 printf
00148 (
00149 "\n"
00150 "N.B.: For this program, the user must specify 1 and only 1 sub-brick \n"
00151 " with each -dset command. That is, if an input dataset contains \n"
00152 " more than 1 sub-brick, a sub-brick selector must be used, e.g.: \n"
00153 " -dset 2 4 'fred+orig[3]' \n"
00154 );
00155
00156 printf("\n" MASTER_SHORTHELP_STRING ) ;
00157
00158 exit(0);
00159 }
00160
00161
00162
00163
00164
00165
00166
00167 void get_options (int argc, char ** argv, anova_options * option_data)
00168 {
00169 int nopt = 1;
00170 int ival, jval;
00171 int i, j, k;
00172 int nij;
00173 float fval;
00174 THD_3dim_dataset * dset=NULL;
00175 char message[MAX_NAME_LENGTH];
00176 int n[MAX_LEVELS][MAX_LEVELS];
00177
00178
00179
00180 if (argc < 2 || strncmp(argv[1], "-help", 5) == 0) display_help_menu();
00181
00182
00183 AFNI_logger (PROGRAM_NAME,argc,argv);
00184
00185
00186 initialize_options (option_data);
00187
00188
00189 for (i = 0; i < MAX_LEVELS; i++)
00190 for (j = 0; j < MAX_LEVELS; j++)
00191 n[i][j] = 0;
00192
00193
00194
00195 while (nopt < argc)
00196 {
00197
00198
00199 if ((option_data->xname == NULL) && (option_data->a > 0) &&
00200 (option_data->b > 0))
00201 {
00202 option_data->xname =
00203 (char *****) malloc (sizeof(char ****) * option_data->a);
00204 for (i = 0; i < option_data->a; i++)
00205 {
00206 option_data->xname[i] =
00207 (char ****) malloc (sizeof(char ***) * option_data->b);
00208 for (j = 0; j < option_data->b; j++)
00209 {
00210 option_data->xname[i][j] =
00211 (char ***) malloc (sizeof(char **) * 1);
00212 for (k = 0; k < 1; k++)
00213 {
00214 option_data->xname[i][j][k] =
00215 (char **) malloc (sizeof(char *) * MAX_OBSERVATIONS);
00216 }
00217 }
00218 }
00219 }
00220
00221
00222
00223 if( strncmp(argv[nopt],"-diskspace",5) == 0 )
00224 {
00225 option_data->diskspace = 1;
00226 nopt++ ; continue ;
00227 }
00228
00229
00230
00231 if( strncmp(argv[nopt],"-datum",5) == 0 ){
00232 if( ++nopt >= argc ) ANOVA_error("need an argument after -datum!") ;
00233
00234 if( strcmp(argv[nopt],"short") == 0 ){
00235 option_data->datum = MRI_short ;
00236 } else if( strcmp(argv[nopt],"float") == 0 ){
00237 option_data->datum = MRI_float ;
00238 } else {
00239 char buf[256] ;
00240 sprintf(buf,"-datum of type '%s' is not supported in 3dANOVA2!",
00241 argv[nopt] ) ;
00242 ANOVA_error(buf) ;
00243 }
00244 nopt++ ; continue ;
00245 }
00246
00247
00248
00249 if( strncmp(argv[nopt],"-session",5) == 0 ){
00250 nopt++ ;
00251 if( nopt >= argc ) ANOVA_error("need argument after -session!") ;
00252 strcpy(option_data->session , argv[nopt++]) ;
00253 continue ;
00254 }
00255
00256
00257
00258 if (strncmp(argv[nopt], "-voxel", 6) == 0)
00259 {
00260 nopt++;
00261 if (nopt >= argc) ANOVA_error ("need argument after -voxel ");
00262 sscanf (argv[nopt], "%d", &ival);
00263 if (ival <= 0)
00264 ANOVA_error ("illegal argument after -voxel ");
00265 option_data->nvoxel = ival;
00266 nopt++;
00267 continue;
00268 }
00269
00270
00271
00272 if (strncmp(argv[nopt], "-type", 5) == 0)
00273 {
00274 nopt++;
00275 if (nopt >= argc) ANOVA_error ("need argument after -type ");
00276 sscanf (argv[nopt], "%d", &ival);
00277 if ((ival < 1) || (ival > 3))
00278 ANOVA_error ("illegal argument after -type ");
00279 option_data->model = ival;
00280 nopt++;
00281 continue;
00282 }
00283
00284
00285
00286 if (strncmp(argv[nopt], "-alevels", 5) == 0)
00287 {
00288 nopt++;
00289 if (nopt >= argc) ANOVA_error ("need argument after -alevels ");
00290 sscanf (argv[nopt], "%d", &ival);
00291 if ((ival <= 0) || (ival > MAX_LEVELS))
00292 ANOVA_error ("illegal argument after -alevels ");
00293 option_data->a = ival;
00294 nopt++;
00295 continue;
00296 }
00297
00298
00299
00300 if (strncmp(argv[nopt], "-blevels", 5) == 0)
00301 {
00302 nopt++;
00303 if (nopt >= argc) ANOVA_error ("need argument after -blevels ");
00304 sscanf (argv[nopt], "%d", &ival);
00305 if ((ival <= 0) || (ival > MAX_LEVELS))
00306 ANOVA_error ("illegal argument after -blevels ");
00307 option_data->b = ival;
00308 nopt++;
00309 continue;
00310 }
00311
00312
00313
00314 if (strncmp(argv[nopt], "-dset", 5) == 0)
00315 {
00316 nopt++;
00317 if (nopt+2 >= argc) ANOVA_error ("need 3 arguments after -dset ");
00318 sscanf (argv[nopt], "%d", &ival);
00319 if ((ival <= 0) || (ival > option_data->a))
00320 ANOVA_error ("illegal argument after -dset ");
00321
00322 nopt++;
00323 sscanf (argv[nopt], "%d", &jval);
00324 if ((jval <= 0) || (jval > option_data->b))
00325 ANOVA_error ("illegal argument after -dset ");
00326
00327 n[ival-1][jval-1] += 1;
00328 nij = n[ival-1][jval-1];
00329 if (nij > MAX_OBSERVATIONS)
00330 ANOVA_error ("too many data files");
00331
00332
00333 nopt++;
00334 dset = THD_open_dataset( argv[nopt] ) ;
00335 if( ! ISVALID_3DIM_DATASET(dset) )
00336 {
00337 sprintf(message,"Unable to open dataset file %s\n",
00338 argv[nopt]);
00339 ANOVA_error (message);
00340 }
00341
00342
00343 if (DSET_NVALS(dset) != 1)
00344 {
00345 sprintf(message,"Must specify exactly 1 sub-brick for file %s\n",
00346 argv[nopt]);
00347 ANOVA_error (message);
00348 }
00349
00350 THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00351
00352 option_data->xname[ival-1][jval-1][0][nij-1]
00353 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00354 strcpy (option_data->xname[ival-1][jval-1][0][nij-1], argv[nopt]);
00355 nopt++;
00356 continue;
00357 }
00358
00359
00360
00361 if (strncmp(argv[nopt], "-ftr", 5) == 0)
00362 {
00363 nopt++;
00364 if (nopt >= argc) ANOVA_error ("need argument after -ftr ");
00365 option_data->nftr = 1;
00366 option_data->ftrname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00367 strcpy (option_data->ftrname, argv[nopt]);
00368 nopt++;
00369 continue;
00370 }
00371
00372
00373
00374 if (strncmp(argv[nopt], "-fa", 5) == 0)
00375 {
00376 nopt++;
00377 if (nopt >= argc) ANOVA_error ("need argument after -fa ");
00378 option_data->nfa = 1;
00379 option_data->faname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00380 strcpy (option_data->faname, argv[nopt]);
00381 nopt++;
00382 continue;
00383 }
00384
00385
00386
00387 if (strncmp(argv[nopt], "-fb", 5) == 0)
00388 {
00389 nopt++;
00390 if (nopt >= argc) ANOVA_error ("need argument after -fb ");
00391 option_data->nfb = 1;
00392 option_data->fbname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00393 strcpy (option_data->fbname, argv[nopt]);
00394 nopt++;
00395 continue;
00396 }
00397
00398
00399
00400 if (strncmp(argv[nopt], "-fab", 5) == 0)
00401 {
00402 nopt++;
00403 if (nopt >= argc) ANOVA_error ("need argument after -fab ");
00404 option_data->nfab = 1;
00405 option_data->fabname = malloc (sizeof(char) * MAX_NAME_LENGTH);
00406 strcpy (option_data->fabname, argv[nopt]);
00407 nopt++;
00408 continue;
00409 }
00410
00411
00412
00413 if (strncmp(argv[nopt], "-amean", 5) == 0)
00414 {
00415 nopt++;
00416 if (nopt+1 >= argc) ANOVA_error ("need 2 arguments after -amean ");
00417
00418 option_data->num_ameans++;
00419 if (option_data->num_ameans > MAX_MEANS)
00420 ANOVA_error ("too many factor A level mean estimates");
00421
00422 sscanf (argv[nopt], "%d", &ival);
00423 if ((ival <= 0) || (ival > option_data->a))
00424 ANOVA_error ("illegal argument after -amean ");
00425 option_data->ameans[option_data->num_ameans-1] = ival - 1;
00426 nopt++;
00427
00428 option_data->amname[option_data->num_ameans-1]
00429 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00430 strcpy (option_data->amname[option_data->num_ameans-1], argv[nopt]);
00431 nopt++;
00432 continue;
00433 }
00434
00435
00436
00437 if (strncmp(argv[nopt], "-bmean", 5) == 0)
00438 {
00439 nopt++;
00440 if (nopt+1 >= argc) ANOVA_error ("need 2 arguments after -bmean ");
00441
00442 option_data->num_bmeans++;
00443 if (option_data->num_bmeans > MAX_MEANS)
00444 ANOVA_error ("too many factor B level mean estimates");
00445
00446 sscanf (argv[nopt], "%d", &ival);
00447 if ((ival <= 0) || (ival > option_data->b))
00448 ANOVA_error ("illegal argument after -bmean ");
00449 option_data->bmeans[option_data->num_bmeans-1] = ival - 1;
00450 nopt++;
00451
00452 option_data->bmname[option_data->num_bmeans-1]
00453 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00454 strcpy (option_data->bmname[option_data->num_bmeans-1], argv[nopt]);
00455 nopt++;
00456 continue;
00457 }
00458
00459
00460
00461 if (strncmp(argv[nopt], "-xmean", 5) == 0)
00462 {
00463 nopt++;
00464 if (nopt+2 >= argc) ANOVA_error ("need 3 arguments after -xmean ");
00465
00466 option_data->num_xmeans++;
00467 if (option_data->num_xmeans > MAX_MEANS)
00468 ANOVA_error ("too many cell mean estimates");
00469
00470 sscanf (argv[nopt], "%d", &ival);
00471 if ((ival <= 0) || (ival > option_data->a))
00472 ANOVA_error ("illegal argument after -xmean ");
00473 option_data->xmeans[option_data->num_xmeans-1][0] = ival - 1;
00474 nopt++;
00475
00476 sscanf (argv[nopt], "%d", &ival);
00477 if ((ival <= 0) || (ival > option_data->b))
00478 ANOVA_error ("illegal argument after -xmean ");
00479 option_data->xmeans[option_data->num_xmeans-1][1] = ival - 1;
00480 nopt++;
00481
00482 option_data->xmname[option_data->num_xmeans-1]
00483 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00484 strcpy (option_data->xmname[option_data->num_xmeans-1], argv[nopt]);
00485 nopt++;
00486 continue;
00487 }
00488
00489
00490
00491 if (strncmp(argv[nopt], "-adiff", 5) == 0)
00492 {
00493 nopt++;
00494 if (nopt+2 >= argc) ANOVA_error ("need 3 arguments after -adiff ");
00495
00496 option_data->num_adiffs++;
00497 if (option_data->num_adiffs > MAX_DIFFS)
00498 ANOVA_error ("too many factor A level differences");
00499
00500 sscanf (argv[nopt], "%d", &ival);
00501 if ((ival <= 0) || (ival > option_data->a))
00502 ANOVA_error ("illegal argument after -adiff ");
00503 option_data->adiffs[option_data->num_adiffs-1][0] = ival - 1;
00504 nopt++;
00505
00506 sscanf (argv[nopt], "%d", &ival);
00507 if ((ival <= 0) || (ival > option_data->a))
00508 ANOVA_error ("illegal argument after -adiff ");
00509 option_data->adiffs[option_data->num_adiffs-1][1] = ival - 1;
00510 nopt++;
00511
00512 option_data->adname[option_data->num_adiffs-1]
00513 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00514 strcpy (option_data->adname[option_data->num_adiffs-1], argv[nopt]);
00515 nopt++;
00516 continue;
00517 }
00518
00519
00520
00521 if (strncmp(argv[nopt], "-bdiff", 5) == 0)
00522 {
00523 nopt++;
00524 if (nopt+2 >= argc) ANOVA_error ("need 3 arguments after -bdiff ");
00525
00526 option_data->num_bdiffs++;
00527 if (option_data->num_bdiffs > MAX_DIFFS)
00528 ANOVA_error ("too many factor B level differences");
00529
00530 sscanf (argv[nopt], "%d", &ival);
00531 if ((ival <= 0) || (ival > option_data->b))
00532 ANOVA_error ("illegal argument after -bdiff ");
00533 option_data->bdiffs[option_data->num_bdiffs-1][0] = ival - 1;
00534 nopt++;
00535
00536 sscanf (argv[nopt], "%d", &ival);
00537 if ((ival <= 0) || (ival > option_data->b))
00538 ANOVA_error ("illegal argument after -bdiff ");
00539 option_data->bdiffs[option_data->num_bdiffs-1][1] = ival - 1;
00540 nopt++;
00541
00542 option_data->bdname[option_data->num_bdiffs-1]
00543 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00544 strcpy (option_data->bdname[option_data->num_bdiffs-1], argv[nopt]);
00545 nopt++;
00546 continue;
00547 }
00548
00549
00550
00551 if (strncmp(argv[nopt], "-xdiff", 5) == 0)
00552 {
00553 nopt++;
00554 if (nopt+4 >= argc) ANOVA_error ("need 5 arguments after -xdiff ");
00555
00556 option_data->num_xdiffs++;
00557 if (option_data->num_xdiffs > MAX_DIFFS)
00558 ANOVA_error ("too many cell means differences");
00559
00560 sscanf (argv[nopt], "%d", &ival);
00561 if ((ival <= 0) || (ival > option_data->a))
00562 ANOVA_error ("illegal argument after -xdiff ");
00563 option_data->xdiffs[option_data->num_xdiffs-1][0][0] = ival - 1;
00564 nopt++;
00565
00566 sscanf (argv[nopt], "%d", &ival);
00567 if ((ival <= 0) || (ival > option_data->b))
00568 ANOVA_error ("illegal argument after -xdiff ");
00569 option_data->xdiffs[option_data->num_xdiffs-1][0][1] = ival - 1;
00570 nopt++;
00571
00572 sscanf (argv[nopt], "%d", &ival);
00573 if ((ival <= 0) || (ival > option_data->a))
00574 ANOVA_error ("illegal argument after -xdiff ");
00575 option_data->xdiffs[option_data->num_xdiffs-1][1][0] = ival - 1;
00576 nopt++;
00577
00578 sscanf (argv[nopt], "%d", &ival);
00579 if ((ival <= 0) || (ival > option_data->b))
00580 ANOVA_error ("illegal argument after -xdiff ");
00581 option_data->xdiffs[option_data->num_xdiffs-1][1][1] = ival - 1;
00582 nopt++;
00583
00584 option_data->xdname[option_data->num_xdiffs-1]
00585 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00586 strcpy (option_data->xdname[option_data->num_xdiffs-1], argv[nopt]);
00587 nopt++;
00588 continue;
00589 }
00590
00591
00592
00593 if (strncmp(argv[nopt], "-acontr", 5) == 0)
00594 {
00595 nopt++;
00596 if (nopt + option_data->a >= argc)
00597 ANOVA_error ("need a+1 arguments after -acontr ");
00598
00599 option_data->num_acontr++;
00600 if (option_data->num_acontr > MAX_CONTR)
00601 ANOVA_error ("too many factor A level contrasts");
00602
00603 for (i = 0; i < option_data->a; i++)
00604 {
00605 sscanf (argv[nopt], "%f", &fval);
00606 option_data->acontr[option_data->num_acontr - 1][i] = fval ;
00607 nopt++;
00608 }
00609
00610 option_data->acname[option_data->num_acontr-1]
00611 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00612 strcpy (option_data->acname[option_data->num_acontr-1], argv[nopt]);
00613 nopt++;
00614 continue;
00615 }
00616
00617
00618
00619 if (strncmp(argv[nopt], "-bcontr", 5) == 0)
00620 {
00621 nopt++;
00622 if (nopt + option_data->b >= argc)
00623 ANOVA_error ("need b+1 arguments after -bcontr ");
00624
00625 option_data->num_bcontr++;
00626 if (option_data->num_bcontr > MAX_CONTR)
00627 ANOVA_error ("too many factor B level contrasts");
00628
00629 for (i = 0; i < option_data->b; i++)
00630 {
00631 sscanf (argv[nopt], "%f", &fval);
00632 option_data->bcontr[option_data->num_bcontr - 1][i] = fval ;
00633 nopt++;
00634 }
00635
00636 option_data->bcname[option_data->num_bcontr-1]
00637 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00638 strcpy (option_data->bcname[option_data->num_bcontr-1], argv[nopt]);
00639 nopt++;
00640 continue;
00641 }
00642
00643
00644
00645 if (strncmp(argv[nopt], "-xcontr", 5) == 0)
00646 {
00647 nopt++;
00648 if (nopt + (option_data->a * option_data->b) >= argc)
00649 ANOVA_error ("need ab+1 arguments after -xcontr ");
00650
00651 option_data->num_xcontr++;
00652 if (option_data->num_xcontr > MAX_CONTR)
00653 ANOVA_error ("too many cell means contrasts");
00654
00655 for (i = 0; i < option_data->a; i++)
00656 for (j = 0; j < option_data->b; j++)
00657 {
00658 sscanf (argv[nopt], "%f", &fval);
00659 option_data->xcontr[option_data->num_xcontr - 1][i][j] = fval ;
00660 nopt++;
00661 }
00662
00663 option_data->xcname[option_data->num_xcontr-1]
00664 = malloc (sizeof(char) * MAX_NAME_LENGTH);
00665 strcpy (option_data->xcname[option_data->num_xcontr-1], argv[nopt]);
00666 nopt++;
00667 continue;
00668 }
00669
00670
00671
00672 if (strncmp(argv[nopt], "-bucket", 4) == 0)
00673 {
00674 nopt++;
00675 if (nopt >= argc) ANOVA_error ("need argument after -bucket ");
00676 option_data->bucket_filename = malloc (sizeof(char)*MAX_NAME_LENGTH);
00677 strcpy (option_data->bucket_filename, argv[nopt]);
00678 nopt++;
00679 continue;
00680 }
00681
00682
00683
00684 sprintf (message,"Unrecognized command line option: %s\n", argv[nopt]);
00685 ANOVA_error (message);
00686 }
00687
00688
00689
00690 option_data->n = n[0][0];
00691 for (i = 0; i < option_data->a; i++)
00692 for (j = 0; j < option_data->b; j++)
00693 if (n[i][j] != option_data->n)
00694 ANOVA_error ("must have equal sample sizes for 3dANOVA2");
00695 }
00696
00697
00698
00699
00700
00701
00702
00703 void check_temporary_files (anova_options * option_data)
00704 {
00705
00706 check_one_temporary_file ("ss0");
00707 check_one_temporary_file ("ssi");
00708 check_one_temporary_file ("ssj");
00709 check_one_temporary_file ("ssij");
00710 check_one_temporary_file ("ssijk");
00711
00712 check_one_temporary_file ("sse");
00713 check_one_temporary_file ("sstr");
00714 check_one_temporary_file ("ssa");
00715 check_one_temporary_file ("ssb");
00716 check_one_temporary_file ("ssab");
00717 }
00718
00719
00720
00721
00722
00723
00724
00725 void check_output_files (anova_options * option_data)
00726 {
00727 int i;
00728
00729 if (option_data->nftr > 0)
00730 check_one_output_file (option_data, option_data->ftrname);
00731
00732 if (option_data->nfa > 0)
00733 check_one_output_file (option_data, option_data->faname);
00734
00735 if (option_data->nfb > 0)
00736 check_one_output_file (option_data, option_data->fbname);
00737
00738 if (option_data->nfab > 0)
00739 check_one_output_file (option_data, option_data->fabname);
00740
00741 if (option_data->num_ameans > 0)
00742 for (i = 0; i < option_data->num_ameans; i++)
00743 check_one_output_file (option_data, option_data->amname[i]);
00744
00745 if (option_data->num_bmeans > 0)
00746 for (i = 0; i < option_data->num_bmeans; i++)
00747 check_one_output_file (option_data, option_data->bmname[i]);
00748
00749 if (option_data->num_xmeans > 0)
00750 for (i = 0; i < option_data->num_xmeans; i++)
00751 check_one_output_file (option_data, option_data->xmname[i]);
00752
00753 if (option_data->num_adiffs > 0)
00754 for (i = 0; i < option_data->num_adiffs; i++)
00755 check_one_output_file (option_data, option_data->adname[i]);
00756
00757 if (option_data->num_bdiffs > 0)
00758 for (i = 0; i < option_data->num_bdiffs; i++)
00759 check_one_output_file (option_data, option_data->bdname[i]);
00760
00761 if (option_data->num_xdiffs > 0)
00762 for (i = 0; i < option_data->num_xdiffs; i++)
00763 check_one_output_file (option_data, option_data->xdname[i]);
00764
00765 if (option_data->num_acontr > 0)
00766 for (i = 0; i < option_data->num_acontr; i++)
00767 check_one_output_file (option_data, option_data->acname[i]);
00768
00769 if (option_data->num_bcontr > 0)
00770 for (i = 0; i < option_data->num_bcontr; i++)
00771 check_one_output_file (option_data, option_data->bcname[i]);
00772
00773 if (option_data->num_xcontr > 0)
00774 for (i = 0; i < option_data->num_xcontr; i++)
00775 check_one_output_file (option_data, option_data->xcname[i]);
00776
00777 if (option_data->bucket_filename != NULL)
00778 check_one_output_file (option_data, option_data->bucket_filename);
00779
00780 }
00781
00782
00783
00784
00785
00786
00787
00788 void check_for_valid_inputs (anova_options * option_data)
00789 {
00790 int a, b;
00791 int n;
00792
00793
00794
00795 a = option_data->a;
00796 b = option_data->b;
00797 n = option_data->n;
00798
00799
00800
00801 if (a < 2) ANOVA_error ("must specify number of factor A levels (a>1) ");
00802 if (b < 2) ANOVA_error ("must specify number of factor B levels (b>1) ");
00803 if (n < 1) ANOVA_error ("sample size is too small");
00804
00805 switch (option_data->model)
00806 {
00807 case 1:
00808 if (n == 1)
00809 ANOVA_error ("sample size is too small for fixed effects model");
00810 break;
00811 case 2:
00812 if (option_data->nftr > 0)
00813 ANOVA_error ("-ftr is inappropriate for random effects model");
00814 if (option_data->num_ameans > 0)
00815 ANOVA_error ("-amean is inappropriate for random effects model");
00816 if (option_data->num_bmeans > 0)
00817 ANOVA_error ("-bmean is inappropriate for random effects model");
00818 if (option_data->num_xmeans > 0)
00819 ANOVA_error ("-xmean is inappropriate for random effects model");
00820 if (option_data->num_adiffs > 0)
00821 ANOVA_error ("-adiff is inappropriate for random effects model");
00822 if (option_data->num_bdiffs > 0)
00823 ANOVA_error ("-bdiff is inappropriate for random effects model");
00824 if (option_data->num_xdiffs > 0)
00825 ANOVA_error ("-xdiff is inappropriate for random effects model");
00826 if (option_data->num_acontr > 0)
00827 ANOVA_error ("-acontr is inappropriate for random effects model");
00828 if (option_data->num_bcontr > 0)
00829 ANOVA_error ("-bcontr is inappropriate for random effects model");
00830 if (option_data->num_xcontr > 0)
00831 ANOVA_error ("-xcontr is inappropriate for random effects model");
00832 if ((n == 1) && (option_data->nfab > 0))
00833 ANOVA_error ("sample size too small to calculate F-interaction");
00834 break;
00835 case 3:
00836 if (option_data->nftr > 0)
00837 ANOVA_error ("-ftr is inappropriate for mixed effects model");
00838 if (option_data->num_bmeans > 0)
00839 ANOVA_error ("-bmean is inappropriate for mixed effects model");
00840 if (option_data->num_xmeans > 0)
00841 ANOVA_error ("-xmean is inappropriate for mixed effects model");
00842 if (option_data->num_bdiffs > 0)
00843 ANOVA_error ("-bdiff is inappropriate for mixed effects model");
00844 if (option_data->num_xdiffs > 0)
00845 ANOVA_error ("-xdiff is inappropriate for mixed effects model");
00846 if (option_data->num_bcontr > 0)
00847 ANOVA_error ("-bcontr is inappropriate for mixed effects model");
00848 if (option_data->num_xcontr > 0)
00849 ANOVA_error ("-xcontr is inappropriate for mixed effects model");
00850 if ((n == 1) && (option_data->nfab > 0))
00851 ANOVA_error ("sample size too small to calculate F-interaction");
00852 if ((n == 1) && (option_data->nfb > 0))
00853 ANOVA_error ("sample size too small to calculate F for B effect");
00854 break;
00855 }
00856
00857 }
00858
00859
00860
00861
00862
00863
00864
00865
00866 int required_data_files (anova_options * option_data)
00867 {
00868 int now;
00869
00870 int nout;
00871 int nmax;
00872
00873
00874 if (option_data->n != 1)
00875 {
00876 nmax = 6; now = 5;
00877 }
00878 else
00879 {
00880 nmax = 5; now = 5;
00881 }
00882
00883 nout = option_data->nftr + option_data->nfab
00884 + option_data->nfa + option_data->nfb
00885 + option_data->num_ameans + option_data->num_bmeans + option_data->num_xmeans
00886 + option_data->num_adiffs + option_data->num_bdiffs + option_data->num_xdiffs
00887 + option_data->num_acontr + option_data->num_bcontr
00888 + option_data->num_xcontr;
00889
00890 now = now + nout;
00891
00892 nmax = max (now, nmax);
00893
00894 if (option_data->bucket_filename != NULL)
00895 nmax = max (nmax, 2*nout);
00896
00897 return (nmax);
00898 }
00899
00900
00901
00902
00903
00904
00905
00906 void initialize (int argc, char ** argv, anova_options ** option_data)
00907 {
00908 int a, b;
00909 int n;
00910
00911
00912
00913 commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00914
00915
00916
00917 *option_data = (anova_options *) malloc(sizeof(anova_options));
00918 if (*option_data == NULL)
00919 ANOVA_error ("memory allocation error");
00920
00921
00922 get_options(argc, argv, *option_data);
00923
00924
00925 (*option_data)->first_dataset = (*option_data)->xname[0][0][0][0];
00926 get_dimensions (*option_data);
00927 printf ("Data set dimensions: nx = %d ny = %d nz = %d nxyz = %d \n",
00928 (*option_data)->nx, (*option_data)->ny,
00929 (*option_data)->nz, (*option_data)->nxyz);
00930 if ((*option_data)->nvoxel > (*option_data)->nxyz)
00931 ANOVA_error ("argument of -voxel is too large");
00932
00933
00934 a = (*option_data)->a;
00935 b = (*option_data)->b;
00936 n = (*option_data)->n;
00937
00938
00939 (*option_data)->nt = n * a * b;
00940
00941
00942 check_for_valid_inputs (*option_data);
00943
00944
00945 check_temporary_files (*option_data);
00946
00947
00948 check_output_files (*option_data);
00949
00950
00951 if ((*option_data)->diskspace) check_disk_space (*option_data);
00952 }
00953
00954
00955
00956
00957
00958
00959
00960
00961 void calculate_sum (anova_options * option_data,
00962 int ii, int jj, float * ysum)
00963 {
00964 float * y = NULL;
00965 int i, itop, ibot;
00966 int j, jtop, jbot;
00967 int m;
00968 int a;
00969 int b;
00970 int n;
00971 int ixyz, nxyz;
00972 int nvoxel;
00973 char sum_label[MAX_NAME_LENGTH];
00974 char str[MAX_NAME_LENGTH];
00975
00976
00977
00978 a = option_data->a;
00979 b = option_data->b;
00980 n = option_data->n;
00981 nxyz = option_data->nxyz;
00982 nvoxel = option_data->nvoxel;
00983
00984
00985 y = (float *) malloc(sizeof(float)*nxyz);
00986 if (y == NULL) ANOVA_error ("unable to allocate sufficient memory");
00987
00988
00989
00990 if (ii < 0)
00991 { ibot = 0; itop = a; }
00992 else
00993 { ibot = ii; itop = ii+1; }
00994
00995 if (jj < 0)
00996 { jbot = 0; jtop = b; }
00997 else
00998 { jbot = jj; jtop = jj+1; }
00999
01000
01001 volume_zero (ysum, nxyz);
01002
01003
01004 for (i = ibot; i < itop; i++)
01005 {
01006
01007 for (j = jbot; j < jtop; j++)
01008 {
01009
01010 for (m = 0; m < n; m++)
01011 {
01012 read_afni_data (option_data,
01013 option_data->xname[i][j][0][m], y);
01014 if (nvoxel > 0)
01015 printf ("y[%d][%d][%d] = %f \n",
01016 i+1, j+1, m+1, y[nvoxel-1]);
01017 for (ixyz = 0; ixyz < nxyz; ixyz++)
01018 ysum[ixyz] += y[ixyz];
01019 }
01020 }
01021 }
01022
01023
01024
01025 if (nvoxel > 0)
01026 {
01027 strcpy (sum_label, "y");
01028 if (ii < 0)
01029 strcat (sum_label, "[.]");
01030 else
01031 {
01032 sprintf (str, "[%d]", ii+1);
01033 strcat (sum_label, str);
01034 }
01035 if (jj < 0)
01036 strcat (sum_label, "[.]");
01037 else
01038 {
01039 sprintf (str, "[%d]", jj+1);
01040 strcat (sum_label, str);
01041 }
01042 printf ("%s[.] = %f \n", sum_label, ysum[nvoxel-1]);
01043 }
01044
01045
01046 free (y); y = NULL;
01047
01048 }
01049
01050
01051
01052
01053
01054
01055
01056
01057
01058
01059
01060
01061
01062
01063
01064
01065 void calculate_t_from_sums(float * result, float * mean, float * sum_sq,
01066 int df, int nvox)
01067 {
01068 const float EPSILON = 1.0e-10;
01069 double dval;
01070 float *res, *mp, *ssp, sval;
01071 int index, df_prod;
01072
01073
01074 df_prod = df * (df + 1);
01075
01076 res = result;
01077 mp = mean;
01078 ssp = sum_sq;
01079 for (index = 0; index < nvox; index++)
01080 {
01081 sval = *mp;
01082 dval = *ssp - (df+1.0) * sval * sval;
01083
01084 if (dval < EPSILON) *res = 0.0;
01085 else *res = sval * sqrt( df_prod / dval );
01086
01087 res++; mp++; ssp++;
01088 }
01089 }
01090
01091
01092
01093
01094
01095
01096
01097
01098
01099
01100 void calculate_sum_sq (anova_options * option_data,
01101 int ii, int jj, float * ysum)
01102 {
01103 double * yd = NULL;
01104 float * y = NULL;
01105 int i, itop, ibot;
01106 int j, jtop, jbot;
01107 int m;
01108 int a;
01109 int b;
01110 int n;
01111 int ixyz, nxyz;
01112 int nvoxel;
01113 char sum_label[MAX_NAME_LENGTH];
01114 char str[MAX_NAME_LENGTH];
01115
01116
01117
01118 a = option_data->a;
01119 b = option_data->b;
01120 n = option_data->n;
01121 nxyz = option_data->nxyz;
01122 nvoxel = option_data->nvoxel;
01123
01124
01125 y = (float *) malloc(sizeof(float)*nxyz);
01126 yd = (double *) malloc(sizeof(double)*nxyz);
01127 if (!y || !yd) ANOVA_error ("sum_sq: unable to allocate sufficient memory");
01128
01129
01130 if (ii < 0) { ibot = 0; itop = a; }
01131 else { ibot = ii; itop = ii+1; }
01132
01133 if (jj < 0) { jbot = 0; jtop = b; }
01134 else { jbot = jj; jtop = jj+1; }
01135
01136 for (ixyz = 0; ixyz < nxyz; ixyz++)
01137 yd[ixyz] = 0.0;
01138
01139
01140 for (i = ibot; i < itop; i++)
01141
01142 for (j = jbot; j < jtop; j++)
01143
01144 for (m = 0; m < n; m++)
01145 {
01146 read_afni_data (option_data, option_data->xname[i][j][0][m], y);
01147 if (nvoxel > 0)
01148 printf ("y[%d][%d][%d] = %f \n", i+1, j+1, m+1, y[nvoxel-1]);
01149 for (ixyz = 0; ixyz < nxyz; ixyz++)
01150 yd[ixyz] += y[ixyz] * y[ixyz];
01151 }
01152
01153
01154 for (ixyz = 0; ixyz < nxyz; ixyz++)
01155 ysum[ixyz] = yd[ixyz];
01156
01157
01158 if (nvoxel > 0)
01159 {
01160 strcpy (sum_label, "y");
01161 if (ii < 0)
01162 strcat (sum_label, "[.]");
01163 else
01164 {
01165 sprintf (str, "[%d]", ii+1);
01166 strcat (sum_label, str);
01167 }
01168 if (jj < 0)
01169 strcat (sum_label, "[.]");
01170 else
01171 {
01172 sprintf (str, "[%d]", jj+1);
01173 strcat (sum_label, str);
01174 }
01175 printf ("%s_squares[.] = %f \n", sum_label, ysum[nvoxel-1]);
01176 }
01177
01178
01179 free (y); y = NULL;
01180 free (yd); yd = NULL;
01181 }
01182
01183
01184
01185
01186
01187
01188
01189
01190
01191
01192
01193
01194 void calc_acontr_mean (anova_options *option_data, float *contr, float *cmean)
01195 {
01196 double * sum = NULL;
01197 int i;
01198 int a, b;
01199 int n;
01200 int ixyz, nxyz;
01201 int nvoxel;
01202
01203
01204
01205 a = option_data->a;
01206 b = option_data->b;
01207 n = option_data->n;
01208 nxyz = option_data->nxyz;
01209 nvoxel = option_data->nvoxel;
01210
01211
01212 sum = (double *) malloc(sizeof(double)*nxyz);
01213 if (sum == NULL)
01214 ANOVA_error ("calc_acontr_mean: unable to allocate sufficient memory");
01215
01216 for (ixyz = 0; ixyz < nxyz; ixyz++)
01217 sum[ixyz] = 0.0;
01218
01219
01220 for (i = 0; i < a; i++)
01221 {
01222 if (contr[i] == 0.0 ) continue;
01223
01224
01225 calculate_sum(option_data, i, -1, cmean);
01226 if (nvoxel > 0)
01227 printf( "acontr[%d] = %f, ymean = %f\n",
01228 i, contr[i], cmean[nvoxel-1] / (n*b) );
01229 for (ixyz = 0; ixyz < nxyz; ixyz++)
01230 sum[ixyz] += (double)cmean[ixyz] / (n*b) * contr[i];
01231 }
01232
01233
01234 for (ixyz = 0; ixyz < nxyz; ixyz++)
01235 cmean[ixyz] = sum[ixyz];
01236
01237
01238 free (sum);
01239 }
01240
01241
01242
01243
01244
01245
01246
01247
01248
01249
01250
01251
01252 void calc_sum_sq_acontr(anova_options *option_data, float *acontr, float *sum)
01253 {
01254 double * dsum = NULL;
01255 double * dcontr = NULL;
01256 int i, j;
01257 int k;
01258 int a, b;
01259 int n;
01260 int ixyz, nxyz;
01261 int nvoxel;
01262
01263
01264
01265 a = option_data->a;
01266 b = option_data->b;
01267 n = option_data->n;
01268 nxyz = option_data->nxyz;
01269 nvoxel = option_data->nvoxel;
01270
01271
01272 dsum = (double *) malloc(sizeof(double)*nxyz);
01273 dcontr = (double *) malloc(sizeof(double)*nxyz);
01274 if (dsum == NULL || dcontr == NULL)
01275 ANOVA_error ("calc_sum_sq_acontr: unable to allocate sufficient memory");
01276
01277 for (ixyz = 0; ixyz < nxyz; ixyz++)
01278 dsum[ixyz] = 0.0;
01279
01280
01281 for ( j = 0; j < b; j++ )
01282 for ( k = 0; k < n; k++ )
01283 {
01284
01285 for (ixyz = 0; ixyz < nxyz; ixyz++)
01286 dcontr[ixyz] = 0.0;
01287
01288 for (i = 0; i < a; i++)
01289 {
01290 if (acontr[i] == 0.0 ) continue;
01291
01292
01293 read_afni_data(option_data, option_data->xname[i][j][0][k], sum);
01294 for (ixyz = 0; ixyz < nxyz; ixyz++)
01295 dcontr[ixyz] += (double)sum[ixyz] * acontr[i];
01296 }
01297 for (ixyz = 0; ixyz < nxyz; ixyz++)
01298 dsum[ixyz] += dcontr[ixyz] * dcontr[ixyz];
01299 }
01300
01301
01302 for (ixyz = 0; ixyz < nxyz; ixyz++)
01303 sum[ixyz] = dsum[ixyz];
01304
01305
01306 free (dsum);
01307 free (dcontr);
01308 }
01309
01310
01311
01312
01313
01314
01315
01316
01317 void calculate_ss0 (anova_options * option_data)
01318 {
01319 float * ss0 = NULL;
01320 float * ysum = NULL;
01321 int a;
01322 int b;
01323 int n;
01324 int ixyz, nxyz;
01325 int nvoxel;
01326 int nval;
01327 char filename[MAX_NAME_LENGTH];
01328
01329
01330
01331 a = option_data->a;
01332 b = option_data->b;
01333 n = option_data->n;
01334 nxyz = option_data->nxyz;
01335 nvoxel = option_data->nvoxel;
01336 nval = a * b * n;
01337
01338
01339 ss0 = (float *) malloc(sizeof(float)*nxyz);
01340 ysum = (float *) malloc(sizeof(float)*nxyz);
01341 if ((ss0 == NULL) || (ysum == NULL))
01342 ANOVA_error ("unable to allocate sufficient memory");
01343
01344
01345 calculate_sum (option_data, -1, -1, ysum);
01346
01347
01348 for (ixyz = 0; ixyz < nxyz; ixyz++)
01349 ss0[ixyz] = ysum[ixyz] * ysum[ixyz] / nval;
01350
01351
01352
01353 if (nvoxel > 0)
01354 printf ("SS0 = %f \n", ss0[nvoxel-1]);
01355 strcpy (filename, "ss0");
01356 volume_write (filename, ss0, nxyz);
01357
01358
01359
01360 free (ysum); ysum = NULL;
01361 free (ss0); ss0 = NULL;
01362
01363 }
01364
01365
01366
01367
01368
01369
01370
01371
01372 void calculate_ssi (anova_options * option_data)
01373 {
01374 float * ssi = NULL;
01375 float * ysum = NULL;
01376 int a;
01377 int b;
01378 int n;
01379 int i;
01380 int ixyz, nxyz;
01381 int nvoxel;
01382 int nval;
01383 char filename[MAX_NAME_LENGTH];
01384
01385
01386
01387 a = option_data->a;
01388 b = option_data->b;
01389 n = option_data->n;
01390 nxyz = option_data->nxyz;
01391 nvoxel = option_data->nvoxel;
01392 nval = b * n;
01393
01394
01395 ssi = (float *) malloc(sizeof(float)*nxyz);
01396 ysum = (float *) malloc(sizeof(float)*nxyz);
01397 if ((ssi == NULL) || (ysum == NULL))
01398 ANOVA_error ("unable to allocate sufficient memory");
01399
01400 volume_zero (ssi, nxyz);
01401
01402
01403 for (i = 0; i < a; i++)
01404 {
01405
01406 calculate_sum (option_data, i, -1, ysum);
01407
01408
01409 for (ixyz = 0; ixyz < nxyz; ixyz++)
01410 ssi[ixyz] += ysum[ixyz] * ysum[ixyz] / nval;
01411 }
01412
01413
01414 if (nvoxel > 0)
01415 printf ("SSI = %f \n", ssi[nvoxel-1]);
01416 strcpy (filename, "ssi");
01417 volume_write (filename, ssi, nxyz);
01418
01419
01420
01421 free (ysum); ysum = NULL;
01422 free (ssi); ssi = NULL;
01423
01424 }
01425
01426
01427
01428
01429
01430
01431
01432
01433 void calculate_ssj (anova_options * option_data)
01434 {
01435 float * ssj = NULL;
01436 float * ysum = NULL;
01437 int a;
01438 int b;
01439 int n;
01440 int j;
01441 int ixyz, nxyz;
01442 int nvoxel;
01443 int nval;
01444 char filename[MAX_NAME_LENGTH];
01445
01446
01447
01448 a = option_data->a;
01449 b = option_data->b;
01450 n = option_data->n;
01451 nxyz = option_data->nxyz;
01452 nvoxel = option_data->nvoxel;
01453 nval = a * n;
01454
01455
01456 ssj = (float *) malloc(sizeof(float)*nxyz);
01457 ysum = (float *) malloc(sizeof(float)*nxyz);
01458 if ((ssj == NULL) || (ysum == NULL))
01459 ANOVA_error ("unable to allocate sufficient memory");
01460
01461 volume_zero (ssj, nxyz);
01462
01463
01464 for (j = 0; j < b; j++)
01465 {
01466
01467 calculate_sum (option_data, -1, j, ysum);
01468
01469
01470 for (ixyz = 0; ixyz < nxyz; ixyz++)
01471 ssj[ixyz] += ysum[ixyz] * ysum[ixyz] / nval;
01472 }
01473
01474
01475 if (nvoxel > 0)
01476 printf ("SSJ = %f \n", ssj[nvoxel-1]);
01477 strcpy (filename, "ssj");
01478 volume_write (filename, ssj, nxyz);
01479
01480
01481
01482 free (ysum); ysum = NULL;
01483 free (ssj); ssj = NULL;
01484
01485 }
01486
01487
01488
01489
01490
01491
01492
01493
01494 void calculate_ssij (anova_options * option_data)
01495 {
01496 float * ssij = NULL;
01497 float * ysum = NULL;
01498 int a;
01499 int b;
01500 int n;
01501 int i, j;
01502 int ixyz, nxyz;
01503 int nvoxel;
01504 int nval;
01505 char filename[MAX_NAME_LENGTH];
01506
01507
01508
01509 a = option_data->a;
01510 b = option_data->b;
01511 n = option_data->n;
01512 nxyz = option_data->nxyz;
01513 nvoxel = option_data->nvoxel;
01514 nval = n;
01515
01516
01517 ssij = (float *) malloc(sizeof(float)*nxyz);
01518 ysum = (float *) malloc(sizeof(float)*nxyz);
01519 if ((ssij == NULL) || (ysum == NULL))
01520 ANOVA_error ("unable to allocate sufficient memory");
01521
01522 volume_zero (ssij, nxyz);
01523
01524
01525 for (i = 0; i < a; i++)
01526 {
01527
01528 for (j = 0; j < b; j++)
01529 {
01530
01531 calculate_sum (option_data, i, j, ysum);
01532
01533
01534 for (ixyz = 0; ixyz < nxyz; ixyz++)
01535 ssij[ixyz] += ysum[ixyz] * ysum[ixyz] / nval;
01536 }
01537 }
01538
01539
01540 if (nvoxel > 0)
01541 printf ("SSIJ = %f \n", ssij[nvoxel-1]);
01542 strcpy (filename, "ssij");
01543 volume_write (filename, ssij, nxyz);
01544
01545
01546
01547 free (ysum); ysum = NULL;
01548 free (ssij); ssij = NULL;
01549
01550 }
01551
01552
01553
01554
01555
01556
01557
01558
01559 void calculate_ssijk (anova_options * option_data)
01560 {
01561 float * ssijk = NULL;
01562 float * y = NULL;
01563 int i;
01564 int j;
01565 int m;
01566 int a;
01567 int b;
01568 int n;
01569 int ixyz, nxyz;
01570 int nvoxel;
01571
01572
01573
01574 a = option_data->a;
01575 b = option_data->b;
01576 n = option_data->n;
01577 nxyz = option_data->nxyz;
01578 nvoxel = option_data->nvoxel;
01579
01580
01581 ssijk = (float *) malloc(sizeof(float)*nxyz);
01582 y = (float *) malloc(sizeof(float)*nxyz);
01583 if ((ssijk == NULL) || (y == NULL))
01584 ANOVA_error ("unable to allocate sufficient memory");
01585
01586
01587 volume_zero (ssijk, nxyz);
01588
01589 for (i = 0; i < a; i++)
01590 {
01591 for (j = 0; j < b; j++)
01592 {
01593 for (m = 0; m < n; m++)
01594 {
01595 read_afni_data (option_data,
01596 option_data->xname[i][j][0][m], y);
01597
01598 for (ixyz = 0; ixyz < nxyz; ixyz++)
01599 ssijk[ixyz] += y[ixyz] * y[ixyz];
01600 }
01601 }
01602 }
01603
01604
01605
01606 if (nvoxel > 0)
01607 printf ("SSIJK = %f \n", ssijk[nvoxel-1]);
01608 volume_write ("ssijk", ssijk, nxyz);
01609
01610
01611 free (y); y = NULL;
01612 free (ssijk); ssijk = NULL;
01613
01614 }
01615
01616
01617
01618
01619
01620
01621
01622
01623 void calculate_sse (anova_options * option_data)
01624 {
01625 float * y = NULL;
01626 float * sse = NULL;
01627 int ixyz, nxyz;
01628 int nvoxel;
01629
01630
01631
01632 nxyz = option_data->nxyz;
01633 nvoxel = option_data->nvoxel;
01634
01635
01636 sse = (float *) malloc (sizeof(float)*nxyz);
01637 y = (float *) malloc (sizeof(float)*nxyz);
01638 if ((y == NULL) || (sse == NULL))
01639 ANOVA_error ("unable to allocate sufficient memory");
01640
01641
01642
01643 volume_read ("ssijk", sse, nxyz);
01644
01645 volume_read ("ssij", y, nxyz);
01646 for (ixyz = 0; ixyz < nxyz; ixyz++)
01647 sse[ixyz] -= y[ixyz];
01648
01649
01650
01651 for (ixyz = 0; ixyz < nxyz; ixyz++)
01652 if (sse[ixyz] < 0.0) sse[ixyz] = 0.0;
01653
01654
01655 if (nvoxel > 0)
01656 printf ("SSE = %f \n", sse[nvoxel-1]);
01657 volume_write ("sse", sse, nxyz);
01658
01659
01660 free (y); y = NULL;
01661 free (sse); sse = NULL;
01662
01663 }
01664
01665
01666
01667
01668
01669
01670
01671
01672 void calculate_sstr (anova_options * option_data)
01673 {
01674 float * y = NULL;
01675 float * sstr = NULL;
01676 int ixyz, nxyz;
01677 int nvoxel;
01678
01679
01680
01681 nxyz = option_data->nxyz;
01682 nvoxel = option_data->nvoxel;
01683
01684
01685 sstr = (float *) malloc (sizeof(float)*nxyz);
01686 y = (float *) malloc (sizeof(float)*nxyz);
01687 if ((y == NULL) || (sstr == NULL))
01688 ANOVA_error ("unable to allocate sufficient memory");
01689
01690
01691
01692 volume_read ("ssij", sstr, nxyz);
01693
01694 volume_read ("ss0", y, nxyz);
01695 for (ixyz = 0; ixyz < nxyz; ixyz++)
01696 sstr[ixyz] -= y[ixyz];
01697
01698
01699
01700 for (ixyz = 0; ixyz < nxyz; ixyz++)
01701 if (sstr[ixyz] < 0.0) sstr[ixyz] = 0.0;
01702
01703
01704 if (nvoxel > 0)
01705 printf ("SSTR = %f \n", sstr[nvoxel-1]);
01706 volume_write ("sstr", sstr, nxyz);
01707
01708
01709 free (y); y = NULL;
01710 free (sstr); sstr = NULL;
01711
01712 }
01713
01714
01715
01716
01717
01718
01719
01720
01721 void calculate_ssa (anova_options * option_data)
01722 {
01723 float * y = NULL;
01724 float * ssa = NULL;
01725 int ixyz, nxyz;
01726 int nvoxel;
01727
01728
01729
01730 nxyz = option_data->nxyz;
01731 nvoxel = option_data->nvoxel;
01732
01733
01734 ssa = (float *) malloc (sizeof(float)*nxyz);
01735 y = (float *) malloc (sizeof(float)*nxyz);
01736 if ((y == NULL) || (ssa == NULL))
01737 ANOVA_error ("unable to allocate sufficient memory");
01738
01739
01740
01741 volume_read ("ssi", ssa, nxyz);
01742
01743 volume_read ("ss0", y, nxyz);
01744 for (ixyz = 0; ixyz < nxyz; ixyz++)
01745 ssa[ixyz] -= y[ixyz];
01746
01747
01748
01749 for (ixyz = 0; ixyz < nxyz; ixyz++)
01750 if (ssa[ixyz] < 0.0) ssa[ixyz] = 0.0;
01751
01752
01753 if (nvoxel > 0)
01754 printf ("SSA = %f \n", ssa[nvoxel-1]);
01755 volume_write ("ssa", ssa, nxyz);
01756
01757
01758 free (y); y = NULL;
01759 free (ssa); ssa = NULL;
01760
01761 }
01762
01763
01764
01765
01766
01767
01768
01769
01770 void calculate_ssb (anova_options * option_data)
01771 {
01772 float * y = NULL;
01773 float * ssb = NULL;
01774 int ixyz, nxyz;
01775 int nvoxel;
01776
01777
01778
01779 nxyz = option_data->nxyz;
01780 nvoxel = option_data->nvoxel;
01781
01782
01783 ssb = (float *) malloc (sizeof(float)*nxyz);
01784 y = (float *) malloc (sizeof(float)*nxyz);
01785 if ((y == NULL) || (ssb == NULL))
01786 ANOVA_error ("unable to allocate sufficient memory");
01787
01788
01789
01790 volume_read ("ssj", ssb, nxyz);
01791
01792 volume_read ("ss0", y, nxyz);
01793 for (ixyz = 0; ixyz < nxyz; ixyz++)
01794 ssb[ixyz] -= y[ixyz];
01795
01796
01797
01798 for (ixyz = 0; ixyz < nxyz; ixyz++)
01799 if (ssb[ixyz] < 0.0) ssb[ixyz] = 0.0;
01800
01801
01802 if (nvoxel > 0)
01803 printf ("SSB = %f \n", ssb[nvoxel-1]);
01804 volume_write ("ssb", ssb, nxyz);
01805
01806
01807 free (y); y = NULL;
01808 free (ssb); ssb = NULL;
01809
01810 }
01811
01812
01813
01814
01815
01816
01817
01818
01819 void calculate_ssab (anova_options * option_data)
01820 {
01821 float * y = NULL;
01822 float * ssab = NULL;
01823 int ixyz, nxyz;
01824 int nvoxel;
01825
01826
01827
01828 nxyz = option_data->nxyz;
01829 nvoxel = option_data->nvoxel;
01830
01831
01832 ssab = (float *) malloc (sizeof(float)*nxyz);
01833 y = (float *) malloc (sizeof(float)*nxyz);
01834 if ((y == NULL) || (ssab == NULL))
01835 ANOVA_error ("unable to allocate sufficient memory");
01836
01837
01838
01839 volume_read ("sstr", ssab, nxyz);
01840
01841 volume_read ("ssa", y, nxyz);
01842 for (ixyz = 0; ixyz < nxyz; ixyz++)
01843 ssab[ixyz] -= y[ixyz];
01844
01845 volume_read ("ssb", y, nxyz);
01846 for (ixyz = 0; ixyz < nxyz; ixyz++)
01847 ssab[ixyz] -= y[ixyz];
01848
01849
01850
01851 for (ixyz = 0; ixyz < nxyz; ixyz++)
01852 if (ssab[ixyz] < 0.0) ssab[ixyz] = 0.0;
01853
01854
01855 if (nvoxel > 0)
01856 printf ("SSAB = %f \n", ssab[nvoxel-1]);
01857 volume_write ("ssab", ssab, nxyz);
01858
01859
01860 free (y); y = NULL;
01861 free (ssab); ssab = NULL;
01862
01863 }
01864
01865
01866
01867
01868
01869
01870
01871
01872
01873
01874
01875
01876
01877 void calculate_ftr (anova_options * option_data)
01878 {
01879 const float EPSILON = 1.0e-10;
01880 float * mstr = NULL;
01881 float * ftr = NULL;
01882 int a;
01883 int b;
01884 int n;
01885 int ixyz, nxyz;
01886 int nvoxel;
01887 float fval;
01888 float mse;
01889
01890
01891
01892
01893 a = option_data->a;
01894 b = option_data->b;
01895 n = option_data->n;
01896 nxyz = option_data->nxyz;
01897 nvoxel = option_data->nvoxel;
01898
01899
01900 ftr = (float *) malloc(sizeof(float)*nxyz);
01901 mstr = (float *) malloc(sizeof(float)*nxyz);
01902 if ((ftr == NULL) || (mstr == NULL))
01903 ANOVA_error ("unable to allocate sufficient memory");
01904
01905
01906 volume_read ("sstr", mstr, nxyz);
01907 fval = 1.0 / (a*b - 1.0);
01908 for (ixyz = 0; ixyz < nxyz; ixyz++)
01909 mstr[ixyz] = mstr[ixyz] * fval;
01910 if (nvoxel > 0)
01911 printf ("MSTR = %f \n", mstr[nvoxel-1]);
01912
01913
01914 volume_read ("sse", ftr, nxyz);
01915 fval = 1.0 / (a * b * (n-1));
01916 for (ixyz = 0; ixyz < nxyz; ixyz++)
01917 {
01918 mse = ftr[ixyz] * fval;
01919 if (mse < EPSILON)
01920 ftr[ixyz] = 0.0;
01921 else
01922 ftr[ixyz] = mstr[ixyz] / mse;
01923 }
01924 if (nvoxel > 0)
01925 printf ("FTR = %f \n", ftr[nvoxel-1]);
01926
01927
01928 for (ixyz = 0; ixyz < nxyz; ixyz++)
01929 mstr[ixyz] = sqrt(mstr[ixyz]);
01930 write_afni_data (option_data, option_data->ftrname,
01931 mstr, ftr, a*b-1, a*b*(n-1));
01932
01933
01934 volume_delete ("sstr");
01935
01936
01937 free (mstr); mstr = NULL;
01938 free (ftr); ftr = NULL;
01939
01940 }
01941
01942
01943
01944
01945
01946
01947
01948
01949
01950
01951
01952
01953
01954
01955
01956
01957
01958
01959
01960
01961
01962 void calculate_fa (anova_options * option_data)
01963 {
01964 const float EPSILON = 1.0e-10;
01965 float * msa = NULL;
01966 float * fa = NULL;
01967 int a;
01968 int b;
01969 int n;
01970 int ixyz, nxyz;
01971 int nvoxel;
01972 int numdf;
01973 int dendf;
01974 float mse;
01975 float msab;
01976
01977
01978
01979 a = option_data->a;
01980 b = option_data->b;
01981 n = option_data->n;
01982 nxyz = option_data->nxyz;
01983 nvoxel = option_data->nvoxel;
01984
01985
01986 fa = (float *) malloc(sizeof(float)*nxyz);
01987 msa = (float *) malloc(sizeof(float)*nxyz);
01988 if ((fa == NULL) || (msa == NULL))
01989 ANOVA_error ("unable to allocate sufficient memory");
01990
01991
01992 volume_read ("ssa", msa, nxyz);
01993 numdf = a - 1;
01994 for (ixyz = 0; ixyz < nxyz; ixyz++)
01995 msa[ixyz] = msa[ixyz] / numdf;
01996 if (nvoxel > 0)
01997 printf ("MSA = %f \n", msa[nvoxel-1]);
01998
01999
02000 if (option_data->model == 1)
02001 {
02002
02003 volume_read ("sse", fa, nxyz);
02004 dendf = a * b * (n-1);
02005 for (ixyz = 0; ixyz < nxyz; ixyz++)
02006 {
02007 mse = fa[ixyz] / dendf;
02008 if (mse < EPSILON)
02009 fa[ixyz] = 0.0;
02010 else
02011 fa[ixyz] = msa[ixyz] / mse;
02012 }
02013 }
02014 else
02015 {
02016
02017 volume_read ("ssab", fa, nxyz);
02018 dendf = (a-1) * (b-1);
02019 for (ixyz = 0; ixyz < nxyz; ixyz++)
02020 {
02021 msab = fa[ixyz] / dendf;
02022 if (msab < EPSILON)
02023 fa[ixyz] = 0.0;
02024 else
02025 fa[ixyz] = msa[ixyz] / msab;
02026 }
02027 }
02028
02029 if (nvoxel > 0)
02030 printf ("FA = %f \n", fa[nvoxel-1]);
02031
02032
02033 for (ixyz = 0; ixyz < nxyz; ixyz++)
02034 msa[ixyz] = sqrt(msa[ixyz]);
02035 write_afni_data (option_data, option_data->faname,
02036 msa, fa, numdf, dendf);
02037
02038
02039 volume_delete ("ssa");
02040
02041
02042 free (msa); msa = NULL;
02043 free (fa); fa = NULL;
02044
02045 }
02046
02047
02048
02049
02050
02051
02052
02053
02054
02055
02056
02057
02058
02059
02060
02061
02062
02063
02064
02065
02066
02067 void calculate_fb (anova_options * option_data)
02068 {
02069 const float EPSILON = 1.0e-10;
02070 float * msb = NULL;
02071 float * fb = NULL;
02072 int a;
02073 int b;
02074 int n;
02075 int ixyz, nxyz;
02076 int nvoxel;
02077 int numdf;
02078 int dendf;
02079 float mse;
02080 float msab;
02081
02082
02083
02084 a = option_data->a;
02085 b = option_data->b;
02086 n = option_data->n;
02087 nxyz = option_data->nxyz;
02088 nvoxel = option_data->nvoxel;
02089
02090
02091 fb = (float *) malloc(sizeof(float)*nxyz);
02092 msb = (float *) malloc(sizeof(float)*nxyz);
02093 if ((fb == NULL) || (msb == NULL))
02094 ANOVA_error ("unable to allocate sufficient memory");
02095
02096
02097 volume_read ("ssb", msb, nxyz);
02098 numdf = b - 1;
02099 for (ixyz = 0; ixyz < nxyz; ixyz++)
02100 msb[ixyz] = msb[ixyz] / numdf;
02101 if (nvoxel > 0)
02102 printf ("MSB = %f \n", msb[nvoxel-1]);
02103
02104
02105 if ((option_data->model == 1) || (option_data->model == 3))
02106 {
02107
02108 volume_read ("sse", fb, nxyz);
02109 dendf = a * b * (n-1);
02110 for (ixyz = 0; ixyz < nxyz; ixyz++)
02111 {
02112 mse = fb[ixyz] / dendf;
02113 if (mse < EPSILON)
02114 fb[ixyz] = 0.0;
02115 else
02116 fb[ixyz] = msb[ixyz] / mse;
02117 }
02118 }
02119 else
02120 {
02121
02122 volume_read ("ssab", fb, nxyz);
02123 dendf = (a-1) * (b-1);
02124 for (ixyz = 0; ixyz < nxyz; ixyz++)
02125 {
02126 msab = fb[ixyz] / dendf;
02127 if (msab < EPSILON)
02128 fb[ixyz] = 0.0;
02129 else
02130 fb[ixyz] = msb[ixyz] / msab;
02131 }
02132 }
02133
02134 if (nvoxel > 0)
02135 printf ("FB = %f \n", fb[nvoxel-1]);
02136
02137
02138 for (ixyz = 0; ixyz < nxyz; ixyz++)
02139 msb[ixyz] = sqrt(msb[ixyz]);
02140 write_afni_data (option_data, option_data->fbname,
02141 msb, fb, numdf, dendf);
02142
02143
02144 volume_delete ("ssb");
02145
02146
02147 free (msb); msb = NULL;
02148 free (fb); fb = NULL;
02149
02150 }
02151
02152
02153
02154
02155
02156
02157
02158
02159
02160
02161
02162
02163
02164 void calculate_fab (anova_options * option_data)
02165 {
02166 const float EPSILON = 1.0e-10;
02167 float * msab = NULL;
02168 float * fab = NULL;
02169 int a;
02170 int b;
02171 int n;
02172 int ixyz, nxyz;
02173 int nvoxel;
02174 float fval;
02175 float mse;
02176
02177
02178
02179 a = option_data->a;
02180 b = option_data->b;
02181 n = option_data->n;
02182 nxyz = option_data->nxyz;
02183 nvoxel = option_data->nvoxel;
02184
02185
02186 fab = (float *) malloc(sizeof(float)*nxyz);
02187 msab = (float *) malloc(sizeof(float)*nxyz);
02188 if ((fab == NULL) || (msab == NULL))
02189 ANOVA_error ("unable to allocate sufficient memory");
02190
02191
02192 volume_read ("ssab", msab, nxyz);
02193 fval = 1.0 / ((a - 1.0)*(b - 1.0));
02194 for (ixyz = 0; ixyz < nxyz; ixyz++)
02195 msab[ixyz] = msab[ixyz] * fval;
02196 if (nvoxel > 0)
02197 printf ("MSAB = %f \n", msab[nvoxel-1]);
02198
02199
02200 volume_read ("sse", fab, nxyz);
02201 fval = 1.0 / (a * b * (n-1));
02202 for (ixyz = 0; ixyz < nxyz; ixyz++)
02203 {
02204 mse = fab[ixyz] * fval;
02205 if (mse < EPSILON)
02206 fab[ixyz] = 0.0;
02207 else
02208 fab[ixyz] = msab[ixyz] / mse;
02209 }
02210 if (nvoxel > 0)
02211 printf ("FAB = %f \n", fab[nvoxel-1]);
02212
02213
02214 for (ixyz = 0; ixyz < nxyz; ixyz++)
02215 msab[ixyz] = sqrt(msab[ixyz]);
02216 write_afni_data (option_data, option_data->fabname,
02217 msab, fab, (a-1)*(b-1), a*b*(n-1));
02218
02219
02220 free (msab); msab = NULL;
02221 free (fab); fab = NULL;
02222
02223 }
02224
02225
02226
02227
02228
02229
02230
02231
02232
02233
02234
02235
02236
02237
02238
02239
02240
02241
02242
02243 void calculate_ameans (anova_options * option_data)
02244 {
02245 float * mean = NULL;
02246 float * tmean = NULL;
02247 int imean;
02248 int level;
02249 int n;
02250 int ixyz, nxyz;
02251 int nvoxel;
02252 int b;
02253 int num_means;
02254 int df;
02255
02256
02257
02258 b = option_data->b;
02259 n = option_data->n;
02260 num_means = option_data->num_ameans;
02261 nxyz = option_data->nxyz;
02262 nvoxel = option_data->nvoxel;
02263
02264
02265 df = b * n - 1;
02266
02267
02268 mean = (float *) malloc(sizeof(float)*nxyz);
02269 tmean = (float *) malloc(sizeof(float)*nxyz);
02270 if ((mean == NULL) || (tmean == NULL))
02271 ANOVA_error ("unable to allocate sufficient memory");
02272
02273
02274 for (imean = 0; imean < num_means; imean++)
02275 {
02276 level = option_data->ameans[imean];
02277
02278
02279 calculate_sum (option_data, level, -1, mean);
02280 calculate_sum_sq (option_data, level, -1, tmean);
02281
02282
02283 for (ixyz = 0; ixyz < nxyz; ixyz++)
02284 mean[ixyz] /= (df + 1.0);
02285
02286
02287 calculate_t_from_sums(tmean, mean, tmean, df, nxyz);
02288
02289 if (nvoxel > 0)
02290 printf ("factor A level %d: mean = %f, t = %f, df = %d\n",
02291 level+1, mean[nvoxel-1], tmean[nvoxel-1], df);
02292
02293 #if 0
02294
02295 calculate_sum (option_data, level, -1, mean);
02296 for (ixyz = 0; ixyz < nxyz; ixyz++)
02297 mean[ixyz] = mean[ixyz] / (n*b);
02298
02299
02300 if (option_data->model == 1)
02301 {
02302 volume_read ("sse", tmean, nxyz);
02303 df = a*b*(n-1);
02304 }
02305 else
02306 {
02307 volume_read ("ssab", tmean, nxyz);
02308 df = (a-1)*(b-1);
02309 }
02310 fval = (1.0 / df) * (1.0 / (b*n));
02311 for (ixyz = 0; ixyz < nxyz; ixyz++)
02312 {
02313 stddev = sqrt(tmean[ixyz] * fval);
02314 if (stddev < EPSILON) tmean[ixyz] = 0.0;
02315 else tmean[ixyz] = mean[ixyz] / stddev;
02316 }
02317 #endif
02318
02319
02320 write_afni_data (option_data, option_data->amname[imean],
02321 mean, tmean, df, 0);
02322
02323 }
02324
02325
02326 free (tmean); tmean = NULL;
02327 free (mean); mean = NULL;
02328 }
02329
02330
02331
02332
02333
02334
02335
02336
02337
02338
02339 void calculate_bmeans (anova_options * option_data)
02340 {
02341 const float EPSILON = 1.0e-10;
02342 float * mean = NULL;
02343 float * tmean = NULL;
02344 int imean;
02345 int level;
02346 int n;
02347 int ixyz, nxyz;
02348 int nvoxel;
02349 int a;
02350 int b;
02351 int nt;
02352 int num_means;
02353 int df;
02354 float fval;
02355 float stddev;
02356
02357
02358
02359 a = option_data->a;
02360 b = option_data->b;
02361 n = option_data->n;
02362 df = a * b * (n-1);
02363 nt = option_data->nt;
02364 num_means = option_data->num_bmeans;
02365 nxyz = option_data->nxyz;
02366 nvoxel = option_data->nvoxel;
02367
02368
02369 mean = (float *) malloc(sizeof(float)*nxyz);
02370 tmean = (float *) malloc(sizeof(float)*nxyz);
02371 if ((mean == NULL) || (tmean == NULL))
02372 ANOVA_error ("unable to allocate sufficient memory");
02373
02374
02375 for (imean = 0; imean < num_means; imean++)
02376 {
02377 level = option_data->bmeans[imean];
02378
02379
02380 calculate_sum (option_data, -1, level, mean);
02381 for (ixyz = 0; ixyz < nxyz; ixyz++)
02382 mean[ixyz] = mean[ixyz] / (n*a);
02383 if (nvoxel > 0)
02384 printf ("Mean of factor B level %d = %f \n", level+1, mean[nvoxel-1]);
02385
02386
02387 volume_read ("sse", tmean, nxyz);
02388 fval = (1.0 / df) * (1.0 / (a*n));
02389 for (ixyz = 0; ixyz < nxyz; ixyz++)
02390 {
02391 stddev = sqrt(tmean[ixyz] * fval);
02392 if (stddev < EPSILON)
02393 tmean[ixyz] = 0.0;
02394 else
02395 tmean[ixyz] = mean[ixyz] / stddev;
02396 }
02397 if (nvoxel > 0)
02398 printf ("t for mean of factor B level %d = %f \n",
02399 level+1, tmean[nvoxel-1]);
02400
02401
02402 write_afni_data (option_data, option_data->bmname[imean],
02403 mean, tmean, df, 0);
02404
02405 }
02406
02407
02408 free (tmean); tmean = NULL;
02409 free (mean); mean = NULL;
02410 }
02411
02412
02413
02414
02415
02416
02417
02418
02419
02420
02421 void calculate_xmeans (anova_options * option_data)
02422 {
02423 const float EPSILON = 1.0e-10;
02424 float * mean = NULL;
02425 float * tmean = NULL;
02426 int imean;
02427 int alevel;
02428 int blevel;
02429 int n;
02430 int ixyz, nxyz;
02431 int nvoxel;
02432 int a;
02433 int b;
02434 int nt;
02435 int num_means;
02436 int df;
02437 float fval;
02438 float stddev;
02439
02440
02441
02442 a = option_data->a;
02443 b = option_data->b;
02444 n = option_data->n;
02445 df = a * b * (n-1);
02446 nt = option_data->nt;
02447 num_means = option_data->num_xmeans;
02448 nxyz = option_data->nxyz;
02449 nvoxel = option_data->nvoxel;
02450
02451
02452 mean = (float *) malloc(sizeof(float)*nxyz);
02453 tmean = (float *) malloc(sizeof(float)*nxyz);
02454 if ((mean == NULL) || (tmean == NULL))
02455 ANOVA_error ("unable to allocate sufficient memory");
02456
02457
02458 for (imean = 0; imean < num_means; imean++)
02459 {
02460 alevel = option_data->xmeans[imean][0];
02461 blevel = option_data->xmeans[imean][1];
02462
02463
02464 calculate_sum (option_data, alevel, blevel, mean);
02465 for (ixyz = 0; ixyz < nxyz; ixyz++)
02466 mean[ixyz] = mean[ixyz] / n;
02467 if (nvoxel > 0)
02468 printf ("Mean of Cell[%d][%d] = %f \n",
02469 alevel+1, blevel+1, mean[nvoxel-1]);
02470
02471
02472 volume_read ("sse", tmean, nxyz);
02473 fval = (1.0 / df) * (1.0 / n);
02474 for (ixyz = 0; ixyz < nxyz; ixyz++)
02475 {
02476 stddev = sqrt(tmean[ixyz] * fval);
02477 if (stddev < EPSILON)
02478 tmean[ixyz] = 0.0;
02479 else
02480 tmean[ixyz] = mean[ixyz] / stddev;
02481 }
02482 if (nvoxel > 0)
02483 printf ("t-stat for Mean of Cell[%d][%d] = %f \n",
02484 alevel+1, blevel+1, tmean[nvoxel-1]);
02485
02486
02487 write_afni_data (option_data, option_data->xmname[imean],
02488 mean, tmean, df, 0);
02489
02490 }
02491
02492
02493 free (tmean); tmean = NULL;
02494 free (mean); mean = NULL;
02495 }
02496
02497
02498
02499
02500
02501
02502
02503
02504
02505
02506
02507
02508
02509 void calculate_adifferences (anova_options * option_data)
02510 {
02511 float * diff = NULL;
02512 float * tdiff = NULL;
02513 float * contrast;
02514 int a;
02515 int b;
02516 int nxyz;
02517 int nvoxel;
02518 int num_diffs;
02519 int idiff;
02520 int i, j;
02521 int n;
02522 int df, df_prod;
02523
02524
02525
02526 a = option_data->a;
02527 b = option_data->b;
02528 n = option_data->n;
02529 num_diffs = option_data->num_adiffs;
02530 nxyz = option_data->nxyz;
02531 nvoxel = option_data->nvoxel;
02532
02533
02534 df = b*n - 1;
02535 df_prod = df * (df+1);
02536
02537
02538 diff = (float *) malloc(sizeof(float)*nxyz);
02539 tdiff = (float *) malloc(sizeof(float)*nxyz);
02540 contrast = (float *) malloc(sizeof(float)*a);
02541 if ((diff == NULL) || (tdiff == NULL) || (contrast == NULL))
02542 ANOVA_error ("calc_adiffs: unable to allocate sufficient memory");
02543
02544
02545 for (idiff = 0; idiff < num_diffs; idiff++)
02546 {
02547 for (i = 0 ; i < a; i++ ) contrast[i] = 0.0;
02548
02549 i = option_data->adiffs[idiff][0];
02550 j = option_data->adiffs[idiff][1];
02551
02552
02553 contrast[i] = 1; contrast[j] = -1;
02554
02555
02556 calc_acontr_mean(option_data, contrast, diff);
02557 calc_sum_sq_acontr(option_data, contrast, tdiff);
02558 calculate_t_from_sums(tdiff, diff, tdiff, df, nxyz);
02559
02560 if (nvoxel > 0)
02561 printf ("Difference of factor A level %d - level %d = %f \n",
02562 i+1, j+1, diff[nvoxel-1]);
02563
02564 #if 0
02565
02566 calculate_sum (option_data, i, -1, diff);
02567 for (ixyz = 0; ixyz < nxyz; ixyz++) diff[ixyz] = diff[ixyz] / (b*n);
02568
02569
02570 calculate_sum (option_data, j, -1, tdiff);
02571 for (ixyz = 0; ixyz < nxyz; ixyz++) diff[ixyz] -= tdiff[ixyz] / (b*n);
02572
02573
02574 if (option_data->model == 1)
02575 {
02576 volume_read ("sse", tdiff, nxyz);
02577 df = a*b*(n-1);
02578 } else {
02579 volume_read ("ssab", tdiff, nxyz);
02580 df = (a-1)*(b-1);
02581 }
02582 fval = (1.0 / df) * (2.0 / (b*n));
02583 for (ixyz = 0; ixyz < nxyz; ixyz++)
02584 {
02585 stddev = sqrt (tdiff[ixyz] * fval);
02586 if (stddev < EPSILON) tdiff[ixyz] = 0.0;
02587 else tdiff[ixyz] = diff[ixyz] / stddev;
02588 }
02589 #endif
02590
02591 if (nvoxel > 0)
02592 printf ("t for difference of factor A level %d - level %d = %f \n",
02593 i+1, j+1, tdiff[nvoxel-1]);
02594
02595
02596 write_afni_data (option_data, option_data->adname[idiff],
02597 diff, tdiff, df, 0);
02598
02599 }
02600
02601
02602 free (tdiff); tdiff = NULL;
02603 free (diff); diff = NULL;
02604 free (contrast); contrast = NULL;
02605 }
02606
02607
02608
02609
02610
02611
02612
02613
02614
02615
02616 void calculate_bdifferences (anova_options * option_data)
02617 {
02618 const float EPSILON = 1.0e-10;
02619 float * diff = NULL;
02620 float * tdiff = NULL;
02621 int a;
02622 int b;
02623 int ixyz, nxyz;
02624 int nvoxel;
02625 int num_diffs;
02626 int idiff;
02627 int i, j;
02628 int n;
02629 int df;
02630 float fval;
02631 float stddev;
02632
02633
02634
02635 a = option_data->a;
02636 b = option_data->b;
02637 n = option_data->n;
02638 df = a*b*(n-1);
02639 num_diffs = option_data->num_bdiffs;
02640 nxyz = option_data->nxyz;
02641 nvoxel = option_data->nvoxel;
02642
02643
02644 diff = (float *) malloc(sizeof(float)*nxyz);
02645 tdiff = (float *) malloc(sizeof(float)*nxyz);
02646 if ((diff == NULL) || (tdiff == NULL))
02647 ANOVA_error ("unable to allocate sufficient memory");
02648
02649
02650 for (idiff = 0; idiff < num_diffs; idiff++)
02651 {
02652
02653
02654 i = option_data->bdiffs[idiff][0];
02655 calculate_sum (option_data, -1, i, diff);
02656 for (ixyz = 0; ixyz < nxyz; ixyz++)
02657 diff[ixyz] = diff[ixyz] / (a*n);
02658
02659
02660 j = option_data->bdiffs[idiff][1];
02661 calculate_sum (option_data, -1, j, tdiff);
02662 for (ixyz = 0; ixyz < nxyz; ixyz++)
02663 diff[ixyz] -= tdiff[ixyz] / (a*n);
02664 if (nvoxel > 0)
02665 printf ("Difference of factor B level %d - level %d = %f \n",
02666 i+1, j+1, diff[nvoxel-1]);
02667
02668
02669 volume_read ("sse", tdiff, nxyz);
02670 fval = (1.0 / df) * (2.0 / (a*n));
02671 for (ixyz = 0; ixyz < nxyz; ixyz++)
02672 {
02673 stddev = sqrt (tdiff[ixyz] * fval);
02674 if (stddev < EPSILON)
02675 tdiff[ixyz] = 0.0;
02676 else
02677 tdiff[ixyz] = diff[ixyz] / stddev;
02678 }
02679
02680 if (nvoxel > 0)
02681 printf ("t for difference of factor B level %d - level %d = %f \n",
02682 i+1, j+1, tdiff[nvoxel-1]);
02683
02684
02685 write_afni_data (option_data, option_data->bdname[idiff],
02686 diff, tdiff, df, 0);
02687
02688 }
02689
02690
02691 free (tdiff); tdiff = NULL;
02692 free (diff); diff = NULL;
02693 }
02694
02695
02696
02697
02698
02699
02700
02701
02702
02703
02704 void calculate_xdifferences (anova_options * option_data)
02705 {
02706 const float EPSILON = 1.0e-10;
02707 float * diff = NULL;
02708 float * tdiff = NULL;
02709 int a;
02710 int b;
02711 int ixyz, nxyz;
02712 int nvoxel;
02713 int num_diffs;
02714 int idiff;
02715 int ia, ib, ja, jb;
02716 int n;
02717 int df;
02718 float fval;
02719 float stddev;
02720
02721
02722
02723 a = option_data->a;
02724 b = option_data->b;
02725 n = option_data->n;
02726 df = a*b*(n-1);
02727 num_diffs = option_data->num_xdiffs;
02728 nxyz = option_data->nxyz;
02729 nvoxel = option_data->nvoxel;
02730
02731
02732 diff = (float *) malloc(sizeof(float)*nxyz);
02733 tdiff = (float *) malloc(sizeof(float)*nxyz);
02734 if ((diff == NULL) || (tdiff == NULL))
02735 ANOVA_error ("unable to allocate sufficient memory");
02736
02737
02738 for (idiff = 0; idiff < num_diffs; idiff++)
02739 {
02740
02741
02742 ia = option_data->xdiffs[idiff][0][0];
02743 ib = option_data->xdiffs[idiff][0][1];
02744 calculate_sum (option_data, ia, ib, diff);
02745 for (ixyz = 0; ixyz < nxyz; ixyz++)
02746 diff[ixyz] = diff[ixyz] / n;
02747
02748
02749 ja = option_data->xdiffs[idiff][1][0];
02750 jb = option_data->xdiffs[idiff][1][1];
02751 calculate_sum (option_data, ja, jb, tdiff);
02752 for (ixyz = 0; ixyz < nxyz; ixyz++)
02753 diff[ixyz] -= tdiff[ixyz] / n;
02754 if (nvoxel > 0)
02755 printf ("Difference of Cell[%d][%d] - Cell[%d][%d] = %f \n",
02756 ia+1, ib+1, ja+1, jb+1, diff[nvoxel-1]);
02757
02758
02759 volume_read ("sse", tdiff, nxyz);
02760 fval = (1.0 / df) * (2.0 / n);
02761 for (ixyz = 0; ixyz < nxyz; ixyz++)
02762 {
02763 stddev = sqrt (tdiff[ixyz] * fval);
02764 if (stddev < EPSILON)
02765 tdiff[ixyz] = 0.0;
02766 else
02767 tdiff[ixyz] = diff[ixyz] / stddev;
02768 }
02769
02770 if (nvoxel > 0)
02771 printf ("t-stat for Cell[%d][%d] - Cell[%d][%d] = %f \n",
02772 ia+1, ib+1, ja+1, jb+1, tdiff[nvoxel-1]);
02773
02774
02775 write_afni_data (option_data, option_data->xdname[idiff],
02776 diff, tdiff, df, 0);
02777
02778 }
02779
02780
02781 free (tdiff); tdiff = NULL;
02782 free (diff); diff = NULL;
02783
02784 }
02785
02786
02787
02788
02789
02790
02791
02792
02793
02794
02795
02796
02797 void calculate_acontrasts (anova_options * option_data)
02798 {
02799 float * contr = NULL;
02800 float * tcontr = NULL;
02801 int b;
02802 int nxyz;
02803 int nvoxel;
02804 int num_contr;
02805 int icontr;
02806 int n;
02807 int df;
02808
02809
02810
02811 b = option_data->b;
02812 n = option_data->n;
02813 num_contr = option_data->num_acontr;
02814 nxyz = option_data->nxyz;
02815 nvoxel = option_data->nvoxel;
02816
02817 df = b*n - 1;
02818
02819
02820 contr = (float *) malloc(sizeof(float)*nxyz);
02821 tcontr = (float *) malloc(sizeof(float)*nxyz);
02822 if ((contr == NULL) || (tcontr == NULL))
02823 ANOVA_error ("unable to allocate sufficient memory");
02824
02825
02826
02827 for (icontr = 0; icontr < num_contr; icontr++)
02828 {
02829
02830 calc_acontr_mean(option_data, option_data->acontr[icontr], contr);
02831 calc_sum_sq_acontr(option_data, option_data->acontr[icontr], tcontr);
02832
02833
02834 calculate_t_from_sums(tcontr, contr, tcontr, df, nxyz);
02835
02836 #if 0
02837 for (level = 0; level < a; level++)
02838 {
02839 c = option_data->acontr[icontr][level];
02840 if (c == 0.0) continue;
02841
02842
02843 calculate_sum (option_data, level, -1, tcontr);
02844 fval += c * c / (b*n);
02845 for (ixyz = 0; ixyz < nxyz; ixyz++)
02846 contr[ixyz] += c * tcontr[ixyz] / (b*n);
02847 }
02848
02849
02850 if (option_data->model == 1) {
02851 volume_read ("sse", tcontr, nxyz);
02852 df = a*b*(n-1);
02853 } else {
02854 volume_read ("ssab", tcontr, nxyz);
02855 df = (a-1)*(b-1);
02856 }
02857
02858
02859 for (ixyz = 0; ixyz < nxyz; ixyz++)
02860 {
02861 stddev = sqrt ((tcontr[ixyz] / df) * fval);
02862 if (stddev < EPSILON) tcontr[ixyz] = 0.0;
02863 else tcontr[ixyz] = contr[ixyz] / stddev;
02864 }
02865 #endif
02866
02867 if (nvoxel > 0)
02868 printf ("No.%d contrast for factor A = %f, t = %f, df = %d\n",
02869 icontr+1, contr[nvoxel-1], tcontr[nvoxel-1], df);
02870
02871
02872 write_afni_data (option_data, option_data->acname[icontr],
02873 contr, tcontr, df, 0);
02874
02875 }
02876
02877
02878 free (tcontr); tcontr = NULL;
02879 free (contr); contr = NULL;
02880
02881 }
02882
02883
02884
02885
02886
02887
02888
02889
02890
02891 void calculate_bcontrasts (anova_options * option_data)
02892 {
02893 const float EPSILON = 1.0e-10;
02894 float * contr = NULL;
02895 float * tcontr = NULL;
02896 int a;
02897 int b;
02898 int ixyz, nxyz;
02899 int nvoxel;
02900 int num_contr;
02901 int icontr;
02902 int level;
02903 int df;
02904 int n;
02905 float fval;
02906 float c;
02907 float stddev;
02908
02909
02910
02911 a = option_data->a;
02912 b = option_data->b;
02913 n = option_data->n;
02914 num_contr = option_data->num_bcontr;
02915 nxyz = option_data->nxyz;
02916 nvoxel = option_data->nvoxel;
02917
02918
02919 contr = (float *) malloc(sizeof(float)*nxyz);
02920 tcontr = (float *) malloc(sizeof(float)*nxyz);
02921 if ((contr == NULL) || (tcontr == NULL))
02922 ANOVA_error ("unable to allocate sufficient memory");
02923
02924
02925
02926 for (icontr = 0; icontr < num_contr; icontr++)
02927 {
02928 volume_zero (contr, nxyz);
02929 fval = 0.0;
02930
02931 for (level = 0; level < b; level++)
02932 {
02933 c = option_data->bcontr[icontr][level];
02934 if (c == 0.0) continue;
02935
02936
02937 calculate_sum (option_data, -1, level, tcontr);
02938 fval += c * c / (a*n);
02939 for (ixyz = 0; ixyz < nxyz; ixyz++)
02940 contr[ixyz] += c * tcontr[ixyz] / (a*n);
02941 }
02942 if (nvoxel > 0)
02943 printf ("No.%d contrast for factor B = %f \n",
02944 icontr+1, contr[nvoxel-1]);
02945
02946
02947 volume_read ("sse", tcontr, nxyz);
02948 df = a * b * (n-1);
02949 for (ixyz = 0; ixyz < nxyz; ixyz++)
02950 {
02951 stddev = sqrt ((tcontr[ixyz] / df) * fval);
02952 if (stddev < EPSILON)
02953 tcontr[ixyz] = 0.0;
02954 else
02955 tcontr[ixyz] = contr[ixyz] / stddev;
02956 }
02957
02958 if (nvoxel > 0)
02959 printf ("t of No.%d contrast for factor B = %f \n",
02960 icontr+1, tcontr[nvoxel-1]);
02961
02962
02963 write_afni_data (option_data, option_data->bcname[icontr],
02964 contr, tcontr, a*b*(n-1), 0);
02965
02966 }
02967
02968
02969 free (tcontr); tcontr = NULL;
02970 free (contr); contr = NULL;
02971 }
02972
02973
02974
02975
02976
02977
02978
02979
02980
02981
02982 void calculate_xcontrasts (anova_options * option_data)
02983 {
02984 const float EPSILON = 1.0e-10;
02985 float * contr = NULL;
02986 float * tcontr = NULL;
02987 int a;
02988 int b;
02989 int ixyz, nxyz;
02990 int nvoxel;
02991 int num_contr;
02992 int icontr;
02993 int ilevel, jlevel;
02994 int df;
02995 int n;
02996 float fval;
02997 float c;
02998 float stddev;
02999
03000
03001
03002 a = option_data->a;
03003 b = option_data->b;
03004 n = option_data->n;
03005 num_contr = option_data->num_xcontr;
03006 nxyz = option_data->nxyz;
03007 nvoxel = option_data->nvoxel;
03008
03009
03010 contr = (float *) malloc(sizeof(float)*nxyz);
03011 tcontr = (float *) malloc(sizeof(float)*nxyz);
03012 if ((contr == NULL) || (tcontr == NULL))
03013 ANOVA_error ("unable to allocate sufficient memory");
03014
03015
03016
03017 for (icontr = 0; icontr < num_contr; icontr++)
03018 {
03019 volume_zero (contr, nxyz);
03020 fval = 0.0;
03021
03022 for (ilevel = 0; ilevel < a; ilevel++)
03023 for (jlevel = 0; jlevel < b; jlevel++)
03024 {
03025 c = option_data->xcontr[icontr][ilevel][jlevel];
03026 if (c == 0.0) continue;
03027
03028
03029 calculate_sum (option_data, ilevel, jlevel, tcontr);
03030 fval += c * c / n;
03031 for (ixyz = 0; ixyz < nxyz; ixyz++)
03032 contr[ixyz] += c * tcontr[ixyz] / n;
03033 }
03034 if (nvoxel > 0)
03035 printf ("No.%d contrast for cell means = %f \n",
03036 icontr+1, contr[nvoxel-1]);
03037
03038
03039 volume_read ("sse", tcontr, nxyz);
03040 df = a * b * (n-1);
03041 for (ixyz = 0; ixyz < nxyz; ixyz++)
03042 {
03043 stddev = sqrt ((tcontr[ixyz] / df) * fval);
03044 if (stddev < EPSILON)
03045 tcontr[ixyz] = 0.0;
03046 else
03047 tcontr[ixyz] = contr[ixyz] / stddev;
03048 }
03049
03050 if (nvoxel > 0)
03051 printf ("t-stat for No.%d contrast for cell means = %f \n",
03052 icontr+1, tcontr[nvoxel-1]);
03053
03054
03055 write_afni_data (option_data, option_data->xcname[icontr],
03056 contr, tcontr, a*b*(n-1), 0);
03057
03058 }
03059
03060
03061 free (tcontr); tcontr = NULL;
03062 free (contr); contr = NULL;
03063
03064 }
03065
03066
03067
03068
03069
03070
03071
03072 void calculate_anova (anova_options * option_data)
03073 {
03074
03075
03076 calculate_ss0 (option_data);
03077 calculate_ssi (option_data);
03078 calculate_ssj (option_data);
03079 calculate_ssij (option_data);
03080 if (option_data->n != 1) calculate_ssijk (option_data);
03081
03082
03083
03084 if (option_data->n != 1)
03085 {
03086 calculate_sse (option_data);
03087 volume_delete ("ssijk");
03088 }
03089
03090
03091 calculate_sstr (option_data);
03092 volume_delete ("ssij");
03093
03094
03095 calculate_ssa (option_data);
03096 volume_delete ("ssi");
03097
03098
03099 calculate_ssb (option_data);
03100 volume_delete ("ssj");
03101
03102 volume_delete ("ss0");
03103
03104
03105 calculate_ssab (option_data);
03106
03107 }
03108
03109
03110
03111
03112
03113
03114
03115 void analyze_results (anova_options * option_data)
03116 {
03117
03118
03119 if (option_data->nftr) calculate_ftr (option_data);
03120
03121
03122 if (option_data->nfa) calculate_fa (option_data);
03123
03124
03125 if (option_data->nfb) calculate_fb (option_data);
03126
03127
03128 if (option_data->nfab) calculate_fab (option_data);
03129
03130
03131 if (option_data->num_ameans) calculate_ameans (option_data);
03132
03133
03134 if (option_data->num_bmeans) calculate_bmeans (option_data);
03135
03136
03137 if (option_data->num_xmeans) calculate_xmeans (option_data);
03138
03139
03140 if (option_data->num_adiffs) calculate_adifferences (option_data);
03141
03142
03143 if (option_data->num_bdiffs) calculate_bdifferences (option_data);
03144
03145
03146 if (option_data->num_xdiffs) calculate_xdifferences (option_data);
03147
03148
03149 if (option_data->num_acontr) calculate_acontrasts (option_data);
03150
03151
03152 if (option_data->num_bcontr) calculate_bcontrasts (option_data);
03153
03154
03155 if (option_data->num_xcontr) calculate_xcontrasts (option_data);
03156
03157 }
03158
03159
03160
03161
03162
03163
03164
03165 void create_bucket (anova_options * option_data)
03166
03167 {
03168 char bucket_str[10000];
03169 char refit_str[10000];
03170 THD_3dim_dataset * dset=NULL;
03171 THD_3dim_dataset * new_dset=NULL;
03172 int i;
03173 int ibrick;
03174 char str[100];
03175
03176
03177
03178 dset = THD_open_dataset (option_data->first_dataset) ;
03179 if( ! ISVALID_3DIM_DATASET(dset) ){
03180 fprintf(stderr,"*** Unable to open dataset file %s\n",
03181 option_data->first_dataset);
03182 exit(1) ;
03183 }
03184
03185 if( DSET_IS_1D(dset) ) USE_1D_filenames(1) ;
03186 else if( DSET_IS_3D(dset) ) USE_1D_filenames(3) ;
03187
03188
03189
03190 new_dset = EDIT_empty_copy( dset ) ;
03191 THD_delete_3dim_dataset (dset , False); dset = NULL;
03192 EDIT_dset_items (new_dset,
03193 ADN_directory_name, option_data->session,
03194 ADN_none);
03195
03196
03197
03198 strcpy (bucket_str, "3dbucket");
03199 strcat (bucket_str, " -prefix ");
03200 strcat (bucket_str, option_data->bucket_filename);
03201
03202
03203
03204 strcpy (refit_str, "3drefit ");
03205 ibrick = -1;
03206
03207
03208
03209 if (option_data->nftr != 0)
03210 {
03211 add_file_name (new_dset, option_data->ftrname, bucket_str);
03212
03213 ibrick++;
03214 sprintf (str, " -sublabel %d %s:Inten ",
03215 ibrick, option_data->ftrname);
03216 strcat (refit_str, str);
03217
03218 ibrick++;
03219 sprintf (str, " -sublabel %d %s:F-stat ",
03220 ibrick, option_data->ftrname);
03221 strcat (refit_str, str);
03222 }
03223
03224
03225
03226 if (option_data->nfa != 0)
03227 {
03228 add_file_name (new_dset, option_data->faname, bucket_str);
03229
03230 ibrick++;
03231 sprintf (str, " -sublabel %d %s:Inten ",
03232 ibrick, option_data->faname);
03233 strcat (refit_str, str);
03234
03235 ibrick++;
03236 sprintf (str, " -sublabel %d %s:F-stat ",
03237 ibrick, option_data->faname);
03238 strcat (refit_str, str);
03239 }
03240
03241
03242
03243 if (option_data->nfb != 0)
03244 {
03245 add_file_name (new_dset, option_data->fbname, bucket_str);
03246
03247 ibrick++;
03248 sprintf (str, " -sublabel %d %s:Inten ",
03249 ibrick, option_data->fbname);
03250 strcat (refit_str, str);
03251
03252 ibrick++;
03253 sprintf (str, " -sublabel %d %s:F-stat ",
03254 ibrick, option_data->fbname);
03255 strcat (refit_str, str);
03256 }
03257
03258
03259
03260 if (option_data->nfab != 0)
03261 {
03262 add_file_name (new_dset, option_data->fabname, bucket_str);
03263
03264 ibrick++;
03265 sprintf (str, " -sublabel %d %s:Inten ",
03266 ibrick, option_data->fabname);
03267 strcat (refit_str, str);
03268
03269 ibrick++;
03270 sprintf (str, " -sublabel %d %s:F-stat ",
03271 ibrick, option_data->fabname);
03272 strcat (refit_str, str);
03273 }
03274
03275
03276
03277 if (option_data->num_ameans > 0)
03278 for (i = 0; i < option_data->num_ameans; i++)
03279 {
03280 add_file_name (new_dset, option_data->amname[i], bucket_str);
03281
03282 ibrick++;
03283 sprintf (str, " -sublabel %d %s:Mean ",
03284 ibrick, option_data->amname[i]);
03285 strcat (refit_str, str);
03286
03287 ibrick++;
03288 sprintf (str, " -sublabel %d %s:t-stat ",
03289 ibrick, option_data->amname[i]);
03290 strcat (refit_str, str);
03291 }
03292
03293
03294
03295 if (option_data->num_bmeans > 0)
03296 for (i = 0; i < option_data->num_bmeans; i++)
03297 {
03298 add_file_name (new_dset, option_data->bmname[i], bucket_str);
03299
03300 ibrick++;
03301 sprintf (str, " -sublabel %d %s:Mean ",
03302 ibrick, option_data->bmname[i]);
03303 strcat (refit_str, str);
03304
03305 ibrick++;
03306 sprintf (str, " -sublabel %d %s:t-stat ",
03307 ibrick, option_data->bmname[i]);
03308 strcat (refit_str, str);
03309 }
03310
03311
03312
03313 if (option_data->num_xmeans > 0)
03314 for (i = 0; i < option_data->num_xmeans; i++)
03315 {
03316 add_file_name (new_dset, option_data->xmname[i], bucket_str);
03317
03318 ibrick++;
03319 sprintf (str, " -sublabel %d %s:Mean ",
03320 ibrick, option_data->xmname[i]);
03321 strcat (refit_str, str);
03322
03323 ibrick++;
03324 sprintf (str, " -sublabel %d %s:t-stat ",
03325 ibrick, option_data->xmname[i]);
03326 strcat (refit_str, str);
03327 }
03328
03329
03330
03331 if (option_data->num_adiffs > 0)
03332 for (i = 0; i < option_data->num_adiffs; i++)
03333 {
03334 add_file_name (new_dset, option_data->adname[i], bucket_str);
03335
03336 ibrick++;
03337 sprintf (str, " -sublabel %d %s:Diff ",
03338 ibrick, option_data->adname[i]);
03339 strcat (refit_str, str);
03340
03341 ibrick++;
03342 sprintf (str, " -sublabel %d %s:t-stat ",
03343 ibrick, option_data->adname[i]);
03344 strcat (refit_str, str);
03345 }
03346
03347
03348
03349 if (option_data->num_bdiffs > 0)
03350 for (i = 0; i < option_data->num_bdiffs; i++)
03351 {
03352 add_file_name (new_dset, option_data->bdname[i], bucket_str);
03353
03354 ibrick++;
03355 sprintf (str, " -sublabel %d %s:Diff ",
03356 ibrick, option_data->bdname[i]);
03357 strcat (refit_str, str);
03358
03359 ibrick++;
03360 sprintf (str, " -sublabel %d %s:t-stat ",
03361 ibrick, option_data->bdname[i]);
03362 strcat (refit_str, str);
03363 }
03364
03365
03366
03367 if (option_data->num_xdiffs > 0)
03368 for (i = 0; i < option_data->num_xdiffs; i++)
03369 {
03370 add_file_name (new_dset, option_data->xdname[i], bucket_str);
03371
03372 ibrick++;
03373 sprintf (str, " -sublabel %d %s:Diff ",
03374 ibrick, option_data->xdname[i]);
03375 strcat (refit_str, str);
03376
03377 ibrick++;
03378 sprintf (str, " -sublabel %d %s:t-stat ",
03379 ibrick, option_data->xdname[i]);
03380 strcat (refit_str, str);
03381 }
03382
03383
03384
03385 if (option_data->num_acontr > 0)
03386 for (i = 0; i < option_data->num_acontr; i++)
03387 {
03388 add_file_name (new_dset, option_data->acname[i], bucket_str);
03389
03390 ibrick++;
03391 sprintf (str, " -sublabel %d %s:Contr ",
03392 ibrick, option_data->acname[i]);
03393 strcat (refit_str, str);
03394
03395 ibrick++;
03396 sprintf (str, " -sublabel %d %s:t-stat ",
03397 ibrick, option_data->acname[i]);
03398 strcat (refit_str, str);
03399 }
03400
03401
03402
03403 if (option_data->num_bcontr > 0)
03404 for (i = 0; i < option_data->num_bcontr; i++)
03405 {
03406 add_file_name (new_dset, option_data->bcname[i], bucket_str);
03407
03408 ibrick++;
03409 sprintf (str, " -sublabel %d %s:Contr ",
03410 ibrick, option_data->bcname[i]);
03411 strcat (refit_str, str);
03412
03413 ibrick++;
03414 sprintf (str, " -sublabel %d %s:t-stat ",
03415 ibrick, option_data->bcname[i]);
03416 strcat (refit_str, str);
03417 }
03418
03419
03420
03421 if (option_data->num_xcontr > 0)
03422 for (i = 0; i < option_data->num_xcontr; i++)
03423 {
03424 add_file_name (new_dset, option_data->xcname[i], bucket_str);
03425
03426 ibrick++;
03427 sprintf (str, " -sublabel %d %s:Contr ",
03428 ibrick, option_data->xcname[i]);
03429 strcat (refit_str, str);
03430
03431 ibrick++;
03432 sprintf (str, " -sublabel %d %s:t-stat ",
03433 ibrick, option_data->xcname[i]);
03434 strcat (refit_str, str);
03435 }
03436
03437
03438
03439
03440 printf("Writing `bucket' dataset ");
03441 printf("into %s\n", option_data->bucket_filename);
03442
03443 fprintf(stderr,"RUNNING COMMAND: %s\n",bucket_str) ;
03444 system (bucket_str);
03445
03446
03447
03448 add_file_name (new_dset, option_data->bucket_filename, refit_str);
03449 fprintf(stderr,"RUNNING COMMAND: %s\n",refit_str) ;
03450 system (refit_str);
03451
03452
03453
03454 THD_delete_3dim_dataset (new_dset , False); new_dset = NULL;
03455
03456 }
03457
03458
03459
03460
03461
03462
03463
03464
03465
03466 void terminate (anova_options * option_data)
03467 {
03468 int i, j;
03469 THD_3dim_dataset * dset=NULL;
03470 THD_3dim_dataset * new_dset=NULL;
03471 char filename[MAX_NAME_LENGTH];
03472
03473
03474
03475 volume_delete ("sstr");
03476 volume_delete ("sse");
03477 volume_delete ("ssa");
03478 volume_delete ("ssb");
03479 volume_delete ("ssab");
03480 for (i = 0; i < option_data->a; i++)
03481 {
03482 sprintf (filename, "ya.%d", i);
03483 volume_delete (filename);
03484 }
03485 for (j = 0; j < option_data->b; j++)
03486 {
03487 sprintf (filename, "yb.%d", j);
03488 volume_delete (filename);
03489 }
03490
03491
03492
03493 if (option_data->bucket_filename != NULL)
03494 create_bucket (option_data);
03495
03496
03497
03498
03499 if (option_data->bucket_filename != NULL)
03500 {
03501
03502
03503 dset = THD_open_dataset (option_data->first_dataset) ;
03504 if( ! ISVALID_3DIM_DATASET(dset) ){
03505 fprintf(stderr,"*** Unable to open dataset file %s\n",
03506 option_data->first_dataset);
03507 exit(1) ;
03508 }
03509
03510
03511 new_dset = EDIT_empty_copy (dset);
03512 THD_delete_3dim_dataset (dset , False); dset = NULL;
03513 EDIT_dset_items (new_dset,
03514 ADN_directory_name, option_data->session,
03515 ADN_none);
03516
03517
03518 if (option_data->nftr != 0)
03519 remove_dataset (new_dset, option_data->ftrname);
03520
03521
03522 if (option_data->nfa != 0)
03523 remove_dataset (new_dset, option_data->faname);
03524
03525
03526 if (option_data->nfb != 0)
03527 remove_dataset (new_dset, option_data->fbname);
03528
03529
03530 if (option_data->nfab != 0)
03531 remove_dataset (new_dset, option_data->fabname);
03532
03533
03534 if (option_data->num_ameans > 0)
03535 for (i = 0; i < option_data->num_ameans; i++)
03536 remove_dataset (new_dset, option_data->amname[i]);
03537
03538
03539 if (option_data->num_bmeans > 0)
03540 for (i = 0; i < option_data->num_bmeans; i++)
03541 remove_dataset (new_dset, option_data->bmname[i]);
03542
03543
03544 if (option_data->num_xmeans > 0)
03545 for (i = 0; i < option_data->num_xmeans; i++)
03546 remove_dataset (new_dset, option_data->xmname[i]);
03547
03548
03549 if (option_data->num_adiffs > 0)
03550 for (i = 0; i < option_data->num_adiffs; i++)
03551 remove_dataset (new_dset, option_data->adname[i]);
03552
03553
03554 if (option_data->num_bdiffs > 0)
03555 for (i = 0; i < option_data->num_bdiffs; i++)
03556 remove_dataset (new_dset, option_data->bdname[i]);
03557
03558
03559 if (option_data->num_xdiffs > 0)
03560 for (i = 0; i < option_data->num_xdiffs; i++)
03561 remove_dataset (new_dset, option_data->xdname[i]);
03562
03563
03564 if (option_data->num_acontr > 0)
03565 for (i = 0; i < option_data->num_acontr; i++)
03566 remove_dataset (new_dset, option_data->acname[i]);
03567
03568
03569 if (option_data->num_bcontr > 0)
03570 for (i = 0; i < option_data->num_bcontr; i++)
03571 remove_dataset (new_dset, option_data->bcname[i]);
03572
03573
03574 if (option_data->num_xcontr > 0)
03575 for (i = 0; i < option_data->num_xcontr; i++)
03576 remove_dataset (new_dset, option_data->xcname[i]);
03577
03578 THD_delete_3dim_dataset (new_dset , False); new_dset = NULL;
03579 }
03580
03581
03582
03583 destroy_anova_options (option_data);
03584
03585 }
03586
03587
03588
03589
03590
03591
03592
03593 int main (int argc, char ** argv)
03594 {
03595 anova_options * option_data = NULL;
03596
03597
03598
03599 printf ("\n\n");
03600 printf ("Program: %s \n", PROGRAM_NAME);
03601 printf ("Author: %s \n", PROGRAM_AUTHOR);
03602 printf ("Initial Release: %s \n", PROGRAM_INITIAL);
03603 printf ("Latest Revision: %s \n", PROGRAM_LATEST);
03604 printf ("\n");
03605
03606
03607
03608 machdep() ;
03609 { int new_argc ; char ** new_argv ;
03610 addto_args( argc , argv , &new_argc , &new_argv ) ;
03611 if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
03612 }
03613
03614
03615 initialize (argc, argv, &option_data);
03616
03617
03618 calculate_anova (option_data);
03619
03620
03621 analyze_results (option_data);
03622
03623
03624 terminate (option_data);
03625 free (option_data); option_data = NULL;
03626
03627 exit(0);
03628 }
