1 2 /***************************************************************************** 3 Major portions of this software are copyrighted by the Medical College 4 of Wisconsin, 1994-2000, and are released under the Gnu General Public 5 License, Version 2. See the file README.Copyright for details. 6 ******************************************************************************/ 7 8 /* 9 This is the header file for the 3dANOVA.lib library routines. 10 11 File: 3dANOVA.h 12 Author: B. D. Ward 13 Date: 20 January 1997 14 15 Mod: Added command -diskspace to print out how much disk space is 16 required to execute the problem. 17 Date: 23 January 1997 18 19 Mod: Change to routine write_afni_data to reduce truncation error. 20 Date: 27 January 1997 21 22 Mod: Added machine specific code for checking disk space. 23 Date: 29 January 1997 24 25 Mod: Extensive changes required to implement the 'bucket' dataset. 26 Date: 30 December 1997 27 28 Mod: Library routines moved to 3dANOVA.lib. 29 Date: 5 January 1998 30 31 Mod: Added software for statistical tests of individual cell means, 32 cell differences, and cell contrasts. 33 Date: 27 October 1998 34 35 Mod: Added changes for incorporating History notes. 36 Date: 09 September 1999 37 38 Mod: Replaced dataset input code with calls to THD_open_dataset, 39 to allow operator selection of individual sub-bricks for input. 40 Date: 02 December 1999 41 42 Mod: Increased maximum number of user requested means, differences, 43 and contrasts. 44 Date: 31 January 2000 45 46 Mod: Set MAX_NAME_LENGTH equal to THD_MAX_NAME. 47 Date: 02 December 2002 48 49 Mod: Added aBcontr/Abcontr fields to the anova_options struct. 50 Date: 14 October 2005 [rickr] 51 52 Mod: Added old_method flag to anova_options struct. 53 Added prototype for contrasts_are_valid(), and ANOVA_MODS_LINK. 54 Date: 01 December 2005 [rickr] 55 56 Mod: Added debug field to anova_options struct. 57 Date: 08 December 2005 [rickr] 58 59 Mod: Added aBdiff, Abdiff and abmean fields to anova_options struct. 60 Added ANOVA_BOUND() macro. 61 Date: 16 December 2005 [rickr] 62 63 Mod: Added mask and nmask fields to anova_options struct. 64 Date: 11 Mar 2009 [RWCox] 65 */ 66 67 /*---------------------------------------------------------------------------*/ 68 69 #include "mrilib.h" 70 #include "3ddata.h" 71 72 static char * commandline = NULL ; /* command line for history notes */ 73 74 /* documentation about modifications to ANOVA computations */ 75 #define ANOVA_MODS_LINK "https://afni.nimh.nih.gov/sscc/gangc/ANOVA_Mod.html" 76 77 /*** HP-UX ***/ 78 #ifdef HP 79 # define DF "bdf ." 80 #endif 81 82 /*** SGI IRIX ***/ 83 #ifdef SGI 84 # define DF "df -k ." 85 #endif 86 87 /*** DEC OSF1 ***/ 88 #ifdef OSF1 89 # define DF "df -k ." 90 #endif 91 92 /*** SunOS or Solaris ***/ 93 #if defined(SOLARIS) || defined(SUN) 94 # define DF "df -k" 95 #endif 96 97 /*** IBM RS6000 ***/ 98 #ifdef RS6000 99 #endif 100 101 /*** Linux 1.2.x ***/ 102 #ifdef LINUX 103 # define DF "df -k ." 104 #endif 105 106 /*** other? ***/ 107 #ifndef DF 108 # define DF "df" 109 #endif 110 111 112 #define MAX_LEVELS 300 /* max. number of factor levels */ 113 #define MAX_OBSERVATIONS 666 /* max. number of observations per cell */ 114 #define MAX_MEANS 50 /* max. number of user requested means */ 115 #define MAX_DIFFS 50 /* max. number of user requested diffs. */ 116 #define MAX_CONTR 75 /* max. number of user requested contrasts */ 117 #define MAX_NAME_LENGTH THD_MAX_NAME /* max. string length for file names */ 118 119 120 typedef struct anova_options 121 { 122 int datum; /* data type for "intensity" data subbrick */ 123 char session[MAX_NAME_LENGTH]; /* name of output directory */ 124 125 int diskspace; /* if positive, print out how much disk space 126 is required for program execution */ 127 128 int nvoxel; /* number of voxel for special output */ 129 130 int model; /* indicates type of ANOVA model to be used: 131 model=1 A,B,C fixed; AxBxC 132 model=2 A,B,C random; AxBxC 133 model=3 A fixed; B,C random; AxBxC 134 model=4 A,B fixed; C random; AxBxC 135 model=5 A,B fixed; C random; AxB,BxC,C(A) 136 */ 137 138 int old_method; /* flag indicating use of old functionality 139 bits: 001 = old_method 140 010 = OK 141 100 = assume sphericity 142 */ 143 144 int a; /* number of levels for factor A */ 145 int b; /* number of levels for factor B */ 146 int c; /* number of levels for factor C */ 147 int n; /* number of observations in each cell */ 148 int nt; /* total number of observations */ 149 150 int na[MAX_LEVELS]; /* number of observations at each level 151 of factor A (for 3dANOVA only) */ 152 153 char ***** xname; /* names of the input data files */ 154 char * first_dataset; /* name of the first data set */ 155 156 int nx, ny, nz; /* data set dimensions */ 157 int nxyz; /* number of voxels per image */ 158 159 int nftr; /* if positive, calculate F for treatment */ 160 char * ftrname; /* name of output file of F for treatment */ 161 162 int nfa; /* if positive, calculate F for A effect */ 163 char * faname; /* name of output file of F for A effect */ 164 int nfb; /* if positive, calculate F for B effect */ 165 char * fbname; /* name of output file of F for B effect */ 166 int nfc; /* if positive, calculate F for C effect */ 167 char * fcname; /* name of output file of F for C effect */ 168 int nfab; /* if pos., calculate F for A*B interaction */ 169 char * fabname; /* name of output file for A*B interaction */ 170 int nfac; /* if pos., calculate F for A*C interaction */ 171 char * facname; /* name of output file for A*C interaction */ 172 int nfbc; /* if pos., calculate F for B*C interaction */ 173 char * fbcname; /* name of output file for B*C interaction */ 174 int nfabc; /* if pos, calculate F for A*B*C interaction */ 175 char * fabcname; /* name of output file for A*B*C interaction */ 176 177 178 int num_ameans; /* number of factor A level means */ 179 int ameans[MAX_MEANS]; /* calc means for these factor A levels */ 180 char * amname[MAX_MEANS]; /* names of output files for factor A means */ 181 182 int num_bmeans; /* number of factor B level means */ 183 int bmeans[MAX_MEANS]; /* calc means for these factor B levels */ 184 char * bmname[MAX_MEANS]; /* names of output files for factor B means */ 185 186 int num_cmeans; /* number of factor C level means */ 187 int cmeans[MAX_MEANS]; /* calc means for these factor C levels */ 188 char * cmname[MAX_MEANS]; /* names of output files for factor C means */ 189 190 int num_xmeans; /* number of cell means */ 191 int xmeans[MAX_MEANS][3]; /* calc means for these cells */ 192 char * xmname[MAX_MEANS]; /* name of output files for cell means */ 193 194 int num_adiffs; /* num of diffs in factor A level means */ 195 int adiffs[MAX_DIFFS][2]; /* calc diffs in these factor A level means */ 196 char * adname[MAX_DIFFS]; /* names of output files for A differences */ 197 198 int num_bdiffs; /* num of diffs in factor B level means */ 199 int bdiffs[MAX_DIFFS][2]; /* calc diffs in these factor B level means */ 200 char * bdname[MAX_DIFFS]; /* names of output files for B differences */ 201 202 int num_cdiffs; /* num of diffs in factor C level means */ 203 int cdiffs[MAX_DIFFS][2]; /* calc diffs in these factor C level means */ 204 char * cdname[MAX_DIFFS]; /* names of output files for C differences */ 205 206 int num_xdiffs; /* num of diffs in cell means */ 207 int xdiffs[MAX_DIFFS][2][3]; /* calc diffs in these cell means */ 208 char * xdname[MAX_DIFFS]; /* names of output files for cell diffs */ 209 210 int num_acontr; /* number of factor A contrasts */ 211 float acontr[MAX_CONTR][MAX_LEVELS]; /* factor A contrast vectors */ 212 char * acname[MAX_CONTR]; /* names of output files for A contrasts */ 213 214 int num_bcontr; /* number of factor B contrasts */ 215 float bcontr[MAX_CONTR][MAX_LEVELS]; /* factor B contrast vectors */ 216 char * bcname[MAX_CONTR]; /* names of output files for B contrasts */ 217 218 int num_ccontr; /* number of factor C contrasts */ 219 float ccontr[MAX_CONTR][MAX_LEVELS]; /* factor C contrast vectors */ 220 char * ccname[MAX_CONTR]; /* names of output files for C contrasts */ 221 222 int num_xcontr; /* number of contrasts of cell means */ 223 float xcontr[MAX_CONTR][MAX_LEVELS][MAX_LEVELS]; 224 /* cell means contrast vectors */ 225 char * xcname[MAX_CONTR]; /* names of output files for cell contrasts */ 226 227 /* second order contrasts for 3dANOVA3 (e.g. -aBcontr) 14 Oct 2005 [rickr] */ 228 int num_aBcontr; /* number of A contrasts at fixed B level */ 229 float aBcontr[MAX_CONTR][MAX_LEVELS]; /* factor contrast vectors */ 230 int aBclevel[MAX_CONTR]; /* factor B level for aB contrasts */ 231 char * aBcname[MAX_CONTR]; /* names of output files for aB contrasts */ 232 233 int num_Abcontr; /* number of B contrasts at fixed A level */ 234 float Abcontr[MAX_CONTR][MAX_LEVELS]; /* factor contrast vectors */ 235 int Abclevel[MAX_CONTR]; /* factor A level for Ab contrasts */ 236 char * Abcname[MAX_CONTR]; /* names of output files for Ab contrasts */ 237 238 /* second order diffs for 3dANOVA3 (e.g. -aBdiff) 16 Dec 2005 [rickr] */ 239 int num_aBdiffs; /* number of A diffs at fixed B level */ 240 float aBdiffs[MAX_DIFFS][2]; /* calc diffs in these factor A level means */ 241 int aBdlevel[MAX_DIFFS]; /* fixed factor B level for each diff */ 242 char * aBdname[MAX_DIFFS]; /* names of output files for aB diffs */ 243 244 int num_Abdiffs; /* number of B diffs at fixed A level */ 245 float Abdiffs[MAX_DIFFS][2]; /* calc diffs in these factor B level means */ 246 int Abdlevel[MAX_DIFFS]; /* fixed factor A level for each diff */ 247 char * Abdname[MAX_DIFFS]; /* names of output files for Ab diffs */ 248 249 /* -abmeans: at fixed A level and fixed B level 16 Dec 2005 [rickr] */ 250 int num_abmeans; /* number of AB level means */ 251 int abmeans[MAX_MEANS][2]; /* calc means at each A-level B-level pair */ 252 char * abmname[MAX_MEANS]; /* names of output files for AB means */ 253 254 char * bucket_filename; /* file name for bucket dataset */ 255 256 int debug; /* for more verbose output */ 257 258 int nmask ; 259 byte *mask ; 260 } anova_options; 261 262 263 /*---------------------------------------------------------------------------*/ 264 /* 265 Routine to initialize input options. 266 */ 267 268 void initialize_options (anova_options * option_data); 269 270 271 /*---------------------------------------------------------------------------*/ 272 273 /** macro to open a dataset and make it ready for processing **/ 274 275 #define DOPEN(ds,name) \ 276 do{ int pv ; (ds) = THD_open_dataset((name)) ; \ 277 if( !ISVALID_3DIM_DATASET((ds)) ){ \ 278 fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \ 279 if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny || \ 280 (ds)->daxes->nzz!=nz ){ \ 281 fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; } \ 282 if( DSET_NVALS((ds)) != 1 ){ \ 283 fprintf(stderr,"*** Must specify 1 sub-brick: %s\n",(name));exit(1);}\ 284 THD_load_datablock( (ds)->dblk ) ; \ 285 pv = DSET_PRINCIPAL_VALUE((ds)) ; \ 286 if( DSET_ARRAY((ds),pv) == NULL ){ \ 287 fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); } \ 288 if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){ \ 289 fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1);\ 290 } \ 291 break ; } while (0) 292 293 294 /*---------------------------------------------------------------------------*/ 295 296 /** macro to return pointer to correct location in brick for current processing **/ 297 298 #define SUB_POINTER(ds,vv,ind,ptr) \ 299 do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){ \ 300 default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1); \ 301 case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ; \ 302 (ptr) = (void *)( fim + (ind) ) ; \ 303 } break ; \ 304 case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ; \ 305 (ptr) = (void *)( fim + (ind) ) ; \ 306 } break ; \ 307 case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ; \ 308 (ptr) = (void *)( fim + (ind) ) ; \ 309 } break ; } break ; } while(0) 310 311 312 /** macro to bound a number **/ 313 #define ANOVA_BOUND(var,bot,top) do{ \ 314 if((var)<(bot)) (var)=(bot); \ 315 else if((var)>(top)) (var)=(top); \ 316 } while (0) 317 318 /*---------------------------------------------------------------------------*/ 319 /* 320 Routine to print error message and stop. 321 */ 322 323 void ANOVA_error (char * message); 324 325 326 /*---------------------------------------------------------------------------*/ 327 /* 328 Routine to write a 3d volume of floating point data to a (temporary) 329 disk file. 330 */ 331 332 int volume_write_count (char * filename, float * fout, int size); 333 334 335 /*---------------------------------------------------------------------------*/ 336 /* 337 Routine to write a 3d volume of floating point data to a (temporary) 338 disk file. Error exit if cannot write entire file. 339 */ 340 341 void volume_write (char * filename, float * fout, int size); 342 343 344 /*---------------------------------------------------------------------------*/ 345 /* 346 Routine to read a 3d volume of floating point data. 347 */ 348 349 void volume_read (char * filename, float * fin, int size); 350 351 352 /*---------------------------------------------------------------------------*/ 353 /* 354 Routine to delete a file containing a 3d volume of floating point data. 355 */ 356 357 void volume_delete (char * filename); 358 359 360 /*---------------------------------------------------------------------------*/ 361 /* 362 Routine to set a 3d volume of floating point data to zero. 363 */ 364 365 void volume_zero (float * f, int size); 366 367 368 /*---------------------------------------------------------------------------*/ 369 /* 370 Routine to get the dimensions of the 3d AFNI data sets. 371 */ 372 373 void get_dimensions (anova_options * option_data); 374 375 376 /*---------------------------------------------------------------------------*/ 377 /* 378 Routine to check whether one temporary file already exists. 379 */ 380 381 void check_one_temporary_file (char * filename); 382 383 384 /*---------------------------------------------------------------------------*/ 385 /* 386 Routine to check whether one output file already exists. 387 */ 388 389 void check_one_output_file (anova_options * option_data, char * filename); 390 391 392 /*---------------------------------------------------------------------------*/ 393 /* 394 Return maximum of two integers. 395 */ 396 397 int max (int n1, int n2); 398 399 400 /*---------------------------------------------------------------------------*/ 401 /* 402 Routine to determine if sufficient disk space exists for storing 403 the temporary data files. 404 */ 405 406 void check_disk_space (anova_options * option_data); 407 408 409 /*---------------------------------------------------------------------------*/ 410 /* 411 Routine to read one AFNI data set from file 'filename'. 412 The data is converted to floating point (in ffim). 413 */ 414 415 void read_afni_data (anova_options * option_data, 416 char * filename, float * ffim); 417 418 419 /*---------------------------------------------------------------------------*/ 420 /* 421 Routine to write one AFNI data set. 422 423 This data set may be either 'fitt' type (intensity + t-statistic) 424 or 'fift' type (intensity + F-statistic). 425 426 The intensity data is in ffim, and the corresponding statistic is in ftr. 427 428 numdof = numerator degrees of freedom 429 dendof = denominator degrees of freedom 430 431 Note: if dendof = 0, then write 'fitt' type data set; 432 if dendof > 0, then write 'fift' type data set. 433 */ 434 435 void write_afni_data (anova_options * option_data, char * filename, 436 float * ffim, float * ftr, int numdof, int dendof); 437 438 439 /*---------------------------------------------------------------------------*/ 440 /* 441 Routine to add file name to command string. 442 */ 443 444 void add_file_name (THD_3dim_dataset * new_dset, char * filename, 445 char * command_str); 446 447 448 /*---------------------------------------------------------------------------*/ 449 /* 450 Routine to get a nice string to use as sub-brick label. 451 */ 452 char *label_from_filename(char *fname); 453 454 /*---------------------------------------------------------------------------*/ 455 /* 456 Routine to remove AFNI .HEAD and .BRIK dataset files. 457 */ 458 459 void remove_dataset (THD_3dim_dataset * new_dset, char * filename); 460 461 462 /*---------------------------------------------------------------------------*/ 463 /* 464 Routine to release memory allocated for anova_options. 465 */ 466 467 void destroy_anova_options (anova_options * option_data); 468 469 470 /*---------------------------------------------------------------------------*/ 471 /* 472 Routine to check for contrasts not summing to zero. 473 */ 474 475 int contrasts_are_valid (anova_options * option_data, int show_errs, int level); 476 477 478 /*---------------------------------------------------------------------------*/ 479 480 #undef ANOVA_FLOAT_HELP 481 #define ANOVA_FLOAT_HELP \ 482 "-------------------------------------------------------------------------\n" \ 483 "STORAGE FORMAT:\n" \ 484 "---------------\n" \ 485 "The default output format is to store the results as scaled short\n" \ 486 "(16-bit) integers. This truncantion might cause significant errors.\n" \ 487 "If you receive warnings that look like this:\n" \ 488 " *+ WARNING: TvsF[0] scale to shorts misfit = 8.09%% -- *** Beware\n" \ 489 "then you can force the results to be saved in float format by\n" \ 490 "defining the environment variable AFNI_FLOATIZE to be YES\n" \ 491 "before running the program. For convenience, you can do this\n" \ 492 "on the command line, as in\n" \ 493 " 3dANOVA -DAFNI_FLOATIZE=YES ... other options ... \n" \ 494 "Also see the following links:\n" \ 495 " https://afni.nimh.nih.gov/pub/dist/doc/program_help/common_options.html\n" \ 496 " https://afni.nimh.nih.gov/pub/dist/doc/program_help/README.environment.html\n" 497 498 /*---------------------------------------------------------------------------*/ 499