do{ int pv ; (ds) = THD_open_dataset((name)) ;                                \
       if( !ISVALID_3DIM_DATASET((ds)) ){                                     \
          fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
       THD_load_datablock( (ds)->dblk ) ;                                     \
       pv = DSET_PRINCIPAL_VALUE((ds)) ;                                      \
       if( DSET_ARRAY((ds),pv) == NULL ){                                     \
          fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); }   \
       if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){                         \
          fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1);\
       }     \
       break ; } while (0)

if((ptr)==NULL) \
( FDR_error ("Cannot allocate memory") )

do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){                                  \
         default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1);       \
            case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ;  \
                            (ptr) = (void *)( fim + (ind) ) ;                 \
            } break ;                                                         \
            case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ;     \
                            (ptr) = (void *)( fim + (ind) ) ;                 \
            } break ;                                                         \
            case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ;  \
                             (ptr) = (void *)( fim + (ind) ) ;                \
            } break ; } break ; } while(0)

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 }

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 }

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 }

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 }

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 }

00116 {
00117   fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00118   exit(1);
00119 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }

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 }