00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 2002, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*---------------------------------------------------------------------------*/
00008 /*
00009   This program implements the False Discovery Rate (FDR) algorithm for
00010   thresholding of voxelwise statistics.  
00011 
00012   File:    3dFDR.c
00013   Author:  B. Douglas Ward
00014   Date:    31 January 2002
00015 
00016 */
00017 /*---------------------------------------------------------------------------*/
00018 
00019 #define PROGRAM_NAME "3dFDR"                         /* name of this program */
00020 #define PROGRAM_AUTHOR "B. Douglas Ward"                   /* program author */
00021 #define PROGRAM_INITIAL "31 January 2002" /* date of initial program release */
00022 #define PROGRAM_LATEST "31 January 2002"    /* date of last program revision */
00023 
00024 
00025 /*---------------------------------------------------------------------------*/
00026 /*
00027   Include header files and source code files.
00028 */
00029 
00030 #include "mrilib.h"
00031 
00032 
00033 /*---------------------------------------------------------------------------*/
00034 /*
00035   Structure declarations 
00036 */
00037 
00038 struct voxel;
00039 
00040 typedef struct voxel
00041 {
00042   int ixyz;                            /* voxel index */
00043   float pvalue;                        /* input p-value or output q-value */
00044   struct voxel * next_voxel;           /* pointer to next voxel in list */
00045 } voxel;
00046 
00047 
00048 
00049 /*-------------------------- global data ------------------------------------*/
00050 
00051 #define  FDR_MAX_LL 10000        /* maximum number of linked lists of voxels */
00052 
00053 
00054 static int    FDR_quiet      = 0;      /* flag for suppress screen output */
00055 static int    FDR_list       = 0;      /* flag for list voxel q-values */
00056 static float  FDR_mask_thr   = 1.0;    /* mask threshold */
00057 static float  FDR_cn         = 1.0;    /* statistical constant c(N) */
00058 static int    FDR_nxyz       = 0;      /* dataset dimensions in voxels */
00059 static int    FDR_nthr       = 0;      /* number of voxels in mask */
00060 
00061 static char * FDR_input_filename   = NULL;   /* input 3d functional dataset */
00062 static char * FDR_input1D_filename = NULL;   /* input list of p-values */
00063 static char * FDR_mask_filename    = NULL;   /* input 3d mask dataset */
00064 static char * FDR_output_prefix    = NULL;   /* name for output 3d dataset */
00065 
00066 static byte  * FDR_mask = NULL;            /* mask for voxels above thr. */
00067 static float * FDR_input1D_data = NULL;    /* input array of p-values */
00068 static voxel * FDR_head_voxel[FDR_MAX_LL]; /* linked lists of voxels */
00069 
00070 static char * commandline = NULL;          /* command line for history notes */
00071 static THD_3dim_dataset * FDR_dset = NULL; /* input dataset */  
00072 
00073 
00074 /*---------------------------------------------------------------------------*/
00075 
00076 /** macro to open a dataset and make it ready for processing **/
00077 
00078 #define DOPEN(ds,name)                                                        \
00079 do{ int pv ; (ds) = THD_open_dataset((name)) ;                                \
00080        if( !ISVALID_3DIM_DATASET((ds)) ){                                     \
00081           fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
00082        THD_load_datablock( (ds)->dblk ) ;                                     \
00083        pv = DSET_PRINCIPAL_VALUE((ds)) ;                                      \
00084        if( DSET_ARRAY((ds),pv) == NULL ){                                     \
00085           fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); }   \
00086        if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){                         \
00087           fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1);\
00088        }     \
00089        break ; } while (0)
00090 
00091 
00092 /*---------------------------------------------------------------------------*/
00093 
00094 /** macro to return pointer to correct location in brick for current processing **/
00095 
00096 #define SUB_POINTER(ds,vv,ind,ptr)                                            \
00097    do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){                                  \
00098          default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1);       \
00099             case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ;  \
00100                             (ptr) = (void *)( fim + (ind) ) ;                 \
00101             } break ;                                                         \
00102             case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ;     \
00103                             (ptr) = (void *)( fim + (ind) ) ;                 \
00104             } break ;                                                         \
00105             case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ;  \
00106                              (ptr) = (void *)( fim + (ind) ) ;                \
00107             } break ; } break ; } while(0)
00108 
00109 
00110 /*---------------------------------------------------------------------------*/
00111 /*
00112    Print error message and stop.
00113 */
00114 
00115 void FDR_error (char * message)
00116 {
00117   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00118   exit(1);
00119 }
00120 
00121 
00122 /*---------------------------------------------------------------------------*/
00123 
00124 /** macro to test a malloc-ed pointer for validity **/
00125 
00126 #define MTEST(ptr) \
00127 if((ptr)==NULL) \
00128 ( FDR_error ("Cannot allocate memory") )
00129      
00130 
00131 /*---------------------------------------------------------------------------*/
00132 /*
00133   Display help file.
00134 */
00135 
00136 void FDR_Syntax(void)
00137 {
00138 printf(
00139 "This program implements the False Discovery Rate (FDR) algorithm for       \n"
00140 "thresholding of voxelwise statistics.                                      \n"
00141 "                                                                           \n"
00142 "Program input consists of a functional dataset containing one (or more)    \n"
00143 "statistical sub-bricks.  Output consists of a bucket dataset with one      \n"
00144 "sub-brick for each input sub-brick.  For non-statistical input sub-bricks, \n"
00145 "the output is a copy of the input.  However, statistical input sub-bricks  \n"
00146 "are replaced by their corresponding FDR values, as follows:                \n"
00147 "                                                                           \n"
00148 "For each voxel, the minimum value of q is determined such that             \n"
00149 "                               E(FDR) <= q                                 \n"
00150 "leads to rejection of the null hypothesis in that voxel. Only voxels inside\n"
00151 "the user specified mask will be considered.  These q-values are then mapped\n"
00152 "to z-scores for compatibility with the AFNI statistical threshold display: \n"
00153 "                                                                           \n"
00154 "               stat ==> p-value ==> FDR q-value ==> FDR z-score            \n"
00155 "                                                                           \n"
00156 );
00157 
00158 printf(
00159 "Usage:                                                                     \n"
00160 "  3dFDR                                                                    \n"
00161 "    -input fname       fname = filename of input 3d functional dataset     \n"
00162 "      OR                                                                   \n"
00163 "    -input1D dname     dname = .1D file containing column of p-values      \n"
00164 "                                                                           \n"
00165 "    -mask_file mname   Use mask values from file mname.                    \n"
00166 "                       Note: If file mname contains more than 1 sub-brick, \n"
00167 "                       the mask sub-brick must be specified!               \n"
00168 "                       Default: No mask                                    \n"
00169 "                                                                           \n"
00170 "    -mask_thr m        Only voxels whose corresponding mask value is       \n"
00171 "                       greater than or equal to m in absolute value will   \n"
00172 "                       be considered.  Default: m=1                        \n"
00173 "                                                                           \n"
00174 "                       Constant c(N) depends on assumption about p-values: \n"
00175 "    -cind              c(N) = 1   p-values are independent across N voxels \n"
00176 "    -cdep              c(N) = sum(1/i), i=1,...,N   any joint distribution \n"
00177 "                       Default:  c(N) = 1                                  \n"
00178 "                                                                           \n"
00179 "    -quiet             Flag to suppress screen output                      \n"
00180 "                                                                           \n"
00181 "    -list              Write sorted list of voxel q-values to screen       \n"
00182 "                                                                           \n"
00183 "    -prefix pname      Use 'pname' for the output dataset prefix name.     \n"
00184 "      OR                                                                   \n"
00185 "    -output pname                                                          \n"
00186 "                                                                           \n"
00187 "\n") ;
00188 
00189    printf("\n" MASTER_SHORTHELP_STRING ) ;
00190 
00191    exit(0) ;
00192 
00193 }
00194 
00195 
00196 /*---------------------------------------------------------------------------*/
00197 /*
00198    Read the arguments, load the global variables
00199 
00200 */
00201 
00202 void read_options ( int argc , char * argv[] )
00203 {
00204   int nopt = 1 ;           /* count of input arguments */
00205   char message[80];        /* error message */
00206 
00207 
00208   
00209   /*----- main loop over input options -----*/
00210   while( nopt < argc )
00211     {
00212 
00213       /*-----   -input fname   -----*/
00214       if (strcmp(argv[nopt], "-input") == 0)
00215         {
00216           nopt++;
00217           if (nopt >= argc)  FDR_error ("need argument after -input ");
00218           FDR_input_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME);
00219           MTEST (FDR_input_filename);
00220           strcpy (FDR_input_filename, argv[nopt]);
00221           nopt++;
00222           continue;
00223         }
00224       
00225       
00226       /*-----   -input1D dname   -----*/
00227       if (strcmp(argv[nopt], "-input1D") == 0)
00228         {
00229           nopt++;
00230           if (nopt >= argc)  FDR_error ("need argument after -input1D ");
00231           FDR_input1D_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME);
00232           MTEST (FDR_input1D_filename);
00233           strcpy (FDR_input1D_filename, argv[nopt]);
00234           nopt++;
00235           continue;
00236         }
00237       
00238       
00239       /*-----   -mask_file mname   -----*/
00240       if (strcmp(argv[nopt], "-mask_file") == 0)
00241         {
00242           nopt++;
00243           if (nopt >= argc)  FDR_error ("need argument after -mask_file ");
00244           FDR_mask_filename = (char *) malloc (sizeof(char)*THD_MAX_NAME);
00245           MTEST (FDR_mask_filename);
00246           strcpy (FDR_mask_filename, argv[nopt]);
00247           nopt++;
00248           continue;
00249         }
00250       
00251 
00252       /*----- -mask_thr m -----*/
00253       if( strcmp(argv[nopt],"-mask_thr") == 0 ){
00254          float fval;
00255          nopt++ ;
00256          if( nopt >= argc ){
00257             FDR_error (" need 1 argument after -mask_thr"); 
00258          }
00259          sscanf (argv[nopt], "%f", &fval); 
00260          if (fval < 0.0){
00261             FDR_error (" Require mask_thr >= 0.0 ");
00262          }
00263          FDR_mask_thr = fval;
00264          nopt++;  continue;
00265       }
00266 
00267 
00268       /*----- -cind -----*/
00269       if( strcmp(argv[nopt],"-cind") == 0 ){
00270          FDR_cn = 1.0;
00271          nopt++ ; continue ;
00272       }
00273 
00274       
00275       /*----- -cdep -----*/
00276       if( strcmp(argv[nopt],"-cdep") == 0 ){
00277          FDR_cn = -1.0;
00278          nopt++ ; continue ;
00279       }
00280 
00281       
00282       /*----- -quiet -----*/
00283       if( strcmp(argv[nopt],"-quiet") == 0 ){
00284          FDR_quiet = 1;
00285          nopt++ ; continue ;
00286       }
00287 
00288       
00289       /*----- -list -----*/
00290       if( strcmp(argv[nopt],"-list") == 0 ){
00291          FDR_list = 1;
00292          nopt++ ; continue ;
00293       }
00294 
00295       
00296       /*----- -prefix prefix -----*/
00297       if( strcmp(argv[nopt],"-prefix") == 0 ||
00298           strcmp(argv[nopt],"-output") == 0   ){
00299          nopt++ ;
00300          if( nopt >= argc ){
00301             FDR_error (" need argument after -prefix!");
00302          }
00303          FDR_output_prefix = (char *) malloc (sizeof(char) * THD_MAX_PREFIX); 
00304          MCW_strncpy( FDR_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
00305          continue ;
00306       }
00307 
00308 
00309       /*----- unknown command -----*/
00310       sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
00311       FDR_error (message);
00312       
00313 
00314    }  /*----- end of loop over command line arguments -----*/
00315 
00316 
00317 }
00318 
00319 
00320 /*---------------------------------------------------------------------------*/
00321 /*
00322   Read time series from specified file name.  This file name may have
00323   a column selector attached.
00324 */
00325 
00326 float * read_time_series 
00327 (
00328   char * ts_filename,          /* time series file name (plus column index) */
00329   int * ts_length              /* output value for time series length */
00330 )
00331 
00332 {
00333   char message[THD_MAX_NAME];    /* error message */
00334   char * cpt;                    /* pointer to column suffix */
00335   char filename[THD_MAX_NAME];   /* time series file name w/o column index */
00336   char subv[THD_MAX_NAME];       /* string containing column index */
00337   MRI_IMAGE * im, * flim;  /* pointers to image structures 
00338                               -- used to read 1D ASCII */
00339   float * far;             /* pointer to MRI_IMAGE floating point data */
00340   int nx;                  /* number of time points in time series */
00341   int ny;                  /* number of columns in time series file */
00342   int iy;                  /* time series file column index */
00343   int ipt;                 /* time point index */
00344   float * ts_data = NULL;  /* input time series data */
00345 
00346 
00347   /*----- First, check for empty filename -----*/
00348   if (ts_filename == NULL)
00349     FDR_error ("Missing input time series file name");
00350 
00351 
00352   /*----- Read the time series file -----*/
00353   flim = mri_read_1D(ts_filename) ;
00354   if (flim == NULL)
00355     {
00356       sprintf (message,  "Unable to read time series file: %s",  ts_filename);
00357       FDR_error (message);
00358     }
00359 
00360   far = MRI_FLOAT_PTR(flim);
00361   nx = flim->nx;
00362   ny = flim->ny; iy = 0 ;
00363   if( ny > 1 ){
00364     fprintf(stderr,"WARNING: time series %s has more than 1 column\n",ts_filename);
00365   }
00366   
00367 
00368   /*----- Save the time series data -----*/
00369   *ts_length = nx;
00370   ts_data = (float *) malloc (sizeof(float) * nx);
00371   MTEST (ts_data);
00372   for (ipt = 0;  ipt < nx;  ipt++)
00373     ts_data[ipt] = far[ipt + iy*nx];   
00374   
00375   
00376   mri_free (flim);  flim = NULL;
00377 
00378   return (ts_data);
00379 }
00380 
00381 
00382 /*---------------------------------------------------------------------------*/
00383 /*
00384   Routine to check whether one output file already exists.
00385 */
00386 
00387 void check_one_output_file 
00388 (
00389   THD_3dim_dataset * dset_time,     /* input 3d+time data set */
00390   char * filename                   /* name of output file */
00391 )
00392 
00393 {
00394   char message[THD_MAX_NAME];      /* error message */
00395   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00396   int ierror;                         /* number of errors in editing data */
00397 
00398   
00399   /*----- make an empty copy of input dataset -----*/
00400   new_dset = EDIT_empty_copy( dset_time ) ;
00401   
00402   
00403   ierror = EDIT_dset_items( new_dset ,
00404                             ADN_prefix , filename ,
00405                             ADN_label1 , filename ,
00406                             ADN_self_name , filename ,
00407                             ADN_type , ISHEAD(dset_time) ? HEAD_FUNC_TYPE : 
00408                                                            GEN_FUNC_TYPE ,
00409                             ADN_none ) ;
00410   
00411   if( ierror > 0 )
00412     {
00413       sprintf (message,
00414                "*** %d errors in attempting to create output dataset!\n", 
00415                ierror);
00416       FDR_error (message);
00417     }
00418   
00419   if( THD_is_file(new_dset->dblk->diskptr->header_name) )
00420     {
00421       sprintf (message,
00422                "Output dataset file %s already exists "
00423                " -- cannot continue! ",
00424                new_dset->dblk->diskptr->header_name);
00425       FDR_error (message);
00426     }
00427   
00428   /*----- deallocate memory -----*/   
00429   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00430   
00431 }
00432 
00433 
00434 /*---------------------------------------------------------------------------*/
00435 /*
00436   Routine to initialize the program: get all operator inputs; create mask
00437   for voxels above mask threshold.
00438 */
00439 
00440 void * initialize_program (int argc, char * argv[])
00441 {
00442   int iv;                  /* index number of sub-brick */
00443   void * vfim = NULL;      /* sub-brick data pointer */
00444   float * ffim = NULL;     /* sub-brick data in floating point format */
00445   int ixyz;                /* voxel index */
00446   int nx, ny, nz, nxyz;    /* numbers of voxels in input dataset */
00447   int mx, my, mz, mxyz;    /* numbers of voxels in mask dataset */
00448   int nthr;                /* number of voxels above mask threshold */
00449   char message[80];        /* error message */
00450   int ibin;                /* p-value bin index */
00451 
00452 
00453   /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/
00454   machdep() ; 
00455   { int new_argc ; char ** new_argv ;
00456   addto_args( argc , argv , &new_argc , &new_argv ) ;
00457   if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
00458   }
00459   
00460 
00461   /*----- Save command line for history notes -----*/
00462   commandline = tross_commandline( PROGRAM_NAME , argc,argv ) ;
00463 
00464 
00465   /*----- Does user request help menu? -----*/
00466   if( argc < 2 || strcmp(argv[1],"-help") == 0 ) FDR_Syntax() ;
00467 
00468   
00469   /*----- Add to program log -----*/
00470   AFNI_logger (PROGRAM_NAME,argc,argv); 
00471 
00472 
00473   /*----- Read input options -----*/
00474   read_options( argc , argv ) ;
00475 
00476 
00477   /*----- Open the mask dataset -----*/
00478   if (FDR_mask_filename != NULL)
00479     {
00480       if (!FDR_quiet) 
00481         printf ("Reading mask dataset: %s \n", FDR_mask_filename);
00482       DOPEN (FDR_dset, FDR_mask_filename);
00483 
00484       if (FDR_dset == NULL)
00485         {
00486           sprintf (message, "Cannot open mask dataset %s", FDR_mask_filename); 
00487           FDR_error (message);
00488         }
00489 
00490       if (DSET_NVALS(FDR_dset) != 1)
00491         FDR_error ("Must specify single sub-brick for mask data");
00492 
00493 
00494       /*----- Get dimensions of mask dataset -----*/
00495       mx   = DSET_NX(FDR_dset);   
00496       my   = DSET_NY(FDR_dset);   
00497       mz   = DSET_NZ(FDR_dset);
00498       mxyz = mx*my*mz;
00499 
00500 
00501       /*----- Allocate memory for float data -----*/
00502       ffim = (float *) malloc (sizeof(float) * mxyz);   MTEST (ffim);
00503 
00504 
00505       /*----- Convert mask dataset sub-brick to floats (in ffim) -----*/
00506       iv = DSET_PRINCIPAL_VALUE (FDR_dset);
00507       SUB_POINTER (FDR_dset, iv, 0, vfim);
00508       EDIT_coerce_scale_type (mxyz, DSET_BRICK_FACTOR(FDR_dset,iv),
00509                               DSET_BRICK_TYPE(FDR_dset,iv), vfim,  /* input  */
00510                               MRI_float                   , ffim); /* output */
00511   
00512       
00513       /*----- Allocate memory for mask volume -----*/
00514       FDR_mask = (byte *) malloc (sizeof(byte) * mxyz);
00515       MTEST (FDR_mask);
00516       
00517       
00518       /*----- Create mask of voxels above mask threshold -----*/
00519       nthr = 0;
00520       for (ixyz = 0;  ixyz < mxyz;  ixyz++)
00521         if (fabs(ffim[ixyz]) >= FDR_mask_thr)  
00522           { 
00523             FDR_mask[ixyz] = 1;
00524             nthr++;
00525           }
00526         else
00527           FDR_mask[ixyz] = 0;
00528 
00529       if (!FDR_quiet)  
00530         printf ("Number of voxels above mask threshold = %d \n", nthr);
00531       if (nthr < 1)  
00532         FDR_error ("No voxels above mask threshold.  Cannot continue.");
00533 
00534 
00535       /*----- Delete floating point sub-brick -----*/
00536       if (ffim != NULL) { free (ffim);   ffim = NULL; }
00537 
00538       /*----- Delete mask dataset -----*/
00539       THD_delete_3dim_dataset (FDR_dset, False);  FDR_dset = NULL ;
00540 
00541     }
00542 
00543 
00544   /*----- Get the input data -----*/
00545 
00546   if (FDR_input1D_filename != NULL)
00547     {
00548       /*----- Read the input .1D file -----*/
00549       if (!FDR_quiet)  printf ("Reading input data: %s \n", 
00550                                FDR_input1D_filename);
00551       FDR_input1D_data = read_time_series (FDR_input1D_filename, &nxyz);
00552 
00553       if (FDR_input1D_data == NULL)  
00554         { 
00555           sprintf (message,  "Unable to read input .1D data file: %s", 
00556                    FDR_input1D_filename);
00557           FDR_error (message);
00558         }
00559       
00560       if (nxyz < 1)  
00561         { 
00562           sprintf (message,  "No p-values in input .1D data file: %s", 
00563                    FDR_input1D_filename);
00564           FDR_error (message);
00565         }
00566 
00567       FDR_nxyz = nxyz;
00568       FDR_nthr = nxyz;
00569     }
00570   
00571   else
00572     {
00573       /*----- Open the input 3D dataset -----*/
00574       if (!FDR_quiet)  printf ("Reading input dataset: %s \n", 
00575                                FDR_input_filename);
00576       FDR_dset = THD_open_dataset(FDR_input_filename);
00577       
00578       if (FDR_dset == NULL)
00579         {
00580           sprintf (message, "Cannot open input dataset %s", 
00581                    FDR_input_filename); 
00582           FDR_error (message);
00583         }
00584       
00585       /*----- Get dimensions of input dataset -----*/
00586       nx   = DSET_NX(FDR_dset);   
00587       ny   = DSET_NY(FDR_dset);   
00588       nz   = DSET_NZ(FDR_dset);
00589       nxyz = nx*ny*nz;
00590       
00591       
00592       /*----- Check for compatible dimensions -----*/
00593       if (FDR_mask != NULL)
00594         {
00595           if ((nx != mx) || (ny != my) || (nz != mz))
00596             FDR_error ("Mask and input dataset have incompatible dimensions");
00597           FDR_nxyz = nxyz;
00598           FDR_nthr = nthr;
00599         }
00600       else
00601         {
00602           FDR_nxyz = nxyz;
00603           FDR_nthr = nxyz;
00604         }
00605 
00606 
00607       /*----- Check whether output dataset already exists -----*/
00608       check_one_output_file (FDR_dset, FDR_output_prefix);
00609     }
00610 
00611 
00612   /*----- Initialize constant c(N) -----*/
00613   if (FDR_cn < 0.0)
00614     {
00615       double cn;
00616       cn = 0.0;
00617       for (ixyz = 1;  ixyz <= FDR_nthr;  ixyz++)
00618         cn += 1.0 / ixyz;
00619       FDR_cn = cn;
00620       if (!FDR_quiet)
00621         printf ("c(N) = %f \n", FDR_cn);
00622     }
00623   
00624   /*----- Initialize voxel pointers -----*/
00625   for (ibin = 0;  ibin < FDR_MAX_LL;  ibin++)
00626     FDR_head_voxel[ibin] = NULL;
00627 
00628 
00629 }
00630 
00631 
00632 /*---------------------------------------------------------------------------*/
00633 /*
00634   Create an empty voxel.
00635 */
00636   
00637 voxel * create_voxel ()
00638 {
00639   voxel * voxel_ptr = NULL;
00640 
00641   voxel_ptr = (voxel *) malloc (sizeof(voxel));
00642   MTEST (voxel_ptr);
00643   
00644   voxel_ptr->ixyz = 0;
00645   voxel_ptr->pvalue = 0.0;
00646   voxel_ptr->next_voxel = NULL;
00647 
00648   return (voxel_ptr);
00649   
00650 }
00651 
00652 
00653 /*---------------------------------------------------------------------------*/
00654 /*
00655   Add a new voxel to the linked list of voxels.
00656 */
00657 
00658 voxel * add_voxel (voxel * new_voxel, voxel * head_voxel)
00659 {
00660   voxel * voxel_ptr = NULL;
00661 
00662   if ((head_voxel == NULL) || (new_voxel->pvalue >= head_voxel->pvalue))
00663     {
00664       new_voxel->next_voxel = head_voxel;
00665       head_voxel = new_voxel;
00666     }
00667 
00668   else
00669     {
00670       voxel_ptr = head_voxel;
00671 
00672       while ((voxel_ptr->next_voxel != NULL) && 
00673              (new_voxel->pvalue < voxel_ptr->next_voxel->pvalue))
00674         voxel_ptr = voxel_ptr->next_voxel;
00675       
00676       new_voxel->next_voxel = voxel_ptr->next_voxel;
00677       voxel_ptr->next_voxel = new_voxel;
00678     }
00679 
00680   return (head_voxel);
00681 }
00682 
00683 
00684 /*---------------------------------------------------------------------------*/
00685 /*
00686   Create and initialize a new voxel, and add to list of voxels.
00687 */
00688 
00689 voxel * new_voxel (int ixyz, float pvalue, voxel * head_voxel)
00690 
00691 {
00692   voxel * voxel_ptr = NULL;
00693 
00694   voxel_ptr = create_voxel ();
00695 
00696   voxel_ptr->ixyz      = ixyz;
00697   voxel_ptr->pvalue    = pvalue;
00698 
00699   head_voxel = add_voxel (voxel_ptr, head_voxel);
00700 
00701   return (head_voxel);
00702   
00703 }
00704 
00705 
00706 /*---------------------------------------------------------------------------*/
00707 /*
00708   Deallocate memory for all voxel lists.
00709 */
00710 
00711 void delete_all_voxels ()
00712 {
00713   int ibin;
00714   voxel * voxel_ptr  = NULL;     /* pointer to current voxel */
00715   voxel * next_voxel = NULL;     /* pointer to next voxel */
00716 
00717 
00718   for (ibin = 0;  ibin < FDR_MAX_LL;  ibin++) 
00719     {
00720       voxel_ptr = FDR_head_voxel[ibin];
00721       while (voxel_ptr != NULL)
00722         {
00723           next_voxel = voxel_ptr->next_voxel;
00724           free (voxel_ptr);
00725           voxel_ptr = next_voxel;
00726         }
00727       FDR_head_voxel[ibin] = NULL;
00728     }
00729   
00730 }
00731 
00732 
00733 /*---------------------------------------------------------------------------*/
00734 /*
00735   Save voxel contents of all voxels into float array (sub-brick).
00736 */
00737 
00738 void save_all_voxels (float * far)
00739 {
00740   int ixyz, ibin;
00741   voxel * voxel_ptr  = NULL;     /* pointer to voxel */
00742   
00743 
00744   /*----- Initialize all voxels to zero -----*/
00745   for (ixyz = 0;  ixyz < FDR_nxyz;  ixyz++)
00746     far[ixyz] = 0.0;
00747 
00748 
00749   for (ibin = 0;  ibin < FDR_MAX_LL;  ibin++) 
00750     {
00751       voxel_ptr = FDR_head_voxel[ibin];
00752   
00753       while (voxel_ptr != NULL)
00754         {
00755           far[voxel_ptr->ixyz] = voxel_ptr->pvalue;
00756           voxel_ptr = voxel_ptr->next_voxel;
00757         }
00758 
00759     }
00760 
00761 }
00762 
00763 
00764 /*---------------------------------------------------------------------------*/
00765 /*
00766   Calculate FDR z-scores for all voxels within one volume.
00767 */
00768 
00769 void process_volume (float * ffim, int statcode, float * stataux)
00770 
00771 {
00772   int ixyz;                      /* voxel index */
00773   int icount;                    /* count of sorted p-values */
00774   float fval;                    /* voxel input statistical value */
00775   float pval;                    /* voxel input stat. p-value */
00776   float qval;                    /* voxel FDR q-value */
00777   float zval;                    /* voxel FDR z-score */
00778   float qval_min;                /* smallest previous q-value */
00779   voxel * head_voxel = NULL;     /* linked list of voxels */
00780   voxel * voxel_ptr  = NULL;     /* pointer to current voxel */
00781   int ibin;                      /* p-value bin */
00782   int   * iarray = NULL;         /* output array of voxel indices */
00783   float * parray = NULL;         /* output array of voxel p-values */
00784   float * qarray = NULL;         /* output array of voxel FDR q-values */
00785   float * zarray = NULL;         /* output array of voxel FDR z-scores */
00786 
00787   
00788   /*----- Allocate memory for screen output arrays -----*/
00789   if (FDR_list)
00790     {
00791       iarray = (int   *) malloc (sizeof(int)   * FDR_nthr);   MTEST(iarray);
00792       parray = (float *) malloc (sizeof(float) * FDR_nthr);   MTEST(parray);
00793       qarray = (float *) malloc (sizeof(float) * FDR_nthr);   MTEST(qarray);
00794       zarray = (float *) malloc (sizeof(float) * FDR_nthr);   MTEST(zarray);
00795     }
00796   
00797   
00798   /*----- Loop over all voxels; sort p-values -----*/
00799   icount = FDR_nthr;
00800   for (ixyz = 0;  ixyz < FDR_nxyz;  ixyz++)
00801     {
00802 
00803       /*----- First, check if voxel is inside the mask -----*/
00804       if (FDR_mask != NULL)
00805         if (!FDR_mask[ixyz]) continue;
00806 
00807 
00808       /*----- Convert stats to p-values -----*/
00809       fval = fabs(ffim[ixyz]);
00810       if (statcode <= 0)
00811         pval = fval;
00812       else
00813         pval = THD_stat_to_pval (fval, statcode, stataux);
00814 
00815       if (pval >= 1.0)  
00816         {
00817           /*----- Count but don't sort voxels with p-value = 1 -----*/
00818           icount--;
00819           if (FDR_list)
00820             {
00821               iarray[icount] = ixyz;
00822               parray[icount] = 1.0;
00823               qarray[icount] = 1.0;
00824               zarray[icount] = 0.0;
00825             }
00826         }
00827       else
00828         { 
00829           /*----- Place voxel in p-value bin -----*/
00830           ibin = (int)  (pval * (FDR_MAX_LL));
00831           if (ibin < 0)  ibin = 0;
00832           if (ibin > FDR_MAX_LL-1)  ibin = FDR_MAX_LL-1;
00833           head_voxel = new_voxel (ixyz, pval, FDR_head_voxel[ibin]);
00834           FDR_head_voxel[ibin] = head_voxel;
00835         }
00836     }
00837 
00838   
00839   /*----- Calculate FDR q-values -----*/
00840   qval_min = 1.0;
00841   ibin = FDR_MAX_LL-1;
00842   while (ibin >= 0) 
00843     {
00844       voxel_ptr = FDR_head_voxel[ibin];
00845   
00846       while (voxel_ptr != NULL)
00847         {
00848           /*----- Convert sorted p-values to FDR q-values -----*/
00849           pval = voxel_ptr->pvalue;
00850           qval = FDR_cn * (pval*FDR_nthr) / icount;
00851           if (qval > qval_min)
00852             qval = qval_min;
00853           else
00854             qval_min = qval;
00855 
00856           /*----- Convert FDR q-value to FDR z-score -----*/
00857           if (qval < 1.0e-20)
00858             zval = 10.0;
00859           else
00860             zval = normal_p2t(qval);
00861 
00862           icount--;
00863 
00864           /*----- Save calculated values -----*/
00865           if (FDR_list)
00866             {
00867               iarray[icount] = voxel_ptr->ixyz;
00868               parray[icount] = pval;
00869               qarray[icount] = qval;
00870               zarray[icount] = zval;
00871             }
00872 
00873           voxel_ptr->pvalue = zval;
00874           voxel_ptr = voxel_ptr->next_voxel;
00875         }
00876 
00877       ibin--;
00878     }
00879 
00880 
00881   /*----- Write out the calculated values -----*/
00882   if (FDR_list)
00883     {
00884       printf ("%12s %12s %12s %12s \n", 
00885               "Index", "p-value", "q-value", "z-score");
00886       for (icount = 0;  icount < FDR_nthr;  icount++)
00887         {
00888           if (FDR_input1D_filename != NULL)
00889             ixyz = iarray[icount] + 1;
00890           else
00891             ixyz = iarray[icount];
00892           printf ("%12d %12.6f %12.6f %12.6f \n",  
00893                   ixyz, parray[icount], qarray[icount], zarray[icount]);
00894         }
00895 
00896       /*----- Deallocate memory for output arrays -----*/
00897       free (iarray);   free (parray);   free (qarray);   free (zarray);
00898     }
00899 
00900 
00901   /*----- Place FDR z-scores into float array -----*/
00902   save_all_voxels (ffim);
00903 
00904 
00905   /*----- Deallocate linked-list memory -----*/
00906   delete_all_voxels();
00907 
00908 }
00909 
00910 
00911 /*---------------------------------------------------------------------------*/
00912 /*
00913   Perform all processing for this array of p-values.
00914 */
00915 
00916 void process_1ddata ()
00917 
00918 {
00919   float * ffim = NULL;     /* input list of p-values */
00920 
00921   
00922   /*----- Set pointer to input array of p-values -----*/
00923   ffim = FDR_input1D_data;
00924 
00925 
00926   /*----- Calculate FDR z-scores for all input p-values  -----*/
00927   process_volume (ffim, -1, NULL);
00928 
00929   if( FDR_output_prefix != NULL && ffim != NULL ){
00930     MRI_IMAGE *im = mri_new_vol_empty( FDR_nxyz,1,1 , MRI_float ) ;
00931     mri_fix_data_pointer( ffim , im ) ;
00932     mri_write_1D( FDR_output_prefix , im ) ;
00933     mri_fix_data_pointer( NULL , im ) ;
00934     mri_free( im ) ;
00935   }
00936 
00937   /*----- Deallocate memory -----*/
00938   if (ffim != NULL) { free (ffim);   ffim = NULL; }
00939 
00940 }
00941 
00942 
00943 /*---------------------------------------------------------------------------*/
00944 /*
00945   Convert one volume to another type, autoscaling:
00946      nxy   = # voxels
00947      itype = input datum type
00948      ivol  = pointer to input volume
00949      otype = output datum type
00950      ovol  = pointer to output volume (again, must be pre-malloc-ed)
00951   Return value is the scaling factor used (0.0 --> no scaling).
00952 */
00953 
00954 float EDIT_coerce_autoscale_new( int nxyz ,
00955                                  int itype,void *ivol , int otype,void *ovol )
00956 {
00957   float fac=0.0 , top ;
00958   
00959   if( MRI_IS_INT_TYPE(otype) ){
00960     top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
00961     if (top == 0.0)  fac = 0.0;
00962     else  fac = MRI_TYPE_maxval[otype]/top;
00963   }
00964   
00965   EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
00966   return ( fac );
00967 }
00968 
00969 
00970 /*---------------------------------------------------------------------------*/
00971 /*
00972   Perform all processing for this sub-brick.
00973 */
00974 
00975 void process_subbrick (THD_3dim_dataset * dset, int ibrick)
00976 
00977 {
00978   const float EPSILON = 1.0e-10;
00979   float factor;            /* factor is new scale factor for this sub-brick */
00980   void * vfim = NULL;      /* sub-brick data pointer */
00981   float * ffim = NULL;     /* sub-brick data in floating point format */
00982   char brick_label[THD_MAX_NAME];       /* sub-brick label */
00983 
00984 
00985   if (!FDR_quiet)  printf ("Processing sub-brick #%d \n", ibrick);
00986 
00987   
00988   /*----- Allocate memory for float data -----*/
00989   ffim = (float *) malloc (sizeof(float) * FDR_nxyz);   MTEST (ffim);
00990 
00991 
00992   /*----- Convert sub-brick to float stats -----*/
00993   SUB_POINTER (dset, ibrick, 0, vfim);
00994   EDIT_coerce_scale_type (FDR_nxyz, DSET_BRICK_FACTOR(dset,ibrick),
00995                           DSET_BRICK_TYPE(dset,ibrick), vfim,   /* input  */
00996                           MRI_float                   , ffim);  /* output */
00997 
00998 
00999   /*----- Calculate FDR z-scores for all voxels within this volume -----*/
01000   process_volume (ffim, DSET_BRICK_STATCODE(dset,ibrick),
01001                         DSET_BRICK_STATAUX (dset,ibrick));
01002 
01003 
01004   /*----- Replace old sub-brick with new z-scores -----*/
01005   SUB_POINTER (dset, ibrick, 0, vfim);
01006   factor = EDIT_coerce_autoscale_new (FDR_nxyz, MRI_float, ffim,
01007                                       DSET_BRICK_TYPE(dset,ibrick), vfim);  
01008   if (factor < EPSILON)  factor = 0.0;
01009   else factor = 1.0 / factor;
01010   
01011 
01012   /*----- edit the sub-brick -----*/
01013   strcpy (brick_label, "FDRz ");
01014   strcat (brick_label, DSET_BRICK_LABEL(dset, ibrick));
01015   EDIT_BRICK_LABEL (dset, ibrick, brick_label);
01016   EDIT_BRICK_FACTOR (dset, ibrick, factor);
01017   EDIT_BRICK_TO_FIZT (dset, ibrick);
01018  
01019 
01020   /*----- Deallocate memory -----*/
01021   if (ffim != NULL) { free (ffim);   ffim = NULL; }
01022 
01023 }
01024 
01025 
01026 /*---------------------------------------------------------------------------*/
01027 /*
01028   Process the input dataset.
01029 */
01030 
01031 THD_3dim_dataset * process_dataset ()
01032 
01033 {
01034   THD_3dim_dataset * new_dset = NULL;     /* output bucket dataset */
01035   int ibrick, nbricks;                    /* sub-brick indices */
01036   int statcode;                           /* type of stat. sub-brick */
01037 
01038 
01039   /*----- Make full copy of input dataset -----*/
01040   new_dset = EDIT_full_copy(FDR_dset, FDR_output_prefix);
01041 
01042 
01043   /*----- Record history of dataset -----*/
01044   tross_Copy_History( FDR_dset , new_dset ) ;
01045 
01046   if( commandline != NULL )
01047     {
01048       tross_Append_History ( new_dset, commandline);
01049       free(commandline) ;
01050     }
01051 
01052 
01053   /*----- Deallocate memory for input dataset -----*/   
01054   THD_delete_3dim_dataset (FDR_dset , False );  FDR_dset = NULL ;
01055 
01056 
01057   /*----- Loop over sub-bricks in the dataset -----*/
01058   nbricks = DSET_NVALS(new_dset);
01059   for (ibrick = 0;  ibrick < nbricks;  ibrick++)
01060     {
01061       statcode = DSET_BRICK_STATCODE(new_dset, ibrick);
01062       if (FUNC_IS_STAT(statcode))
01063         {
01064           /*----- Process the statistical sub-bricks -----*/
01065           if (!FDR_quiet)  
01066             printf ("ibrick = %3d   statcode = %5s \n", 
01067                     ibrick, FUNC_prefixstr[statcode]);
01068           process_subbrick (new_dset, ibrick);
01069         }
01070     }
01071 
01072 
01073   return (new_dset);
01074 }
01075 
01076 
01077 /*---------------------------------------------------------------------------*/
01078 /*
01079   Output the results as a bucket dataset.
01080 */
01081 
01082 void output_results (THD_3dim_dataset * new_dset)
01083 {
01084   int ierror;     /* flag for errors in editing dataset */
01085 
01086 
01087   /*----- Make sure that output is a bucket dataset -----*/
01088   ierror = EDIT_dset_items( new_dset ,
01089                             ADN_func_type , FUNC_BUCK_TYPE,
01090                             ADN_none ) ;
01091   if (ierror > 0)  
01092     FDR_error ("Errors in attempting to create output dataset.");
01093 
01094 
01095   /*----- Output the FDR dataset -----*/
01096   if( !FDR_quiet ) fprintf(stderr,"Computing sub-brick statistics\n") ;
01097   THD_load_statistics( new_dset ) ;
01098 
01099   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
01100   if( !FDR_quiet ) fprintf(stderr,"Wrote output to %s\n", DSET_BRIKNAME(new_dset) );
01101   
01102 
01103   /*----- Deallocate memory for output dataset -----*/   
01104   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
01105   
01106 }
01107 
01108 
01109 /*---------------------------------------------------------------------------*/
01110 
01111 int main( int argc , char * argv[] )
01112 
01113 {
01114   THD_3dim_dataset * new_dset = NULL;      /* output bucket dataset */
01115   
01116 
01117   /*----- Identify software -----*/
01118   printf ("\n\n");
01119   printf ("Program:          %s \n", PROGRAM_NAME);
01120   printf ("Author:           %s \n", PROGRAM_AUTHOR); 
01121   printf ("Initial Release:  %s \n", PROGRAM_INITIAL);
01122   printf ("Latest Revision:  %s \n", PROGRAM_LATEST);
01123   printf ("\n");
01124 
01125 
01126   /*----- Initialize program:  get all operator inputs; 
01127     create mask for voxels above mask threshold -----*/
01128   initialize_program (argc, argv);
01129 
01130 
01131   if (FDR_input1D_filename != NULL)
01132     {
01133       /*----- Process list of p-values -----*/
01134       process_1ddata ();
01135     }
01136   else
01137     {
01138       /*----- Process 3d dataset -----*/
01139       new_dset = process_dataset ();
01140 
01141       /*----- Output the results as a bucket dataset -----*/
01142       output_results (new_dset);
01143     }
01144   
01145   exit(0) ;
01146 }
01147 
01148 
01149 /*---------------------------------------------------------------------------*/
01150 
01151