00001 /*****************************************************************************
00002    Major portions of this software are copyrighted by the Medical College
00003    of Wisconsin, 1994-2000, and are released under the Gnu General Public
00004    License, Version 2.  See the file README.Copyright for details.
00005 ******************************************************************************/
00006 
00007 /*
00008   This file contains routines common to the nonparametric statistical
00009   analysis programs.
00010 
00011   File:    NPstats.c
00012   Author:  B. Douglas Ward
00013   Date:    23 July 1997
00014 
00015   Mod:     Added changes for incorporating History notes.
00016   Date:    08 September 1999
00017 
00018   Mod:     Replaced dataset input code with calls to THD_open_dataset,
00019            to allow operator selection of individual sub-bricks for input.
00020   Date:    03 December 1999
00021 
00022   Mod:     Moved routines for sorting numbers and determining ranks to
00023            separate file ranks.c.
00024   Date:    31 March 2000
00025 
00026   Mod:     Modified routines write_afni_fizt and write_afni_fict so that 
00027            all output subbricks will now have the scaled short integer format.
00028   Date:    14 March 2002
00029 */
00030 
00031 
00032 /*---------------------------------------------------------------------------*/
00033 /*
00034    Print error message and stop.
00035 */
00036 
00037 void NP_error (char * message)
00038 {
00039    fprintf (stderr, "%s Error: %s \n", PROGRAM_NAME, message);
00040    exit(1);
00041 }
00042 
00043 
00044 /*---------------------------------------------------------------------------*/
00045 
00046 /** macro to test a malloc-ed pointer for validity **/
00047 
00048 #define MTEST(ptr) \
00049 if((ptr)==NULL) \
00050 ( NP_error ("Cannot allocate memory") )
00051      
00052 
00053 /*---------------------------------------------------------------------------*/
00054 /*
00055   Include routines for sorting numbers and determining ranks.
00056 */
00057 
00058 #include "ranks.c"
00059 
00060 
00061 /*---------------------------------------------------------------------------*/
00062 /*
00063    Get the dimensions of the 3d AFNI data sets.
00064 */
00065 
00066 void get_dimensions (NP_options * option_data)
00067 {
00068   
00069    THD_3dim_dataset * dset=NULL;
00070 
00071    /*----- read first dataset to get dimensions, etc. -----*/
00072 
00073    dset = THD_open_dataset( option_data->first_dataset ) ;
00074    if( ! ISVALID_3DIM_DATASET(dset) ){
00075       fprintf(stderr,"*** Unable to open dataset file %s\n", 
00076               option_data->first_dataset);
00077       exit(1) ;
00078    }
00079 
00080    /*----- data set dimensions in voxels -----*/
00081    option_data->nx = dset->daxes->nxx ;
00082    option_data->ny = dset->daxes->nyy ;
00083    option_data->nz = dset->daxes->nzz ;       
00084    option_data->nxyz = option_data->nx * option_data->ny * option_data->nz ;
00085 
00086    THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00087 
00088 }
00089 
00090 
00091 /*---------------------------------------------------------------------------*/
00092 /*
00093   Check whether one output file already exists.
00094 */
00095 
00096 void check_one_output_file (NP_options * option_data, char * filename)
00097 {
00098   THD_3dim_dataset * dset=NULL;       /* input afni data set pointer */
00099   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00100   int ierror;                         /* number of errors in editing data */
00101   
00102   
00103   /*----- read first dataset -----*/
00104   dset = THD_open_dataset (option_data->first_dataset ) ;
00105   if( ! ISVALID_3DIM_DATASET(dset) ){
00106     fprintf(stderr,"*** Unable to open dataset file %s\n",
00107             option_data->first_dataset);
00108     exit(1) ;
00109   }
00110   
00111   /*-- make an empty copy of this dataset, for eventual output --*/
00112   new_dset = EDIT_empty_copy( dset ) ;
00113   
00114   
00115   ierror = EDIT_dset_items( new_dset ,
00116                             ADN_prefix , filename ,
00117                             ADN_label1 , filename ,
00118                             ADN_directory_name , option_data->session ,
00119                             ADN_self_name , filename ,
00120                             ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : 
00121                                                       GEN_FUNC_TYPE ,
00122                             ADN_none ) ;
00123   
00124   if( ierror > 0 ){
00125     fprintf(stderr,
00126             "*** %d errors in attempting to create output dataset!\n", ierror ) ;
00127     exit(1) ;
00128   }
00129   
00130   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
00131     fprintf(stderr,
00132             "*** Output dataset file %s already exists--cannot continue!\a\n",
00133             new_dset->dblk->diskptr->header_name ) ;
00134     exit(1) ;
00135   }
00136   
00137   /*----- deallocate memory -----*/   
00138   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00139   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00140   
00141 }
00142 
00143 
00144 /*---------------------------------------------------------------------------*/
00145 
00146 /** macro to open a dataset and make it ready for processing **/
00147 
00148 #define DOPEN(ds,name)                                                        \
00149 do{ int pv ; (ds) = THD_open_dataset((name)) ;                                \
00150        if( !ISVALID_3DIM_DATASET((ds)) ){                                     \
00151           fprintf(stderr,"*** Can't open dataset: %s\n",(name)) ; exit(1) ; } \
00152        if( (ds)->daxes->nxx!=nx || (ds)->daxes->nyy!=ny ||                    \
00153            (ds)->daxes->nzz!=nz ){                                            \
00154           fprintf(stderr,"*** Axes mismatch: %s\n",(name)) ; exit(1) ; }      \
00155        if( DSET_NVALS((ds)) != 1 ){                                           \
00156          fprintf(stderr,"*** Must specify 1 sub-brick: %s\n",(name));exit(1);}\
00157        THD_load_datablock( (ds)->dblk ) ;                                     \
00158        pv = DSET_PRINCIPAL_VALUE((ds)) ;                                      \
00159        if( DSET_ARRAY((ds),pv) == NULL ){                                     \
00160           fprintf(stderr,"*** Can't access data: %s\n",(name)) ; exit(1); }   \
00161        if( DSET_BRICK_TYPE((ds),pv) == MRI_complex ){                         \
00162           fprintf(stderr,"*** Can't use complex data: %s\n",(name)) ; exit(1);\
00163        }                                                                      \
00164        break ; } while (0)
00165 
00166 
00167 /*---------------------------------------------------------------------------*/
00168 
00169 /** macro to return pointer to correct location in brick for current processing **/
00170 
00171 #define SUB_POINTER(ds,vv,ind,ptr)                                            \
00172    do{ switch( DSET_BRICK_TYPE((ds),(vv)) ){                                  \
00173          default: fprintf(stderr,"\n*** Illegal datum! ***\n");exit(1);       \
00174             case MRI_short:{ short * fim = (short *) DSET_ARRAY((ds),(vv)) ;  \
00175                             (ptr) = (void *)( fim + (ind) ) ;                 \
00176             } break ;                                                         \
00177             case MRI_byte:{ byte * fim = (byte *) DSET_ARRAY((ds),(vv)) ;     \
00178                             (ptr) = (void *)( fim + (ind) ) ;                 \
00179             } break ;                                                         \
00180             case MRI_float:{ float * fim = (float *) DSET_ARRAY((ds),(vv)) ;  \
00181                              (ptr) = (void *)( fim + (ind) ) ;                \
00182             } break ; } break ; } while(0)
00183 
00184 
00185 /*---------------------------------------------------------------------------*/
00186 /*
00187   Read one AFNI data set from file 'filename'. 
00188   The data is converted to floating point (in ffim).
00189 */
00190 
00191 void read_afni_data (NP_options * option_data, char * filename, 
00192                      int piece_len, int fim_offset,  float * ffim)
00193 {
00194   int iv;                          /* index number of intensity sub-brick */
00195   THD_3dim_dataset * dset=NULL;    /* data set pointer */
00196   void * vfim = NULL;              /* image data pointer */
00197   int nx, ny, nz, nxyz;            /* data set dimensions in voxels */
00198   
00199   nx = option_data->nx;
00200   ny = option_data->ny;
00201   nz = option_data->nz;
00202   nxyz = option_data->nxyz;
00203   
00204     
00205   /*----- read in the data -----*/
00206   DOPEN (dset, filename) ;
00207   iv = DSET_PRINCIPAL_VALUE(dset) ;
00208   
00209   /*----- convert it to floats (in ffim) -----*/
00210   SUB_POINTER (dset, iv, fim_offset, vfim) ;
00211   EDIT_coerce_scale_type (piece_len, DSET_BRICK_FACTOR(dset,iv),
00212                           DSET_BRICK_TYPE(dset,iv), vfim,      /* input  */
00213                           MRI_float               ,ffim  ) ;   /* output */
00214   
00215   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00216 }
00217 
00218 
00219 /*---------------------------------------------------------------------------*/
00220 /*
00221   Convert one volume to another type, autoscaling:
00222      nxy   = # voxels
00223      itype = input datum type
00224      ivol  = pointer to input volume
00225      otype = output datum type
00226      ovol  = pointer to output volume (again, must be pre-malloc-ed)
00227   Return value is the scaling factor used (0.0 --> no scaling).
00228 */
00229 
00230 float EDIT_coerce_autoscale_new( int nxyz ,
00231                                  int itype,void *ivol , int otype,void *ovol )
00232 {
00233   float fac=0.0 , top ;
00234   
00235   if( MRI_IS_INT_TYPE(otype) ){
00236     top = MCW_vol_amax( nxyz,1,1 , itype,ivol ) ;
00237     if (top == 0.0)  fac = 0.0;
00238     else  fac = MRI_TYPE_maxval[otype]/top;
00239   }
00240   
00241   EDIT_coerce_scale_type( nxyz , fac , itype,ivol , otype,ovol ) ;
00242   return ( fac );
00243 }
00244 
00245 
00246 /*---------------------------------------------------------------------------*/
00247 /*
00248   Routine to write one AFNI fizt data set.
00249   
00250   This data set is 'fizt' type (intensity + Standard Normal statistic)
00251   
00252   The intensity data is in ffim, and the corresponding statistic is in ftr.
00253   
00254 */
00255 
00256 void write_afni_fizt (int argc, char ** argv, NP_options * option_data,  
00257                       char * filename, float * ffim,  float * ftr)
00258 {
00259   int nxyz;                           /* number of voxels */
00260   int ii;                             /* voxel index */
00261   THD_3dim_dataset * dset=NULL;       /* input afni data set pointer */
00262   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00263   int ierror;                         /* number of errors in editing data */
00264   int ibuf[32];                       /* integer buffer */
00265   float fbuf[MAX_STAT_AUX];           /* float buffer */
00266   float fimfac;                       /* scale factor for short data */
00267   int output_datum;                   /* data type for output data */
00268   short * tsp;                        /* 2nd sub-brick data pointer */
00269   void  * vdif;                       /* 1st sub-brick data pointer */
00270   int func_type;                      /* afni data set type */
00271   float top, func_scale_short;        /* parameters for scaling data */
00272   
00273   
00274   /*----- initialize local variables -----*/
00275   nxyz = option_data->nxyz;
00276   
00277   /*----- read first dataset -----*/
00278   dset = THD_open_dataset (option_data->first_dataset) ;
00279   if( ! ISVALID_3DIM_DATASET(dset) ){
00280     fprintf(stderr,"*** Unable to open dataset file %s\n",
00281             option_data->first_dataset);
00282     exit(1) ;
00283   }
00284   
00285   /*-- make an empty copy of this dataset, for eventual output --*/
00286   new_dset = EDIT_empty_copy( dset ) ;
00287   
00288   
00289   output_datum = MRI_short ;
00290   ibuf[0] = output_datum ; ibuf[1] = MRI_short ;
00291 
00292   
00293   func_type = FUNC_ZT_TYPE;
00294   
00295   ierror = EDIT_dset_items( new_dset ,
00296                             ADN_prefix , filename ,
00297                             ADN_label1 , filename ,
00298                             ADN_directory_name , option_data->session ,
00299                             ADN_self_name , filename ,
00300                             ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : 
00301                                                       GEN_FUNC_TYPE ,
00302                             ADN_func_type , func_type ,
00303                             ADN_nvals , FUNC_nvals[func_type] ,
00304                             ADN_datum_array , ibuf ,
00305                             ADN_malloc_type, DATABLOCK_MEM_MALLOC ,  
00306                             ADN_none ) ;
00307   
00308   if( ierror > 0 ){
00309     fprintf(stderr,
00310           "*** %d errors in attempting to create output dataset!\n", ierror ) ;
00311     exit(1) ;
00312   }
00313   
00314   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
00315     fprintf(stderr,
00316             "*** Output dataset file %s already exists--cannot continue!\a\n",
00317             new_dset->dblk->diskptr->header_name ) ;
00318     exit(1) ;
00319   }
00320   
00321   /*----- deleting exemplar dataset -----*/ 
00322   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00323   
00324   
00325   /*----- allocate memory for output data -----*/
00326   vdif = (void *)  malloc( mri_datum_size(output_datum) * nxyz ) ;
00327   tsp  = (short *) malloc( sizeof(short) * nxyz )                ;
00328   
00329   /*----- attach bricks to new data set -----*/
00330   mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0)); 
00331   mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1));  
00332   
00333   
00334   /*----- convert data type to output specification -----*/
00335   fimfac = EDIT_coerce_autoscale_new (nxyz, 
00336                                       MRI_float, ffim, 
00337                                       output_datum, vdif);
00338   if (fimfac != 0.0)  fimfac = 1.0 / fimfac;
00339   
00340 #define TOP_SS  32700
00341   
00342   top = TOP_SS/FUNC_ZT_SCALE_SHORT;
00343   func_scale_short = FUNC_ZT_SCALE_SHORT;
00344   
00345   for (ii = 0;  ii < nxyz;  ii++)
00346     {
00347       if (ftr[ii] > top)
00348         tsp[ii] = TOP_SS;
00349       else  if (ftr[ii] < -top)
00350         tsp[ii] = -TOP_SS;
00351       else  if (ftr[ii] >= 0.0)
00352         tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5);
00353       else
00354         tsp[ii] = (short) (func_scale_short * ftr[ii] - 0.5);
00355     }
00356   
00357 
00358   /*----- write afni fizt data set -----*/
00359   printf("--- Writing 'fizt' dataset into %s\n",DSET_BRIKNAME(new_dset)) ;
00360   
00361   for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
00362   (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
00363   
00364   fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
00365   fbuf[1] = 1.0 / func_scale_short ;
00366   (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
00367 
00368 
00369   /*----- Record history of dataset -----*/  
00370   tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
00371   
00372   THD_load_statistics( new_dset ) ;
00373   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
00374   
00375   /*----- deallocate memory -----*/   
00376   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00377   
00378 }
00379 
00380 
00381 /*---------------------------------------------------------------------------*/
00382 /*
00383   Routine to write one AFNI fict data set.
00384   
00385   This data set is 'fict' type (intensity + Chi-Squared statistic)
00386   
00387   The intensity data is in ffim, and the corresponding statistic is in ftr.
00388   
00389   dof = degrees of freedom
00390 
00391 */
00392 
00393 void write_afni_fict (int argc, char ** argv, NP_options * option_data,  
00394                       char * filename, float * ffim,  float * ftr,  int dof)
00395 {
00396   int nxyz;                           /* number of voxels */
00397   int ii;                             /* voxel index */
00398   THD_3dim_dataset * dset=NULL;       /* input afni data set pointer */
00399   THD_3dim_dataset * new_dset=NULL;   /* output afni data set pointer */
00400   int ierror;                         /* number of errors in editing data */
00401   int ibuf[32];                       /* integer buffer */
00402   float fbuf[MAX_STAT_AUX];           /* float buffer */
00403   float fimfac;                       /* scale factor for short data */
00404   int output_datum;                   /* data type for output data */
00405   short * tsp;                        /* 2nd sub-brick data pointer */
00406   void  * vdif;                       /* 1st sub-brick data pointer */
00407   int func_type;                      /* afni data set type */
00408   float top, func_scale_short;        /* parameters for scaling data */
00409   
00410   
00411   /*----- initialize local variables -----*/
00412   nxyz = option_data->nxyz;
00413   
00414   /*----- read first dataset -----*/
00415   dset = THD_open_dataset (option_data->first_dataset) ;
00416   if( ! ISVALID_3DIM_DATASET(dset) ){
00417     fprintf(stderr,"*** Unable to open dataset file %s\n",
00418             option_data->first_dataset);
00419     exit(1) ;
00420   }
00421   
00422   /*-- make an empty copy of this dataset, for eventual output --*/
00423   new_dset = EDIT_empty_copy( dset ) ;
00424   
00425   
00426   output_datum = MRI_short ;
00427   ibuf[0] = output_datum ; ibuf[1] = MRI_short ;
00428 
00429   
00430   func_type = FUNC_CT_TYPE;
00431   
00432   ierror = EDIT_dset_items( new_dset ,
00433                             ADN_prefix , filename ,
00434                             ADN_label1 , filename ,
00435                             ADN_directory_name , option_data->session ,
00436                             ADN_self_name , filename ,
00437                             ADN_type , ISHEAD(dset) ? HEAD_FUNC_TYPE : 
00438                                                       GEN_FUNC_TYPE ,
00439                             ADN_func_type , func_type ,
00440                             ADN_nvals , FUNC_nvals[func_type] ,
00441                             ADN_datum_array , ibuf ,
00442                             ADN_malloc_type, DATABLOCK_MEM_MALLOC ,  
00443                             ADN_none ) ;
00444   
00445   if( ierror > 0 ){
00446     fprintf(stderr,
00447           "*** %d errors in attempting to create output dataset!\n", ierror ) ;
00448     exit(1) ;
00449   }
00450   
00451   if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
00452     fprintf(stderr,
00453             "*** Output dataset file %s already exists--cannot continue!\a\n",
00454             new_dset->dblk->diskptr->header_name ) ;
00455     exit(1) ;
00456   }
00457   
00458   /*----- deleting exemplar dataset -----*/ 
00459   THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
00460   
00461   
00462   /*----- allocate memory for output data -----*/
00463   vdif = (void *)  malloc( mri_datum_size(output_datum) * nxyz ) ;
00464   tsp  = (short *) malloc( sizeof(short) * nxyz )                ;
00465   
00466   /*----- attach bricks to new data set -----*/
00467   mri_fix_data_pointer (vdif, DSET_BRICK(new_dset,0)); 
00468   mri_fix_data_pointer (tsp, DSET_BRICK(new_dset,1));  
00469   
00470   
00471   /*----- convert data type to output specification -----*/
00472   fimfac = EDIT_coerce_autoscale_new (nxyz, 
00473                                       MRI_float, ffim, 
00474                                       output_datum, vdif);
00475   if (fimfac != 0.0)  fimfac = 1.0 / fimfac;
00476   
00477 #define TOP_SS  32700
00478   
00479   top = TOP_SS/FUNC_CT_SCALE_SHORT;
00480   func_scale_short = FUNC_CT_SCALE_SHORT;
00481   
00482   for (ii = 0;  ii < nxyz;  ii++)
00483     {
00484       if (ftr[ii] > top)
00485         tsp[ii] = TOP_SS;
00486       else  if (ftr[ii] < -top)
00487         tsp[ii] = -TOP_SS;
00488       else  if (ftr[ii] >= 0.0)
00489         tsp[ii] = (short) (func_scale_short * ftr[ii] + 0.5);
00490       else
00491         tsp[ii] = (short) (func_scale_short * ftr[ii] - 0.5);
00492 
00493     }
00494   
00495 
00496   /*----- write afni fict data set -----*/
00497   printf("--- Writing 'fict' dataset into %s\n", DSET_BRIKNAME(new_dset) ) ;
00498   
00499   fbuf[0] = dof;   
00500   for( ii=1 ; ii < MAX_STAT_AUX ; ii++ ) fbuf[ii] = 0.0 ;
00501   (void) EDIT_dset_items( new_dset , ADN_stat_aux , fbuf , ADN_none ) ;
00502   
00503   fbuf[0] = (output_datum == MRI_short && fimfac != 1.0 ) ? fimfac : 0.0 ;
00504   fbuf[1] = 1.0 / func_scale_short ;
00505   (void) EDIT_dset_items( new_dset , ADN_brick_fac , fbuf , ADN_none ) ;
00506 
00507 
00508   /*----- Record history of dataset -----*/  
00509   tross_Make_History( PROGRAM_NAME , argc , argv , new_dset ) ;
00510 
00511   THD_load_statistics( new_dset ) ;
00512   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
00513   
00514   /*----- deallocate memory -----*/   
00515   THD_delete_3dim_dataset( new_dset , False ) ; new_dset = NULL ;
00516   
00517 }
00518 
00519 
00520