/***************************************************************************** Major portions of this software are copyrighted by the Medical College of Wisconsin, 1998-2001, and are released under the Gnu General Public License, Version 2. See the file README.Copyright for details. ******************************************************************************/ /*---------------------------------------------------------------------------*/ /* This file contains routines used by the three ANOVA programs: 3dANOVA, 3dANOVA2, and 3dANOVA3. File: 3dANOVA.lib Author: B. D. Ward Date: 5 January 1998 Mod: Added command -diskspace to print out how much disk space is required to execute the problem. Date: 23 January 1997 Mod: Change to routine write_afni_data to reduce truncation error. Date: 27 January 1997 Mod: Added machine specific code for checking disk space. Date: 29 January 1997 Mod: Extensive changes required to implement the 'bucket' dataset. Date: 30 December 1997 Mod: Divided file 3dANOVA.h into: 3dANOVA.h and 3dANOVA.lib. Date: 5 January 1998 Mod: Added software for statistical tests of individual cell means, cell differences, and cell contrasts. Date: 27 October 1998 Mod: Added changes for incorporating History notes. Date: 09 September 1999 Mod: Replaced dataset input code with calls to THD_open_dataset, to allow operator selection of individual sub-bricks for input. Date: 02 December 1999 Mod: Correction to routine destroy_anova_options. Removed memory deallocation for option_data->first_dataset. Date: 01 August 2001 Mod: Modified routine write_afni_data so that all output subbricks will now have the scaled short integer format. Date: 14 March 2002 Mod: Modify add_file_name() to use .1D filenames if so ordered. Date: 14 March 2003 - RWCox. */ /*---------------------------------------------------------------------------*/ /* Routine to initialize input options. */ void initialize_options (anova_options * option_data) { int i, j, k, m; option_data->datum = ILLEGAL_TYPE; strcpy (option_data->session, "./"); option_data->diskspace = 0; option_data->nvoxel = -1; option_data->model = 1; option_data->a = 0; option_data->b = 0; option_data->c = 0; option_data->n = 0; option_data->nt = 0; for (i = 0; i < MAX_LEVELS; i++) option_data->na[i] = 0; option_data->xname = NULL; option_data->first_dataset = NULL; option_data->nx = 0; option_data->ny = 0; option_data->nz = 0; option_data->nxyz = 0; option_data->nftr = 0; option_data->ftrname = NULL; option_data->nfa = 0; option_data->faname = NULL; option_data->nfb = 0; option_data->fbname = NULL; option_data->nfc = 0; option_data->fcname = NULL; option_data->nfab = 0; option_data->fabname = NULL; option_data->nfac = 0; option_data->facname = NULL; option_data->nfbc = 0; option_data->fbcname = NULL; option_data->nfabc = 0; option_data->fabcname = NULL; option_data->num_ameans = 0; option_data->num_bmeans = 0; option_data->num_cmeans = 0; option_data->num_xmeans = 0; for (i = 0; i < MAX_MEANS; i++) { option_data->amname[i] = NULL; option_data->bmname[i] = NULL; option_data->cmname[i] = NULL; option_data->xmname[i] = NULL; option_data->ameans[i] = 0; option_data->bmeans[i] = 0; option_data->cmeans[i] = 0; for (j = 0; j < 3; j++) option_data->xmeans[i][j] = 0; } option_data->num_adiffs = 0; option_data->num_bdiffs = 0; option_data->num_cdiffs = 0; option_data->num_xdiffs = 0; for (i = 0; i < MAX_DIFFS; i++) { option_data->adname[i] = NULL; option_data->bdname[i] = NULL; option_data->cdname[i] = NULL; option_data->xdname[i] = NULL; for (j = 0; j < 2; j++) { option_data->adiffs[i][j] = 0; option_data->bdiffs[i][j] = 0; option_data->cdiffs[i][j] = 0; for (k = 0; k < 3; k++) option_data->xdiffs[i][j][k] = 0; } } option_data->num_acontr = 0; option_data->num_bcontr = 0; option_data->num_ccontr = 0; option_data->num_xcontr = 0; for (i = 0; i < MAX_CONTR; i++) { option_data->acname[i] = NULL; option_data->bcname[i] = NULL; option_data->ccname[i] = NULL; option_data->xcname[i] = NULL; for (j = 0; j < MAX_LEVELS; j++) { option_data->acontr[i][j] = 0.0; option_data->bcontr[i][j] = 0.0; option_data->ccontr[i][j] = 0.0; for (k = 0; k < MAX_LEVELS; k++) option_data->xcontr[i][j][k] = 0.0; } } option_data->bucket_filename = NULL; } /*---------------------------------------------------------------------------*/ /* Routine to print error message and stop. */ void ANOVA_error (char * message) { fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message); exit(1); } /*---------------------------------------------------------------------------*/ /* Routine to write a 3d volume of floating point data to a (temporary) disk file. */ int volume_write_count (char * filename, float * fout, int size) { char sfilename[MAX_NAME_LENGTH]; /* output file name */ char message[MAX_NAME_LENGTH]; /* error message */ FILE * far; /* floating point output file */ int count; /* number of data items written to disk */ /*----- output file name -----*/ strcpy (sfilename, filename); strcat (sfilename, SUFFIX); /*----- first, see if file already exists -----*/ far = fopen (sfilename, "r"); if (far != NULL) { sprintf (message, "file %s already exists. ", sfilename); ANOVA_error (message); } /*----- open temporary data file for output -----*/ far = fopen (sfilename, "w"); if (far == NULL) ANOVA_error ("unable to write file "); /*----- write 3d data set -----*/ count = fwrite (fout, sizeof(float), size, far); fclose (far); return (count); } /*---------------------------------------------------------------------------*/ /* Routine to write a 3d volume of floating point data to a (temporary) disk file. Error exit if cannot write entire file. */ void volume_write (char * filename, float * fout, int size) { int count; /* number of data items written to disk */ /*----- attempt to write 3d volume of data -----*/ count = volume_write_count (filename, fout, size); /*----- error in writing file? -----*/ if (count != size) ANOVA_error ("error in writing data file "); } /*---------------------------------------------------------------------------*/ /* Routine to read a 3d volume of floating point data. */ void volume_read (char * filename, float * fin, int size) { char sfilename[MAX_NAME_LENGTH]; /* input file name */ char message[MAX_NAME_LENGTH]; /* error message */ FILE * far; /* floating point input file */ int count; /* number of data items read from disk */ /*----- input file name -----*/ strcpy (sfilename, filename); strcat (sfilename, SUFFIX); /*----- open temporary data file for input -----*/ far = fopen (sfilename, "r"); if (far == NULL) { sprintf (message, "Unable to open temporary file %s", sfilename); ANOVA_error (message); } /*----- read 3d data file -----*/ count = fread (fin, sizeof(float), size, far); fclose (far); /*----- error in reading file? -----*/ if (count != size) ANOVA_error ("error in reading data file "); } /*---------------------------------------------------------------------------*/ /* Routine to delete a file containing a 3d volume of floating point data. */ void volume_delete (char * filename) { char sfilename[MAX_NAME_LENGTH]; /* file name */ /*----- build file name -----*/ strcpy (sfilename, filename); strcat (sfilename, SUFFIX); /*----- delete 3d data file -----*/ remove (sfilename); } /*---------------------------------------------------------------------------*/ /* Routine to set a 3d volume of floating point data to zero. */ void volume_zero (float * f, int size) { int i; for (i = 0; i < size; i++) f[i] = 0.0; } /*---------------------------------------------------------------------------*/ /* Routine to get the dimensions of the 3d AFNI data sets. */ void get_dimensions (anova_options * option_data) { THD_3dim_dataset * dset=NULL; /*----- read first dataset to get dimensions, etc. -----*/ dset = THD_open_dataset( option_data->first_dataset ) ; if( ! ISVALID_3DIM_DATASET(dset) ){ fprintf(stderr,"*** Unable to open dataset file %s\n", option_data->first_dataset); exit(1) ; } /*----- data set dimensions in voxels -----*/ option_data->nx = dset->daxes->nxx ; option_data->ny = dset->daxes->nyy ; option_data->nz = dset->daxes->nzz ; option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ; THD_delete_3dim_dataset( dset , False ) ; dset = NULL ; } /*---------------------------------------------------------------------------*/ /* Routine to check whether one temporary file already exists. */ void check_one_temporary_file (char * filename) { FILE * far; /* temporary file pointer */ char sfilename[MAX_NAME_LENGTH]; /* temporary file name */ char message[MAX_NAME_LENGTH]; /* error message */ /*----- see if file already exists -----*/ strcpy (sfilename, filename); strcat (sfilename, SUFFIX); far = fopen (sfilename, "r"); if (far != NULL) { sprintf (message, "temporary file %s already exists. ", sfilename); ANOVA_error (message); } } /*---------------------------------------------------------------------------*/ /* Routine to check whether one output file already exists. */ void check_one_output_file (anova_options * option_data, char * filename) { THD_3dim_dataset * dset=NULL; /* input afni data set pointer */ THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ int ierror; /* number of errors in editing data */ /*----- read first dataset -----*/ dset = THD_open_dataset (option_data->first_dataset ) ; if( ! ISVALID_3DIM_DATASET(dset) ){ fprintf(stderr,"*** Unable to open dataset file %s\n", option_data->first_dataset); exit(1) ; } /*-- make an empty copy of this dataset --*/ new_dset = EDIT_empty_copy( dset ) ; THD_delete_3dim_dataset( dset , False ) ; dset = NULL ; ierror = EDIT_dset_items( new_dset , ADN_prefix , filename , ADN_label1 , filename , ADN_directory_name , option_data->session , ADN_self_name , filename , ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE , ADN_none ) ; if( ierror > 0 ){ fprintf(stderr, "*** %d errors in attempting to create output dataset!\n", ierror ) ; exit(1) ; } if( THD_is_file(new_dset->dblk->diskptr->header_name) ){ fprintf(stderr, "*** Output dataset file %s already exists--cannot continue!\a\n", new_dset->dblk->diskptr->header_name ) ; exit(1) ; } /*----- deallocate memory -----*/ THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; } /*---------------------------------------------------------------------------*/ /* Return maximum of two integers. */ int max (int n1, int n2) { if (n1 > n2) return n1; else return n2; } /*---------------------------------------------------------------------------*/ /* Routine to determine if sufficient disk space exists for storing the temporary data files. */ void check_disk_space (anova_options * option_data) { char ch; /* user response */ int nxyz; /* number of voxels per image */ int nmax; /* maximum number of disk files */ char filename[MAX_NAME_LENGTH]; /* output file name */ /*----- initialize local variables -----*/ nxyz = option_data->nxyz; /*----- first, determine the maximum number of files required -----*/ nmax = required_data_files (option_data); printf("\nThis problem requires approximately %d MB of free disk space.\n", (4*nmax*nxyz/1000000) + 1); /*----- Determine how much disk space is available. -----*/ printf ("Summary of available disk space: \n\n"); system (DF); printf ("\nDo you wish to proceed? "); ch = getchar(); if ((ch != 'Y') && (ch != 'y')) exit(0); printf ("\n"); } /*---------------------------------------------------------------------------*/ /* Routine to read one AFNI data set from file 'filename'. The data is converted to floating point (in ffim). */ void read_afni_data (anova_options * option_data, char * filename, float * ffim) { int iv; /* index number of intensity sub-brick */ THD_3dim_dataset * dset=NULL; /* data set pointer */ void * vfim = NULL; /* image data pointer */ int nx, ny, nz, nxyz; /* data set dimensions in voxels */ nx = option_data->nx; ny = option_data->ny; nz = option_data->nz; nxyz = option_data->nxyz; /*----- read in the data -----*/ DOPEN(dset,filename) ; iv = DSET_PRINCIPAL_VALUE(dset) ; /*----- convert it to floats (in ffim) -----*/ SUB_POINTER(dset,iv,0,vfim) ; EDIT_coerce_scale_type( nxyz , DSET_BRICK_FACTOR(dset,iv) , DSET_BRICK_TYPE(dset,iv),vfim , /* input */ MRI_float ,ffim ) ; /* output */ THD_delete_3dim_dataset( dset , False ) ; dset = NULL ; } /*---------------------------------------------------------------------------*/ /* Convert one volume to another type, autoscaling: nxy = # voxels itype = input datum type ivol = pointer to input volume otype = output datum type ovol = pointer to output volume (again, must be pre-malloc-ed) Return value is the scaling factor used (0.0 --> no scaling). */ float EDIT_coerce_autoscale_new( int nxyz , int itype,void *ivol , int otype,void *ovol ) { float fac=0.0 , top ; if( MRI_IS_INT_TYPE(otype) ){ top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ; if (top == 0.0) fac = 0.0; else fac = MRI_TYPE_maxval[otype]/top; } EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ; return ( fac ); } /*---------------------------------------------------------------------------*/ /* Routine to write one AFNI data set. This data set may be either 'fitt' type (intensity + t-statistic) or 'fift' type (intensity + F-statistic). The intensity data is in ffim, and the corresponding statistic is in ftr. numdof = numerator degrees of freedom dendof = denominator degrees of freedom Note: if dendof = 0, then write 'fitt' type data set; if dendof > 0, then write 'fift' type data set. */ void write_afni_data (anova_options * option_data, char * filename, float * ffim, float * ftr, int numdof, int dendof) { int nxyz; /* number of voxels */ int ii; /* voxel index */ THD_3dim_dataset * dset=NULL; /* input afni data set pointer */ THD_3dim_dataset * new_dset=NULL; /* output afni data set pointer */ int ierror; /* number of errors in editing data */ int ibuf[32]; /* integer buffer */ float fbuf[MAX_STAT_AUX]; /* float buffer */ float fimfac; /* scale factor for short data */ int output_datum; /* data type for output data */ short * tsp; /* 2nd sub-brick data pointer */ void * vdif; /* 1st sub-brick data pointer */ int func_type; /* afni data set type */ float top, func_scale_short; /* parameters for scaling data */ char label[80]; /* label for output file history */ /*----- initialize local variables -----*/ nxyz = option_data->nxyz; /*----- read first dataset -----*/ dset = THD_open_dataset (option_data->first_dataset) ; if( ! ISVALID_3DIM_DATASET(dset) ){ fprintf(stderr,"*** Unable to open dataset file %s\n", option_data->first_dataset); exit(1) ; } /*-- make an empty copy of this dataset, for eventual output --*/ new_dset = EDIT_empty_copy( dset ) ; /*----- Record history of dataset -----*/ sprintf (label, "Output prefix: %s", filename); if( commandline != NULL ) tross_multi_Append_History( new_dset , commandline,label,NULL ) ; else tross_Append_History ( new_dset, label); output_datum = MRI_short ; ibuf[0] = output_datum ; ibuf[1] = MRI_short ; if (dendof == 0) func_type = FUNC_TT_TYPE; else func_type = FUNC_FT_TYPE; ierror = EDIT_dset_items( new_dset , ADN_prefix , filename , ADN_label1 , filename , ADN_directory_name , option_data->session , ADN_self_name , filename , ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : GEN_FUNC_TYPE , ADN_func_type , func_type , ADN_nvals , FUNC_nvals[func_type] , ADN_datum_array , ibuf , ADN_malloc_type, DATABLOCK_MEM_MALLOC , ADN_none ) ; if( ierror > 0 ){ fprintf(stderr, "*** %d errors in attempting to create output dataset!\n", ierror ) ; exit(1) ; } if( THD_is_file(new_dset->dblk->diskptr->header_name) ){ fprintf(stderr, "*** Output dataset file %s already exists--cannot continue!\a\n", new_dset->dblk->diskptr->header_name ) ; exit(1) ; } /*----- deleting exemplar dataset -----*/ THD_delete_3dim_dataset( dset , False ) ; dset = NULL ; /*----- allocate memory for output data -----*/ vdif = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ; tsp = (short *) malloc( sizeof(short) * nxyz ) ; /*----- attach bricks to new data set -----*/ mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0)); mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1)); /*----- convert data type to output specification -----*/ fimfac = EDIT_coerce_autoscale_new (nxyz, MRI_float, ffim, output_datum, vdif); if (fimfac != 0.0) fimfac = 1.0 / fimfac; #define TOP_SS 32700 if (dendof == 0) /* t-statistic */ { top = TOP_SS/FUNC_TT_SCALE_SHORT; func_scale_short = FUNC_TT_SCALE_SHORT; } else /* F-statistic */ { top = TOP_SS/FUNC_FT_SCALE_SHORT; func_scale_short = FUNC_FT_SCALE_SHORT; } for (ii = 0; ii < nxyz; ii++) { if (ftr[ii] > top) tsp[ii] = TOP_SS; else if (ftr[ii] < -top) tsp[ii] = -TOP_SS; else if (ftr[ii] >= 0.0) tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5); else tsp[ii] = (short) (func_scale_short * ftr[ii] - 0.5); } /*----- write afni data set -----*/ printf("--- Writing combined dataset into %s\n", new_dset->dblk->diskptr->header_name) ; fbuf[0] = numdof; fbuf[1] = dendof; for( ii=2 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ; (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ; fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ; fbuf[1] = 1.0 / func_scale_short ; (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ; THD_load_statistics( new_dset ) ; THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ; /*----- deallocate memory -----*/ THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ; } /*---------------------------------------------------------------------------*/ static int use_1D_filenames = 0 ; /* 14 Mar 2003 */ static void USE_1D_filenames( int uu ){ use_1D_filenames=uu; } /*---------------------------------------------------------------------------*/ /* Routine to add file name to command string. */ void add_file_name (THD_3dim_dataset * new_dset, char * filename, char * command_str) { EDIT_dset_items (new_dset, ADN_prefix, filename, ADN_none); if (! ISVALID_3DIM_DATASET(new_dset)) { fprintf (stderr,"*** Unable to open dataset file %s\n", DSET_HEADNAME(new_dset)); exit(1); } strcat (command_str, " "); strcat (command_str, DSET_HEADNAME(new_dset)); switch( use_1D_filenames ){ /* 14 Mar 2003 - RWCox */ default: break ; case 1: case 3:{ int ll = strlen(command_str) ; command_str[ll-10] = '\0' ; strcat(command_str,".") ; if( use_1D_filenames == 1 ) strcat(command_str,"1") ; else strcat(command_str,"3") ; strcat(command_str,"D") ; } break ; } } /*---------------------------------------------------------------------------*/ /* Routine to remove AFNI .HEAD and .BRIK dataset files. */ void remove_dataset (THD_3dim_dataset * new_dset, char * filename) { EDIT_dset_items (new_dset, ADN_prefix, filename, ADN_none); if( !ISVALID_DSET(new_dset) ) return ; #if 0 fprintf(stderr,"remove_dataset(%s): ",filename) ; #endif switch( use_1D_filenames ){ /* 14 Mar 2003 - RWCox */ default: remove (DSET_HEADNAME(new_dset)); remove (DSET_BRIKNAME(new_dset)); #if 0 fprintf(stderr,"%s ",DSET_HEADNAME(new_dset)) ; fprintf(stderr,"%s\n",DSET_BRIKNAME(new_dset)) ; #endif case 1: case 3:{ char *fname=strdup(DSET_BRIKNAME(new_dset)); int ll=strlen(fname); if( ll > 10 && strstr(fname,".BRIK") != NULL ){ fname[ll-10] = '\0' ; strcat(fname,".") ; if( use_1D_filenames == 1 ) strcat(fname,"1") ; else strcat(fname,"3") ; strcat(fname,"D") ; } #if 0 fprintf(stderr,"%s\n",fname); #endif remove(fname) ; free(fname) ; } break ; } } /*---------------------------------------------------------------------------*/ /* Routine to release memory allocated for anova_options. */ void destroy_anova_options (anova_options * option_data) { int i, j, k, m; /*----- deallocate memory -----*/ for (i = 0; i < option_data->a; i++) for (j = 0; j < option_data->b; j++) for (k = 0; k < option_data->c; k++) for (m = 0; m < option_data->n; m++) { free (option_data->xname[i][j][k][m]); option_data->xname[i][j][k][m] = NULL; } if (option_data->nftr) { free (option_data->ftrname); option_data->ftrname = NULL; } if (option_data->nfa) { free (option_data->faname); option_data->faname = NULL; } if (option_data->nfb) { free (option_data->fbname); option_data->fbname = NULL; } if (option_data->nfc) { free (option_data->fcname); option_data->fcname = NULL; } if (option_data->nfab) { free (option_data->fabname); option_data->fabname = NULL; } if (option_data->nfac) { free (option_data->facname); option_data->facname = NULL; } if (option_data->nfbc) { free (option_data->fbcname); option_data->fbcname = NULL; } if (option_data->nfabc) { free (option_data->fabcname); option_data->fabcname = NULL; } if (option_data->num_ameans) for (i=0; i < option_data->num_ameans; i++) { free (option_data->amname[i]); option_data->amname[i] = NULL; } if (option_data->num_bmeans) for (i=0; i < option_data->num_bmeans; i++) { free (option_data->bmname[i]); option_data->bmname[i] = NULL; } if (option_data->num_cmeans) for (i=0; i < option_data->num_cmeans; i++) { free (option_data->cmname[i]); option_data->cmname[i] = NULL; } if (option_data->num_xmeans) for (i=0; i < option_data->num_xmeans; i++) { free (option_data->xmname[i]); option_data->xmname[i] = NULL; } if (option_data->num_adiffs) for (i=0; i < option_data->num_adiffs; i++) { free (option_data->adname[i]); option_data->adname[i] = NULL; } if (option_data->num_bdiffs) for (i=0; i < option_data->num_bdiffs; i++) { free (option_data->bdname[i]); option_data->bdname[i] = NULL; } if (option_data->num_cdiffs) for (i=0; i < option_data->num_cdiffs; i++) { free (option_data->cdname[i]); option_data->cdname[i] = NULL; } if (option_data->num_xdiffs) for (i=0; i < option_data->num_xdiffs; i++) { free (option_data->xdname[i]); option_data->xdname[i] = NULL; } if (option_data->num_acontr) for (i=0; i < option_data->num_acontr; i++) { free (option_data->acname[i]); option_data->acname[i] = NULL; } if (option_data->num_bcontr) for (i=0; i < option_data->num_bcontr; i++) { free (option_data->bcname[i]); option_data->bcname[i] = NULL; } if (option_data->num_ccontr) for (i=0; i < option_data->num_ccontr; i++) { free (option_data->ccname[i]); option_data->ccname[i] = NULL; } if (option_data->num_xcontr) for (i=0; i < option_data->num_xcontr; i++) { free (option_data->xcname[i]); option_data->xcname[i] = NULL; } }