/*****************************************************************************
   Major portions of this software are copyrighted by the Medical College
   of Wisconsin, 1994-2000, and are released under the Gnu General Public
   License, Version 2.  See the file README.Copyright for details.
******************************************************************************/

/*---------------------------------------------------------------------------*/
/*
  This program performs editing and/or merging of 3D datasets.

  File:    3dmerge.c
  Author:  R. W. Cox

  Mod:     Changes to implement "-doall" command option.
  Author:  B. D. Ward
  Date:    04 February 1998

  Mod:     Changes to implement "-1dindex" and "-1tindex" options.
  Author:  RWCox
  Date:    12 Nov 1998

  Mod:     Allow use of THD_open_dataset instead of THD_open_one_dataset.
  Author:  RWCox
  Date:    21 Feb 2001

  Mod:     Check for loaded dataset in merge operations
  Author:  RWCox
  Date:    02 Nov 2001
-----------------------------------------------------------------------------*/

#define PROGRAM_NAME "3dmerge"              /* name of this program */
#define LAST_MOD_DATE "02 Nov 2001"         /* date of last program mod */

#include "mrilib.h"
#include "parser.h"    /* 09 Aug 2000 */

#define MAIN
#define MEGA  1048576  /* 2^20 */

#define ALLOW_SUBV  /* 21 Feb 2001 */

#ifdef ALLOW_SUBV
# define DSET_OPEN THD_open_dataset
#else
# define DSET_OPEN THD_open_one_dataset
#endif

/*-------------------------- global data --------------------------*/

/*** combination flags:  mean, mean of nonzeros, max, ... ***/

#define CFLAG_MEAN   1  /* the default */
#define CFLAG_NZMEAN 2
#define CFLAG_MMAX   3
#define CFLAG_COUNT  4
#define CFLAG_AMAX   5
#define CFLAG_SMAX   6
#define CFLAG_ORDER  7
#define CFLAG_FISHER 8  /* 08 Aug 1996 */

#define THFLAG_NONE  0    /* 29 Aug 1996: added the ability to merge */
#define THFLAG_FICO  71   /*  threshold data as well as intensities  */

#define TANH(z)   tanh(z)  /* for the Fisher transformations */
#define ATANH(z) atanh(z)

/*** variables to hold status from command line options ***/

static EDIT_options MRG_edopt ;
static int MRG_have_edopt = 0 ;

static int   MRG_hits_g       = 0 ;
static int   MRG_cflag_g      = CFLAG_MEAN ;
static int   MRG_keepthr      = 0 ;
static float MRG_clust_rmm_g  = 0.0 ;
static float MRG_clust_vmul_g = 0.0 ;
static int   MRG_datum        = ILLEGAL_TYPE ;
static int   MRG_thdatum      = ILLEGAL_TYPE ;
static int   MRG_be_quiet     = 0 ;
static int   MRG_cflag_gthr   = THFLAG_NONE ;  /* 29 Aug 1996 */
static int   MRG_doall        = 0;             /* 02 Feb 1998 */
static int   MRG_verbose      = 0;             /* 29 Jul 2003 */
static int   MRG_nozero       = 0;             /* 18 Jan 2017 */

static char  MRG_output_session[THD_MAX_NAME]   = "./" ;
static char  MRG_output_prefix [THD_MAX_PREFIX] = "mrg" ;

#if 0
static char  MRG_output_label  [THD_MAX_LABEL]  = "\0" ;
#endif

static int   MRG_ivfim = -1 ;
static int   MRG_ivthr = -1 ;

static int   MRG_nscale = 0 ; /* 15 Sep 2000 */

static int   MRG_rank   = 0 ; /* 13 Nov 2007: ZSS */

/*--------------------------- prototypes ---------------------------*/
int MRG_read_opts( int , char ** ) ;
void MRG_Syntax(void) ;

/*--------------------------------------------------------------------
   read the arguments, load the global MRG_ variables, and return
   the index of the argument that we stopped at (that is, the first
   input filename)
----------------------------------------------------------------------*/

#ifdef MERGE_DEBUG
#  define DUMP1 fprintf(stderr,"ARG: %s\n",argv[nopt])
#  define DUMP2 fprintf(stderr,"ARG: %s %s\n",argv[nopt],argv[nopt+1])
#  define DUMP3 fprintf(stderr,"ARG: %s %s %s\n",argv[nopt],argv[nopt+1],argv[nopt+2])
#else
#  define DUMP1
#  define DUMP2
#  define DUMP3
#endif

int MRG_read_opts( int argc , char * argv[] )
{
   int nopt = 1 ;
   float val ;
   int  ival ;

   INIT_EDOPT( &MRG_edopt ) ;

   while( nopt < argc && argv[nopt][0] == '-' ){

      /**** check editing options ****/

      ival = EDIT_check_argv( argc , argv , nopt , &MRG_edopt ) ;
      if( ival > 0 ){
         nopt += ival ; MRG_have_edopt = 1 ;
         continue ;
      }

      /*** 18 Jan 2017: -nozero ***/

      if( strcasecmp(argv[nopt],"-nozero") == 0 ){
        MRG_nozero = 1 ; nopt++ ; continue ;
      }

      /*** 29 Jul 2003: -verbose ***/

      if( strncmp(argv[nopt],"-verb",5) == 0 ){
        MRG_verbose = MRG_edopt.verbose = 1 ;
        nopt++ ; continue ;
      }

      /**** Nov 1998: -1tindex and -1dindex ****/

      if( strncmp(argv[nopt],"-1dindex",5) == 0 ){
         if( ++nopt >= argc ){
            fprintf(stderr,"need an argument after -1dindex!\n") ; exit(1) ;
         }
         MRG_ivfim = MRG_edopt.iv_fim = (int) strtod( argv[nopt++] , NULL ) ;
         continue ;
      }

      if( strncmp(argv[nopt],"-1tindex",5) == 0 ){
         if( ++nopt >= argc ){
            fprintf(stderr,"need an argument after -1tindex!\n") ; exit(1) ;
         }
         MRG_ivthr = MRG_edopt.iv_thr = (int) strtod( argv[nopt++] , NULL ) ;
         continue ;
      }

      if( strncmp(argv[nopt],"-1rank", 5) == 0) {
         MRG_rank = MRG_edopt.rank = 1;
         nopt++ ; continue;
      }


      /**** 09 Aug 2000: -1fmask dset ****/

      if( strncmp(argv[nopt],"-1fmask",6) == 0 ){
         THD_3dim_dataset *mset ; int nn ;

         if( MRG_edopt.nfmask > 0 ){
            fprintf(stderr,"Can't have 2 -fmask options!\n") ;
            exit(1) ;
         }
         if( ++nopt >= argc ){
            fprintf(stderr,"No argument after %s?\n",argv[nopt-1]) ;
            exit(1);
         }

         mset = THD_open_dataset( argv[nopt++] ) ;
         if( mset == NULL ){
            fprintf(stderr,"*** Can't open -fmask dataset\n") ;
            exit(1) ;
         }
         if( DSET_BRICK_TYPE(mset,0) == MRI_complex ){
            fprintf(stderr,"*** Can't deal with complex-valued -fmask dataset\n");
            exit(1) ;
         }

         MRG_edopt.fmask  = THD_makemask( mset , 0 , 666.0,-666.0 ) ;
         MRG_edopt.nfmask = DSET_NVOX(mset) ;

         nn = THD_countmask(MRG_edopt.nfmask,MRG_edopt.fmask) ;
         if( nn < 2 ){
            fprintf(stderr,"*** Too few (%d) nonzero voxels in -fmask!\n",nn) ;
            exit(1) ;
         }

         DSET_delete(mset) ;

         if( MRG_edopt.filter_opt == FCFLAG_AVER )
            MRG_edopt.filter_opt = FCFLAG_MEAN ;

         if( MRG_edopt.thrfilter_opt == FCFLAG_AVER )
            MRG_edopt.thrfilter_opt = FCFLAG_MEAN ;

         MRG_edopt.fmclip = 1 ; /* 11 Oct 2007 */
         continue ;
      }

      /**** 11 Oct 2007: -1fm_noclip ****/

      if( strncmp(argv[nopt],"-1fm_noclip",9) == 0 ){
        MRG_edopt.fmclip = 0 ; nopt++ ; continue ;
      }

      /**** 11 Sep 2000: -1filter_winsor rmm nw ****/

      if( strcmp(argv[nopt],"-1filter_winsor") == 0 ){
         int nwin ;

         nopt++ ;
         if( nopt+1 >= argc ){
            fprintf(stderr,"*** Need 2 arguments after -1filter_winsor\n") ;
            exit(1) ;
         }
         MRG_edopt.filter_rmm  = strtod( argv[nopt++] , NULL ) ;
         if( MRG_edopt.filter_rmm <= 0.0 ){
            fprintf(stderr,"*** Illegal rmm value after -1filter_winsor\n");
            exit(1) ;
         }
         nwin = (int) strtod( argv[nopt++] , NULL ) ;
         if( nwin <= 0 || nwin >= FCFLAG_ONE_STEP ){
            fprintf(stderr,"*** Illegal nw value after -1filter_winsor\n");
            exit(1) ;
         }
         MRG_edopt.filter_opt = FCFLAG_WINSOR + nwin ;
         MRG_have_edopt = 1 ; continue ;
      }

      /**** 11 Oct 2007: -1filter_blur fwhm (cheaply) ****/
      /**** 09 Aug 2000: -1filter_expr rmm expr ****/

      if( strncmp(argv[nopt],"-1filter_expr",13) == 0 ||
          strncmp(argv[nopt],"-1filter_blur",13) == 0   ){

         PARSER_code *pcode ; int hsym[26] , aa , naa , nee ;
         double fval ;
         int doblur=(strncmp(argv[nopt],"-1filter_blur",13) == 0) ;

         nopt++ ;
         if( doblur ){
           if( nopt >= argc ) ERROR_exit("Need 1 argument after -1filter_blur");
         } else if( nopt+1 > argc ){
           ERROR_exit("Need 2 arguments after -1filter_expr") ;
         }

         MRG_edopt.filter_opt = FCFLAG_EXPR;
         fval = strtod( argv[nopt++] , NULL ) ;  /* rmm or fwhm */
         if( !doblur ){
           MRG_edopt.filter_rmm = fval ;
           if( MRG_edopt.filter_rmm <= 0.0 )
             ERROR_exit("Illegal rmm value after -1filter_expr");
           MRG_edopt.fexpr = argv[nopt++] ;
         } else {                              /* Make up blurring expr */
           MRG_edopt.filter_rmm = 1.3*fval ;
           if( MRG_edopt.filter_rmm <= 0.0 )
             ERROR_exit("Illegal fwhm value after -1filter_blur");
           MRG_edopt.fexpr = (char *)malloc(128) ;
           fval = 1.0 / (0.36067 * fval * fval) ;
           sprintf( MRG_edopt.fexpr , "exp(-r*r*%g)" , fval ) ;
         }
         pcode = PARSER_generate_code( MRG_edopt.fexpr ) ;  /* compile */

         if( pcode == NULL )
           ERROR_exit("Illegal expr after -1filter_%s",(doblur)?"blur":"expr");

#undef   MASK
#undef   PREDEFINED_MASK
#define  MASK(x)          (1 << (x-'a'))
#define  PREDEFINED_MASK  ( MASK('r') | MASK('x') | MASK('y') | MASK('z')   \
                                      | MASK('i') | MASK('j') | MASK('k') )

         /* check if only legal symbols are used */

         PARSER_mark_symbols( pcode , hsym ) ;  /* find symbols used */

         for( nee=naa=aa=0 ; aa < 26 ; aa++ ){
           if( hsym[aa] ){
             if( ((1<<aa) & PREDEFINED_MASK) == 0 ){
               nee++ ;  /* count of illegal symbols */
               ERROR_message("Symbol %c in -1filter_expr is illegal",'a'+aa) ;
             } else {
               naa++ ;  /* count of legal symbols */
             }
           }
         }
         if( nee > 0 ) exit(1) ;  /* can't use this expression! */
         if( naa == 0 ){          /* no symbols?  check if identically 0 */
           double atoz[26] , val ;
           val = PARSER_evaluate_one( pcode , atoz ) ;
           if( val != 0.0 )
             WARNING_message("-1filter_expr is constant = %g",val) ;
           else
             ERROR_exit("-1filter_expr is identically zero") ;
         }

         free(pcode) ;  /* will recompile when it is used */
         MRG_have_edopt = 1 ; continue ;
      }

      /**** -quiet ****/

      if( strncmp(argv[nopt],"-quiet",6) == 0 ){
         MRG_be_quiet = 1 ;
         nopt++ ; continue ;
      }


      /**** -keepthr ****/

      if( strncmp(argv[nopt],"-keepthr",6) == 0 ){
         MRG_keepthr = 1 ;
         nopt++ ; continue ;
      }

      /**** -doall ****/   /* 02 Feb 1998 */

      if( strncmp(argv[nopt],"-doall",6) == 0 ){
         MRG_doall = 1 ;
         nopt++ ; continue ;
      }

      /**** -datum type ****/

      if( strncmp(argv[nopt],"-datum",6) == 0 ){
DUMP2 ;
         if( ++nopt >= argc ){
            fprintf(stderr,"need an argument after -datum!\n") ; exit(1) ;
         }

         if( strcmp(argv[nopt],"short") == 0 ){
            MRG_datum = MRI_short ;
         } else if( strcmp(argv[nopt],"float") == 0 ){
            MRG_datum = MRI_float ;
         } else if( strcmp(argv[nopt],"byte") == 0 ){
            MRG_datum = MRI_byte ;
         } else if( strcmp(argv[nopt],"complex") == 0 ){  /* not listed help */
            MRG_datum = MRI_complex ;
         } else {
            fprintf(stderr,"-datum of type '%s' is not supported in 3dmerge!\n",
                    argv[nopt] ) ;
            exit(1) ;
         }
         nopt++ ; continue ;  /* go to next arg */
      }

      /**** -thdatum type ****/

      if( strncmp(argv[nopt],"-thdatum",6) == 0 ){
DUMP2 ;
         if( ++nopt >= argc ){
            fprintf(stderr,"need an argument after -thdatum!\n") ; exit(1) ;
         }

         if( strcmp(argv[nopt],"short") == 0 ){
            MRG_thdatum = MRI_short ;
         } else if( strcmp(argv[nopt],"float") == 0 ){
            MRG_thdatum = MRI_float ;
         } else if( strcmp(argv[nopt],"byte") == 0 ){
            MRG_thdatum = MRI_byte ;
         } else {
            fprintf(stderr,"-thdatum of type '%s' is not supported in 3dmerge!\n",
                    argv[nopt] ) ;
            exit(1) ;
         }
         MRG_keepthr = 1 ;    /* -thdatum is meaningless without this */
         nopt++ ; continue ;  /* go to next arg */
      }

      /**** -ghits count ****/

      if( strncmp(argv[nopt],"-ghits",6) == 0 ){
DUMP2 ;
         nopt++ ;
         if( nopt >= argc ){
            fprintf(stderr,"need argument after -ghits!\n") ; exit(1) ;
         }
         MRG_hits_g = strtod( argv[nopt++] , NULL ) ;
         if( MRG_hits_g <= 0 ){
            fprintf(stderr,"illegal value after -ghits\n") ;
            exit(1) ;
         }
         continue ;
      }

      /**** -gclust rmm vmul ****/

      if( strncmp(argv[nopt],"-gclust",6) == 0 ){
DUMP3 ;
         nopt++ ;
         if( nopt+1 >= argc ){
            fprintf(stderr,"need 2 arguments after -gclust!\n") ;
            exit(1) ;
         }
         MRG_clust_rmm_g  = strtod( argv[nopt++] , NULL ) ;
         MRG_clust_vmul_g = strtod( argv[nopt++] , NULL ) ;
         if( MRG_clust_rmm_g <= 0 || MRG_clust_vmul_g <= 0 ){
            fprintf(stderr,"illegal value after -gclust\n") ;
            exit(1) ;
         }
         continue ;
      }

      /**** -session dirname ****/

      if( strncmp(argv[nopt],"-session",6) == 0 ){
DUMP2 ;
         nopt++ ;
         if( nopt >= argc ){
            fprintf(stderr,"need argument after -session!\n") ;
            exit(1) ;
         }
         MCW_strncpy( MRG_output_session , argv[nopt++] , THD_MAX_NAME ) ;
         continue ;
      }

      /**** -prefix prefix ****/

      if( strncmp(argv[nopt],"-prefix",6) == 0 ){
DUMP2 ;
         nopt++ ;
         if( nopt >= argc ){
            fprintf(stderr,"need argument after -prefix!\n") ;
            exit(1) ;
         }
         MCW_strncpy( MRG_output_prefix , argv[nopt++] , THD_MAX_PREFIX ) ;
         continue ;
      }

#if 0
      /**** -label string ****/

      if( strncmp(argv[nopt],"-label",6) == 0 ){
DUMP2 ;
         nopt++ ;
         if( nopt >= argc ){
            fprintf(stderr,"need argument after -label!\n") ;
            exit(1) ;
         }
         MCW_strncpy( MRG_output_label , argv[nopt++] , THD_MAX_LABEL ) ;
         continue ;
      }
#endif

      /**** -nscale [15 Sep 2000] ****/

      if( strcmp(argv[nopt],"-nscale") == 0 ){
DUMP1 ;
         MRG_nscale = 1 ;
         nopt++ ; continue ;
      }


      /**** -gmean ****/

      if( strncmp(argv[nopt],"-gmean",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_MEAN ;
         nopt++ ; continue ;
      }

      /**** -gfisher ****/

      if( strncmp(argv[nopt],"-gfisher",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_FISHER ;
         nopt++ ; continue ;
      }

      /**** -tgfisher (29 Aug 1996) ****/

      if( strncmp(argv[nopt],"-tgfisher",6) == 0 ){
DUMP1 ;
         MRG_cflag_gthr = THFLAG_FICO ;
         nopt++ ; continue ;
      }

      /**** -gnzmean ****/

      if( strncmp(argv[nopt],"-gnzmean",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_NZMEAN ;
         nopt++ ; continue ;
      }

      /**** -gmax ****/

      if( strncmp(argv[nopt],"-gmax",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_MMAX ;
         nopt++ ; continue ;
      }

      /**** -gamax ****/

      if( strncmp(argv[nopt],"-gamax",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_AMAX ;
         nopt++ ; continue ;
      }

      /**** -gsmax ****/

      if( strncmp(argv[nopt],"-gsmax",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_SMAX ;
         nopt++ ; continue ;
      }

      /*** -gcount ****/

      if( strncmp(argv[nopt],"-gcount",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_COUNT ;
         nopt++ ; continue ;
      }

      /*** -gorder ****/

      if( strncmp(argv[nopt],"-gorder",6) == 0 ){
DUMP1 ;
         MRG_cflag_g = CFLAG_ORDER ;
         nopt++ ; continue ;
      }

      /**** unknown switch ****/

      ERROR_message("unrecognized option %s",argv[nopt]) ;
      suggest_best_prog_option(argv[0], argv[nopt]);
      exit(1);

   }  /* end of loop over options */

   /*** cleanup ***/

#if 0
   if( strlen(MRG_output_label) == 0 ){
      MCW_strncpy( MRG_output_label , MRG_output_prefix , THD_MAX_LABEL ) ;
   }
#endif

   /* if -1rank is used, pass the prefix */
   if (MRG_rank) {
      char prefix[THD_MAX_PREFIX];
      FILENAME_TO_PREFIX(MRG_output_prefix, prefix);
      if (prefix[0] == '\0') strcpy(prefix, MRG_output_prefix);
      /*fprintf(stderr,"Prefix >>>%s<<< >>>%s<<<\n",
                       MRG_output_prefix, prefix); */
      if (PREFIX_IS_NIFTI(prefix)) {
         char *p2=NULL;
         if (STRING_HAS_SUFFIX( prefix,".nii")) {
            p2 = strstr( prefix, ".nii");
            prefix[p2- prefix]='\0';
         } else if (STRING_HAS_SUFFIX( prefix,".nii.gz")) {
            p2 = strstr( prefix, ".nii.gz");
            prefix[p2- prefix]='\0';
         }
      }
      sprintf(MRG_edopt.rankmapname,"%s/%s.rankmap.1D",
                                 MRG_output_session, prefix);
   }

   return( nopt );
}

/*------------------------------------------------------------------*/

void MRG_Syntax(void)
{
   printf(
   "This program has 2 different functions:\n"
   " (1) To edit 3D datasets in various ways (threshold, blur, cluster, ...);\n"
   " (2) To merge multiple datasets in various ways (average, max, ...).\n"
   "Either or both of these can be applied.\n"
   "\n"
   "The 'editing' operations are controlled by options that start with '-1',\n"
   "which indicates that they apply to individual datasets\n"
   "(e.g., '-1blur_fwhm').\n"
   "\n"
   "The 'merging' operations are controlled by options that start with '-g',\n"
   "which indicate that they apply to the entire group of input datasets\n"
   "(e.g., '-gmax').\n"
   "\n"
   "----------------------------------------------------------------------\n"
   "Usage: 3dmerge [options] datasets ...\n\n"
   ) ;

   printf(
    "Examples:\n"
    "\n"
    "  1. Apply a 4.0mm FWHM Gaussian blur to EPI run 7.\n"
    "\n"
    "       3dmerge -1blur_fwhm 4.0 -doall -prefix e1.run7_blur run7+orig\n"
    "\n"
    "* These examples are based on a data grid of 3.75 x 3.75 x 3.5, in mm.\n"
    "  So a single voxel has a volume of ~49.22 mm^3 (mvul), and a 40 voxel\n"
    "  cluster has a volume of ~1969 mm^3 (as used in some examples).\n"
    "\n"
    "  2. F-stat only:\n"
    "\n"
    "     Cluster based on a threshold of F=10 (F-stats are in sub-brick #0),\n"
    "     and require a volume of 40 voxels (1969 mm^3).  The output will be\n"
    "     the same F-stats as in the input, but subject to the threshold and\n"
    "     clustering.\n"
    "\n"
    "       3dmerge -1clust 3.76 1969 -1thresh 10.0    \\\n"
    "               -prefix e2.f10 stats+orig'[0]'\n"
    "\n"
    "  3. F-stat only:\n"
    "\n"
    "     Perform the same clustering (as in #2), but apply the radius and\n"
    "     cluster size in terms of cubic millimeter voxels (as if the voxels\n"
    "     were 1x1x1).  So add '-dxyz=1', and adjust rmm and mvul.\n"
    "\n"
    "       3dmerge -dxyz=1 -1clust 1 40 -1thresh 10.0    \\\n"
    "               -prefix e3.f10 stats+orig'[0]'\n"
    "\n"
    "  4. t-stat and beta weight:\n"
    "\n"
    "     For some condition, our beta weight is in sub-brick #4, with the\n"
    "     corresponding t-stat in sub-brick #5.  Cluster based on 40 voxels\n"
    "     and a t-stat threshold of 3.25.  Output the data from the beta\n"
    "     weights, not the t-stats.\n"
    "\n"
    "       3dmerge -dxyz=1 -1clust 1 40 -1thresh 3.25    \\\n"
    "               -1tindex 5 -1dindex 4                 \\\n"
    "               -prefix e4.t3.25 stats+orig\n"
    "\n"
    "  5. t-stat mask:\n"
    "\n"
    "     Apply the same threshold and cluster as in #4, but output a mask.\n"
    "     Since there are 5 clusters found in this example, the values in\n"
    "     the mask will be from 1 to 5, representing the largest cluster to\n"
    "     the smallest.  Use -1clust_order on sub-brick 5.\n"
    "\n"
    "       3dmerge -dxyz=1 -1clust_order 1 40 -1thresh 3.25    \\\n"
    "               -prefix e5.mask5 stats+orig'[5]'\n"
    "\n"
    "     Note: this should match the 3dclust output from:\n"
    "\n"
    "       3dclust -1thresh 3.25 -dxyz=1 1 40 stats+orig'[5]'\n"
    "\n"
    "----------------------------------------------------------------------\n"
   ) ;

   printf( "%s\n" , EDIT_options_help() ) ;

   printf(
    "OTHER OPTIONS:\n"
    "  -nozero     = Do NOT write the output dataset if it would be all zero.\n"
    "\n"
    "  -datum type = Coerce the output data to be stored as the given type,\n"
    "                  which may be byte, short, or float.\n"
    "          N.B.: Byte data cannot be negative.  If this datum type is chosen,\n"
    "                  any negative values in the edited and/or merged dataset\n"
    "                  will be set to zero.\n"
    "\n"
    "  -keepthr    = When using 3dmerge to edit exactly one dataset of a\n"
    "                  functional type with a threshold statistic attached,\n"
    "                  normally the resulting dataset is of the 'fim'\n"
    "                  (intensity only) type.  This option tells 3dmerge to\n"
    "                  copy the threshold data (unedited in any way) into\n"
    "                  the output dataset.\n"
    "          N.B.: This option is ignored if 3dmerge is being used to\n"
    "                  combine 2 or more datasets.\n"
    "          N.B.: The -datum option has no effect on the storage of the\n"
    "                  threshold data.  Instead use '-thdatum type'.\n"
    "\n"
    "  -doall      = Apply editing and merging options to ALL sub-bricks \n"
    "                  uniformly in a dataset.\n"
    "          N.B.: All input datasets must have the same number of sub-bricks\n"
    "                  when using the -doall option. \n"
    "          N.B.: The threshold specific options (such as -1thresh, \n"
    "                  -keepthr, -tgfisher, etc.) are not compatible with \n"
    "                  the -doall command.  Neither are the -1dindex or\n"
    "                  the -1tindex options.\n"
    "          N.B.: All labels and statistical parameters for individual \n"
    "                  sub-bricks are copied from the first dataset.  It is \n"
    "                  the responsibility of the user to verify that these \n"
    "                  are appropriate.  Note that sub-brick auxiliary data \n"
    "                  can be modified using program 3drefit. \n"
    "\n"
    "  -quiet      = Reduce the number of messages shown\n"
    "\n"
    "  -1dindex j  = Uses sub-brick #j as the data source , and uses sub-brick\n"
    "  -1tindex k  = #k as the threshold source.  With these, you can operate\n"
    "                  on any given sub-brick of the inputs dataset(s) to produce\n"
    "                  as output a 1 brick dataset.  If desired, a collection\n"
    "                  of 1 brick datasets can later be assembled into a\n"
    "                  multi-brick bucket dataset using program '3dbucket'\n"
    "                  or into a 3D+time dataset using program '3dTcat'.\n"
    "          N.B.: If these options aren't used, j=0 and k=1 are the defaults\n"
    "\n"
    "  The following option allows you to specify a mask dataset that\n"
    "  limits the action of the 'filter' options to voxels that are\n"
    "  nonzero in the mask:\n"
    "\n"
    "  -1fmask mset = Read dataset 'mset' (which can include a\n"
    "                  sub-brick specifier) and use the nonzero\n"
    "                  voxels as a mask for the filter options.\n"
    "                  Filtering calculations will not use voxels\n"
    "                  that are outside the mask.  If an output\n"
    "                  voxel does not have ANY masked voxels inside\n"
    "                  the rmm radius, then that output voxel will\n"
    "                  be set to 0.\n"
    "         N.B.: * Only the -1filter_* and -t1filter_* options are\n"
    "                 affected by -1fmask.\n"
    "               * Voxels NOT in the fmask will be set to zero in the\n"
    "                 output when the filtering occurs.  THIS IS NEW BEHAVIOR,\n"
    "                 as of 11 Oct 2007.  Previously, voxels not in the fmask,\n"
    "                 but within 'rmm' of a voxel in the mask, would get a\n"
    "                 nonzero output value, as those nearby voxels would be\n"
    "                 combined (via whatever '-1f...' option was given).\n"
    "               * If you wish to restore this old behavior, where non-fmask\n"
    "                 voxels can get nonzero output, then use the new option\n"
    "                 '-1fm_noclip' in addition to '-1fmask'. The two comments\n"
    "                 below apply to the case where '-1fm_noclip' is given!\n"
    "                 * In the linear averaging filters (_mean, _nzmean,\n"
    "                   and _expr), voxels not in the mask will not be used\n"
    "                   or counted in either the numerator or denominator.\n"
    "                   This can give unexpected results if you use '-1fm_noclip'.\n"
    "                   For example, if the mask is designed to exclude the volume\n"
    "                   outside the brain, then voxels exterior to the brain,\n"
    "                   but within 'rmm', will have a few voxels inside the brain\n"
    "                   included in the filtering.  Since the sum of weights (the\n"
    "                   denominator) is only over those few intra-brain\n"
    "                   voxels, the effect will be to extend the significant\n"
    "                   part of the result outward by rmm from the surface\n"
    "                   of the brain.  In contrast, without the mask, the\n"
    "                   many small-valued voxels outside the brain would\n"
    "                   be included in the numerator and denominator sums,\n"
    "                   which would barely change the numerator (since the\n"
    "                   voxel values are small outside the brain), but would\n"
    "                   increase the denominator greatly (by including many\n"
    "                   more weights).  The effect in this case (no -1fmask)\n"
    "                   is to make the filtering taper off gradually in the\n"
    "                   rmm-thickness shell around the brain.\n"
    "                 * Thus, if the -1fmask is intended to clip off non-brain\n"
    "                   data from the filtering, its use should be followed by\n"
    "                   masking operation using 3dcalc:\n"
    "   3dmerge -1filter_aver 12 -1fm_noclip -1fmask mask+orig -prefix x input+orig\n"
    "   3dcalc  -a x -b mask+orig -prefix y -expr 'a*step(b)'\n"
    "   rm -f x+orig.*\n"
    "                 The desired result is y+orig - filtered using only\n"
    "                 brain voxels (as defined by mask+orig), and with\n"
    "                 the output confined to the brain voxels as well.\n"
    "\n"
    "  The following option allows you to specify an almost arbitrary\n"
    "  weighting function for 3D linear filtering:\n"
    "\n"
    "  -1filter_expr rmm expr\n"
    "     Defines a linear filter about each voxel of radius 'rmm' mm.\n"
    "     The filter weights are proportional to the expression evaluated\n"
    "     at each voxel offset in the rmm neighborhood.  You can use only\n"
    "     these symbols in the expression:\n"
    "         r = radius from center\n"
    "         x = dataset x-axis offset from center\n"
    "         y = dataset y-axis offset from center\n"
    "         z = dataset z-axis offset from center\n"
    "         i = x-axis index offset from center\n"
    "         j = y-axis index offset from center\n"
    "         k = z-axis index offset from center\n"
    "     Example:\n"
    "       -1filter_expr 12.0 'exp(-r*r/36.067)'\n"
    "     This does a Gaussian filter over a radius of 12 mm.  In this\n"
    "     example, the FWHM of the filter is 10 mm. [in general, the\n"
    "     denominator in the exponent would be 0.36067 * FWHM * FWHM.\n"
    "     This is one way to get a Gaussian blur combined with the\n"
    "     -1fmask option.  The radius rmm=12 is chosen where the weights\n"
    "     get smallish.]  Another example:\n"
    "       -1filter_expr 20.0 'exp(-(x*x+16*y*y+z*z)/36.067)'\n"
    "     which is a non-spherical Gaussian filter.\n"
    "\n"
    "  ** For shorthand, you can also use the new option (11 Oct 2007)\n"
    "  -1filter_blur fwhm\n"
    "        which is equivalent to\n"
    "  -1filter_expr 1.3*fwhm 'exp(-r*r/(.36067*fwhm*fwhm)'\n"
    "        and will implement a Gaussian blur.  The only reason to do\n"
    "        Gaussian blurring this way is if you also want to use -1fmask!\n"
    "\n"
    "  The following option lets you apply a 'Winsor' filter to the data:\n"
    "\n"
    "  -1filter_winsor rmm nw\n"
    "     The data values within the radius rmm of each voxel are sorted.\n"
    "     Suppose there are 'N' voxels in this group.  We index the\n"
    "     sorted voxels as s[0] <= s[1] <= ... <= s[N-1], and we call the\n"
    "     value of the central voxel 'v' (which is also in array s[]).\n"
    "                 If v < s[nw]    , then v is replaced by s[nw]\n"
    "       otherwise If v > s[N-1-nw], then v is replace by s[N-1-nw]\n"
    "       otherwise v is unchanged\n"
    "     The effect is to increase 'too small' values up to some\n"
    "     middling range, and to decrease 'too large' values.\n"
    "     If N is odd, and nw=(N-1)/2, this would be a median filter.\n"
    "     In practice, I recommend that nw be about N/4; for example,\n"
    "       -dxyz=1 -1filter_winsor 2.5 19\n"
    "     is a filter with N=81 that gives nice results.\n"
    "   N.B.: This option is NOT affected by -1fmask\n"
    "   N.B.: This option is slow! and experimental.\n"
    "\n"
    "  The following option returns a rank value at each voxel in \n"
    "  the input dataset.\n"
    "  -1rank \n"
    "     If the input voxels were, say, 12  45  9  0  9  12  0\n"
    "     the output would be             2   3  1  0  1   2  0\n"
    "     This option is handy for turning FreeSurfer's segmentation\n"
    "     volumes to ROI volumes that can be easily colorized with AFNI.\n"
    "     For example:\n"
    "     3dmerge -1rank -prefix aparc+aseg_rank aparc+aseg.nii \n"
    "     To view aparc+aseg_rank+orig, use the ROI_128 colormap\n"
    "     and set the colorbar range to 128.\n"
    "     The -1rank option also outputs a 1D file that contains \n"
    "     the mapping from the input dataset to the ranked output.\n"
    "     Sub-brick float factors are ignored.\n"
    "\n"
    "     This option only works on datasets of integral values or \n"
    "     of integral data types. 'float' values are typecast to 'int' \n"
    "     before being ranked.\n"
    "\n"
    "     See also program 3dRank\n"
    "\n"
    "MERGING OPTIONS APPLIED TO FORM THE OUTPUT DATASET:\n"
    " [That is, different ways to combine results. The]\n"
    " [following '-g' options are mutually exclusive! ]\n"
    "  -gmean     = Combine datasets by averaging intensities\n"
    "                 (including zeros) -- this is the default\n"
    "  -gnzmean   = Combine datasets by averaging intensities\n"
    "                 (not counting zeros)\n"
    "  -gmax      = Combine datasets by taking max intensity\n"
    "                 (e.g., -7 and 2 combine to 2)\n"
    "  -gamax     = Combine datasets by taking max absolute intensity\n"
    "                 (e.g., -7 and 2 combine to 7)\n"
    "  -gsmax     = Combine datasets by taking max signed intensity\n"
    "                 (e.g., -7 and 2 combine to -7)\n"
    "  -gcount    = Combine datasets by counting number of 'hits' in\n"
    "                  each voxel (see below for definition of 'hit')\n"
    "  -gorder    = Combine datasets in order of input:\n"
    "                * If a voxel is nonzero in dataset #1, then\n"
    "                    that value goes into the voxel.\n"
    "                * If a voxel is zero in dataset #1 but nonzero\n"
    "                    in dataset #2, then the value from #2 is used.\n"
    "                * And so forth: the first dataset with a nonzero\n"
    "                    entry in a given voxel 'wins'\n"
    "  -gfisher   = Takes the arctanh of each input, averages these,\n"
    "                  and outputs the tanh of the average.  If the input\n"
    "                  datum is 'short', then input values are scaled by\n"
    "                  0.0001 and output values by 10000.  This option\n"
    "                  is for merging bricks of correlation coefficients.\n"
    "\n"
    "  -nscale    = If the output datum is shorts, don't do the scaling\n"
    "                  to the max range [similar to 3dcalc's -nscale option]\n"
    "\n"
    "MERGING OPERATIONS APPLIED TO THE THRESHOLD DATA:\n"
    " [That is, different ways to combine the thresholds.  If none of these ]\n"
    " [are given, the thresholds will not be merged and the output dataset  ]\n"
    " [will not have threshold data attached.  Note that the following '-tg']\n"
    " [command line options are mutually exclusive, but are independent of  ]\n"
    " [the '-g' options given above for merging the intensity data values.  ]\n"
    "  -tgfisher  = This option is only applicable if each input dataset\n"
    "                  is of the 'fico' or 'fith' types -- functional\n"
    "                  intensity plus correlation or plus threshold.\n"
    "                  (In the latter case, the threshold values are\n"
    "                  interpreted as correlation coefficients.)\n"
    "                  The correlation coefficients are averaged as\n"
    "                  described by -gfisher above, and the output\n"
    "                  dataset will be of the fico type if all inputs\n"
    "                  are fico type; otherwise, the output datasets\n"
    "                  will be of the fith type.\n"
    "         N.B.: The difference between the -tgfisher and -gfisher\n"
    "                  methods is that -tgfisher applies to the threshold\n"
    "                  data stored with a dataset, while -gfisher\n"
    "                  applies to the intensity data.  Thus, -gfisher\n"
    "                  would normally be applied to a dataset created\n"
    "                  from correlation coefficients directly, or from\n"
    "                  the application of the -1thtoin option to a fico\n"
    "                  or fith dataset.\n"
    "\n"
    "OPTIONAL WAYS TO POSTPROCESS THE COMBINED RESULTS:\n"
    " [May be combined with the above methods.]\n"
    " [Any combination of these options may be used.]\n"
    "  -ghits count     = Delete voxels that aren't !=0 in at least\n"
    "                       count datasets (!=0 is a 'hit')\n"
    "  -gclust rmm vmul = Form clusters with connection distance rmm\n"
    "                       and clip off data not in clusters of\n"
    "                       volume at least vmul microliters\n"
    "\n"
    "The '-g' and '-tg' options apply to the entire group of input datasets.\n"
    "\n"

    "OPTIONS THAT CONTROL THE NAMES OF THE OUTPUT DATASET:\n"
    "  -session dirname  = write output into given directory (default=./)\n"
    "  -prefix  pname    = use 'pname' for the output dataset prefix\n"
    "                       (default=mrg)\n"
#if 0
    "  -label   string   = use 'string' for the label in the output\n"
    "                       dataset (the label is used for switching\n"
    "                       between datasets in AFNI)\n"
#endif
    "\n"

    "NOTES:\n"
    " **  If only one dataset is read into this program, then the '-g'\n"
    "       options do not apply, and the output dataset is simply the\n"
    "       '-1' options applied to the input dataset (i.e., edited).\n"
    " **  A merged output dataset is ALWAYS of the intensity-only variety.\n"
    " **  You can combine the outputs of 3dmerge with other sub-bricks\n"
    "       using the program 3dbucket.\n"
    " **  Complex-valued datasets cannot be merged.\n"
    " **  This program cannot handle time-dependent datasets without -doall.\n"
    " **  Note that the input datasets are specified by their .HEAD files,\n"
    "       but that their .BRIK files must exist also!\n"

#ifdef ALLOW_SUBV
    "\n"
    MASTER_SHORTHELP_STRING
    "\n"
    " ** Input datasets using sub-brick selectors are treated as follows:\n"
    "      - 3D+time if the dataset is 3D+time and more than 1 brick is chosen\n"
    "      - otherwise, as bucket datasets (-abuc or -fbuc)\n"
    "       (in particular, fico, fitt, etc. datasets are converted to fbuc)\n"
    " ** If you are NOT using -doall, and choose more than one sub-brick\n"
    "     with the selector, then you may need to use -1dindex to further\n"
    "     pick out the sub-brick on which to operate (why you would do this\n"
    "     I cannot fathom).  If you are also using a thresholding operation\n"
    "     (e.g., -1thresh), then you also MUST use -1tindex to choose which\n"
    "     sub-brick counts as the 'threshold' value.  When used with sub-brick\n"
    "     selection, 'index' refers the dataset AFTER it has been read in:\n"
    "          -1dindex 1 -1tindex 3 'dset+orig[4..7]'\n"
    "     means to use the #5 sub-brick of dset+orig as the data for merging\n"
    "     and the #7 sub-brick of dset+orig as the threshold values.\n"
    " ** The above example would better be done with\n"
    "          -1tindex 1 'dset+orig[5,7]'\n"
    "     since the default data index is 0. (You would only use -1tindex if\n"
    "     you are actually using a thresholding operation.)\n"
    " ** -1dindex and -1tindex apply to all input datasets.\n"
#endif
   ) ;
   PRINT_COMPILE_DATE ; exit(0) ;
}

/*------------------------------------------------------------------*/

int main( int argc , char * argv[] )
{
   int file_num , first_file , nx,ny,nz , nxyz , ii , num_dset ,
       file_count , ptmin , iclu,nclu , edit_type , ival,ivout , tval ;
   float dx,dy,dz , fac , dxyz , rmm,vmul ;
   THD_3dim_dataset * dset=NULL , * new_dset=NULL ;
   THD_3dim_dataset ** dsetar=NULL ;                 /* Nov 1998 */
   short * gnum=NULL ;
   float * gfim=NULL , * tfim=NULL , * ggfim=NULL , * ttfim=NULL ;
   int     datum ;
   MCW_cluster_array * clar ;
   float fimfac , fimfacinv , first_fimfac , thrfac ;
   int   output_datum , output_thdatum ;
   int   input_datum  , input_thdatum , first_datum ;

   float thr_stataux[MAX_STAT_AUX] ;
   int   num_fico=0;
   int   is_int=1 ;   /* 08 Jan 1998 */

   int iv, iv_bot, iv_top;      /* dataset sub-brick indices    02 Feb 1998 */

   /*----- identify program -----*/
   printf ("Program %s \n", PROGRAM_NAME);

   /*** read input options ***/

   if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) MRG_Syntax() ;

   /*-- 20 Apr 2001: addto the arglist, if user wants to [RWCox] --*/

   mainENTRY("3dmerge main") ; machdep() ; PRINT_VERSION("3dmerge") ;

   { int new_argc ; char ** new_argv ;
     addto_args( argc , argv , &new_argc , &new_argv ) ;
     if( new_argv != NULL ){ argc = new_argc ; argv = new_argv ; }
   }

   AFNI_logger("3dmerge",argc,argv) ;

   first_file = MRG_read_opts( argc , argv ) ;
   file_count = argc - first_file ;            /* number of datasets input */

   if( ! MRG_be_quiet )
      printf("3dmerge: edit and combine 3D datasets, by RW Cox\n") ;

   if( first_file < 1 || first_file >= argc ){
      fprintf(stderr,"*** ILLEGAL COMMAND LINE ***\n") ; exit(1) ;
   }

   /*----- check for compatibility of user options -----*/ /* 02 Feb 1998 */
   if (MRG_doall)
     { int nerr = 0 ;

       if( MRG_ivfim >= 0 ){   /* Nov 1998 */
         fprintf(stderr,"-1dindex is not compatible with -doall option \n");
         nerr++ ;  }
       if( MRG_ivthr >= 0 ){   /* Nov 1998 */
         fprintf(stderr,"-1dindex is not compatible with -doall option \n");
         nerr++ ;  }

       if (MRG_edopt.thtoin > 0)  {
	 fprintf (stderr, "-1thtoin is not compatible with -doall option \n");
	 nerr++ ;  }
       if (MRG_edopt.thresh > 0.0)  {
	 fprintf (stderr, "-1thresh is not compatible with -doall option \n");
	 nerr++ ;  }
       if (MRG_edopt.thrfilter_opt > 0)  {
	 fprintf (stderr, "-t1filter is not compatible with -doall option \n");
	 nerr++ ;  }
       if (MRG_edopt.thrblur > 0.0)  {
	 fprintf (stderr, "-t1blur is not compatible with -doall option \n");
	 nerr++ ;  }
       if (MRG_thdatum >= 0)  {
	 fprintf (stderr, "-thdatum is not compatible with -doall option \n");
	 nerr++ ;  }
       if (MRG_keepthr)  {
	 fprintf (stderr, "-keepthr is not compatible with -doall option \n");
	 nerr++ ;  }
       if (MRG_cflag_gthr > 0)  {
	 fprintf (stderr, "-tgfisher is not compatible with -doall option \n");
	 nerr++ ;  }

       if( nerr > 0 ) exit(1) ;
     }

   /* Nov 1998: other checks */

   if( (MRG_ivfim >= 0 || MRG_ivthr >=0) && MRG_keepthr ){
      fprintf(stderr,"-keepthr is not compatible with -1dindex or -1tindex\n") ;
      exit(1) ;
   }

   if( (MRG_ivfim >= 0 || MRG_ivthr >=0) && MRG_cflag_gthr ){
      fprintf(stderr,"-tgfisher is not compatible with -1dindex or -1tindex\n") ;
      exit(1) ;
   }

   /* check for existence of each input data set. */   /* 09 December 1996 */
   for (file_num = first_file;  file_num < argc;  file_num++)
     {
       dset = DSET_OPEN( argv[file_num] ) ;
       if( ! ISVALID_3DIM_DATASET(dset) )
	 {
	   fprintf(stderr,"*** cannot open dataset %s\n",argv[file_num]) ;
	   exit(1) ;
         }

       /* Nov 1998: check user-controlled brick indexes */

       if( MRG_ivfim >= DSET_NVALS(dset) ){
          fprintf(stderr,
                  "*** Dataset %s does not have enough bricks for -1dindex %d\n" ,
                  argv[file_num],MRG_ivfim) ;
          exit(1) ;
       }
       if( MRG_ivthr >= DSET_NVALS(dset) ){
          fprintf(stderr,
                  "*** Dataset %s does not have enough bricks for -1tindex %d\n" ,
                  argv[file_num],MRG_ivthr) ;
          exit(1) ;
       }

       THD_delete_3dim_dataset( dset , False ) ; dset = NULL ;
     }

   /* 02 Mar 2001: print a message about dataset indices */

   if( MRG_ivfim < 0 ) fprintf(stderr,"++ default -1dindex = 0\n") ;
   if( MRG_ivthr < 0 ) fprintf(stderr,"++ default -1tindex = 1\n") ;

   /* read first dataset */

   dset = DSET_OPEN( argv[first_file] ) ;
   if( ! ISVALID_3DIM_DATASET(dset) ){
      fprintf(stderr,"*** Unable to open first dataset %s\n",argv[first_file]) ;
      exit(1) ;
   }

   if( DSET_NUM_TIMES(dset) > 1 && (!MRG_doall && MRG_ivfim < 0) ){
      fprintf(stderr, "*** Unable to merge time-dependent datasets"
                      " without -doall or -1dindex\n") ;
      exit(1) ;
   }

   /* get the dimensions */

   nx = dset->daxes->nxx ;
   ny = dset->daxes->nyy ;
   nz = dset->daxes->nzz ; nxyz = nx*ny*nz ;

   dx = fabs(dset->daxes->xxdel) ;
   dy = fabs(dset->daxes->yydel) ;
   dz = fabs(dset->daxes->zzdel) ;

   if( MRG_edopt.fake_dxyz ) dx = dy = dz = 1.0 ;  /* 11 Sep 2000 */

#if 0
   nice(1) ;  /* slow us down, a little */
#endif

   if( MRG_edopt.nfmask > 0 && MRG_edopt.nfmask != nxyz ){
      fprintf(stderr,
              "*** -1fmask and 1st dataset don't have same number of voxels\n\a");
      exit(1) ;
   }

   /*******************************************************************/
   /****      if only one file, edit it, modify its names ...      ****/
   /****      then write the modified dataset to disk              ****/

   if( file_count == 1 ){

      ival        = DSET_PRINCIPAL_VALUE(dset) ;
      input_datum = DSET_BRICK_TYPE(dset,ival) ;
      if( MRG_datum >= 0 ) output_datum = MRG_datum ;
      else                 output_datum = input_datum ;

#if 1
      /** 17 Sep 1998: Move the creation of the new dataset
                       to AFTER the editing operations, so that any
                       dataset parameter changes in EDIT_one_dataset
                       will properly be propagated to the new dataset.
                       The creation shown here is just to check if the
                       output dataset exists already.                   **/

     if( THD_deathcon() ){
      new_dset = EDIT_empty_copy( dset ) ;

      EDIT_dset_items( new_dset ,
                          ADN_prefix         , MRG_output_prefix ,
                          ADN_directory_name , MRG_output_session ,
                       ADN_none ) ;

      if( THD_is_file(new_dset->dblk->diskptr->header_name) ){
         fprintf(stderr,
                 "*** Output file %s already exists -- cannot continue!\n",
                 new_dset->dblk->diskptr->header_name ) ;
         exit(1) ;
      }

      THD_delete_3dim_dataset( new_dset , False ) ;  /* toss this junk */
     }
#endif

      /* 05 Jul 2002: illegal uses of -keepthr? */

      if( MRG_keepthr && !ISFUNC(dset) ){
        fprintf(stderr,"*** -keepthr can't be used on non-functional dataset!\n");
        exit(1) ;
      }
      if( MRG_keepthr && dset->func_type == FUNC_FIM_TYPE ){
        fprintf(stderr,"*** -keepthr can't be used on 'fim' type dataset!\n") ;
        exit(1) ;
      }
      if( MRG_keepthr && dset->func_type == FUNC_BUCK_TYPE ){
        fprintf(stderr,"*** -keepthr can't be used on 'fbuc' type dataset!\n"
                       "    You can use '3dbuc2fim' first, if needed.\n"    );
        exit(1) ;
      }

      /** get ready to go **/

      if( ! MRG_be_quiet ){
         printf("++ editing input dataset in memory (%.1f MB)\n",
                ((double)dset->dblk->total_bytes) / MEGA ) ;
         fflush(stdout) ;
      }

  /* 02 Feb 1998 */
  if (MRG_doall)
    {  iv_bot = 0;  iv_top = DSET_NVALS(dset);  }
  else if( MRG_ivfim >= 0 )                       /* Nov 1998 */
    {  iv_bot = MRG_ivfim ; iv_top = iv_bot+1 ; }
  else
    {  iv_bot = DSET_PRINCIPAL_VALUE(dset);  iv_top = iv_bot + 1;  }

  /*----- Iterate over sub-bricks -----*/
  for (iv = iv_bot;  iv < iv_top;  iv++)
    {
      /* show if requested */
      if ((MRG_verbose) && MRG_doall) printf ("Editing sub-brick %d\n", iv);

      MRG_edopt.iv_fim = iv;
      EDIT_one_dataset( dset , &MRG_edopt ) ;  /* all the real work */

      if( !MRG_be_quiet && !MRG_verbose ) { printf(".") ; fflush(stdout) ; }

    }  /* end of loop over sub-bricks */

    /* newline after the pacifier dots */
    if( !MRG_be_quiet && !MRG_verbose ) { printf("\n") ; fflush(stdout) ; }

    if( MRG_edopt.nfmask > 0 ){
      free(MRG_edopt.fmask) ; MRG_edopt.fmask = NULL ; MRG_edopt.nfmask = 0 ;
    }

    if( !MRG_be_quiet ) printf("\n") ;

      /** 17 Sep 1998: NOW create the new dataset **/

      new_dset = EDIT_empty_copy( dset ) ;

      tross_Copy_History( dset , new_dset ) ;
      tross_Make_History( "3dmerge" , argc , argv , new_dset ) ;

      EDIT_dset_items( new_dset ,
                          ADN_prefix , MRG_output_prefix ,
                          ADN_label1 , MRG_output_prefix ,
                          ADN_directory_name , MRG_output_session ,
                       ADN_none ) ;
      strcat( new_dset->self_name , "(ED)" ) ;
                                                           /* 02 Feb 1998 */
      if( (! MRG_keepthr) && (new_dset->dblk->nvals > 1) && (! MRG_doall) )
         EDIT_dset_items( new_dset ,
                             ADN_nvals , 1 ,
                             ADN_ntt   , 0 ,                 /* Nov 1998 */
                             ADN_func_type , FUNC_FIM_TYPE ,
                          ADN_none ) ;

      if( MRG_keepthr && ISFUNC(new_dset) && FUNC_HAVE_THR(new_dset->func_type) ){
         ii            = FUNC_ival_thr[dset->func_type] ;
         input_thdatum = DSET_BRICK_TYPE(dset,ii) ;
         if( MRG_thdatum >= 0 ) output_thdatum = MRG_thdatum ;
         else                   output_thdatum = input_thdatum ;
      } else {
         output_thdatum = input_thdatum = ILLEGAL_TYPE ;
      }

      /** Coerce the output data type into a new brick, if needed **/

  /*----- Iterate over sub-bricks [again] -----*/
  for (iv = iv_bot;  iv < iv_top;  iv++)
    {
  /* 02 Feb 1998 */
  if (MRG_doall)
    {
      ival  = iv;
      ivout = iv;
    }
  else if( MRG_ivfim >= 0 )
    {
      ival  = MRG_ivfim ;
      ivout = DSET_PRINCIPAL_VALUE(new_dset) ;
    }
  else
    {
      ival  = DSET_PRINCIPAL_VALUE(dset) ;
      ivout = DSET_PRINCIPAL_VALUE(new_dset) ;
    }

      if( input_datum == output_datum ){

         /** Attach the brick of the input dataset to the brick of the output.  **/
         /** (This isn't exactly kosher, but we are exiting almost immediately) **/

         mri_fix_data_pointer( DSET_ARRAY(dset,ival) , DSET_BRICK(new_dset,ivout) ) ;

#if 1
#   if 0
         if( ivout != ival )
#   endif
            DSET_BRICK_FACTOR(new_dset,ivout) = DSET_BRICK_FACTOR(dset,ival) ;
#endif

      } else {

         /** Must create a new brick and do the conversion **/

         void * dfim , * efim ;
         float efac = DSET_BRICK_FACTOR(dset,ival) ;

         if( ! MRG_be_quiet ){
            printf("-- coercing output datum to be %s\n",
                   MRI_TYPE_name[output_datum]);
         }

         efim = DSET_ARRAY(dset,ival) ;
         dfim = (void *) malloc( mri_datum_size(output_datum) * nxyz ) ;
         if( dfim == NULL ){
            fprintf(stderr,"*** Can't malloc output brick #%d\n",ivout); exit(1);
         }

         /** 03 Dec 1998: scale to integer and float types separately **/

         if( MRI_IS_INT_TYPE(output_datum) ){
            fimfac = EDIT_coerce_autoscale( nxyz , input_datum  , efim ,
                                                   output_datum , dfim  ) ;
            if( fimfac == 0.0 ) fimfac  = 1.0 ;
            if( efac   != 0.0 ) fimfac /= efac ;

            DSET_BRICK_FACTOR(new_dset,ivout) = (fimfac != 0.0 && fimfac != 1.0)
                                                ? 1.0/fimfac : 0.0 ;
         } else {

            EDIT_coerce_scale_type( nxyz , efac , input_datum  , efim ,
                                                  output_datum , dfim  ) ;
            DSET_BRICK_FACTOR(new_dset,ivout) = 0.0 ;
         }

         mri_free( DSET_BRICK(dset,ival) ) ;
         EDIT_substitute_brick( new_dset , ivout , output_datum , dfim ) ;
      }

      /** Now do the threshold data [won't happen if doall is also happening] **/

      if( output_thdatum >= 0 ){

         ival  = FUNC_ival_thr[    dset->func_type] ;
         ivout = FUNC_ival_thr[new_dset->func_type] ;

         if( input_thdatum == output_thdatum ){

            mri_fix_data_pointer( DSET_ARRAY(dset,ival),DSET_BRICK(new_dset,ivout) ) ;

#if 0
            DSET_BRICK_FACTOR(new_dset,ivout) = DSET_BRICK_FACTOR(dset,ival) ;
#endif

         } else {
            void * dfim , * efim ;

            if( ! MRG_be_quiet ){
               printf("-- coercing threshold datum to be %s\n",
                      MRI_TYPE_name[output_thdatum]);
            }

            efim = DSET_ARRAY(dset,ival) ;
            dfim = (void *) RwcMalloc( mri_datum_size(output_thdatum) * nxyz ) ;

            switch( output_thdatum ){
               default: fprintf(stderr,"** illegal output_thdatum = %d\n",
                                output_thdatum);
               exit(1) ;

               case MRI_float:
                  fimfacinv = 0.0 ;
                  fimfac    = DSET_BRICK_FACTOR(dset,ival) ;
                  if( fimfac == 0.0 ){
                     fimfac = (input_thdatum == MRI_short)
                               ? 1.0/FUNC_scale_short[dset->func_type]
                               : (input_thdatum == MRI_byte)
                               ? 1.0/FUNC_scale_byte[dset->func_type] : 0.0 ;
                  }
               break ;

               case MRI_short:
                  if( input_datum == MRI_float ){
                     fimfac    = FUNC_scale_short[new_dset->func_type] ;
                     fimfacinv = 1.0 / fimfac ;
                  } else if( input_datum == MRI_byte ){
                     fimfac    = ((float)FUNC_scale_short[new_dset->func_type])
                                / FUNC_scale_byte[new_dset->func_type] ;
                     fimfacinv = 1.0 / FUNC_scale_short[new_dset->func_type] ;
                  } else {
                     fprintf(stderr,"** illegal input_thdatum = %d\n",input_thdatum);
                     exit(1) ;
                  }
               break ;

               case MRI_byte:
                  if( input_datum == MRI_float ){
                     fimfac    = FUNC_scale_byte[new_dset->func_type] ;
                     fimfacinv = 1.0 / fimfac ;
                  } else if( input_datum == MRI_short ){
                     fimfac    = ((float)FUNC_scale_byte[new_dset->func_type])
                                / FUNC_scale_short[new_dset->func_type] ;
                     fimfacinv = 1.0 / FUNC_scale_byte[new_dset->func_type] ;
                  } else {
                     fprintf(stderr,"** illegal input_thdatum = %d\n",input_thdatum);
                     exit(1) ;
                  }
               break ;
            }

            EDIT_coerce_scale_type( nxyz , fimfac ,
                                    DSET_BRICK_TYPE(dset,ival),efim ,
                                    output_thdatum,dfim ) ;

            DSET_BRICK_FACTOR(new_dset,ivout) = fimfacinv ;
            EDIT_substitute_brick( new_dset , ivout , output_thdatum , dfim ) ;
            mri_free( DSET_BRICK(dset,ival) ) ;
         }
      }

    }  /* iv    End of iteration over sub-bricks */


      if( MRG_nozero && THD_count_nonzero_bricks(new_dset) == 0 ){ /* 17 Jan 2017 */
        INFO_message("output dataset is all zero ==> not writing because of -allzero") ;
        exit(0) ;
      }

      THD_load_statistics( new_dset ) ;
      THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
      if( ! MRG_be_quiet )
        fprintf(stderr,"-- Wrote edited dataset: %s\n" , DSET_BRIKNAME(new_dset) ) ;
      exit(0) ;
   }

   /************************************************************************/
   /********         more than one input dataset --> merger         ********/

   /* Nov 1998: make an array of input datasets, load 1st element */

   dsetar = (THD_3dim_dataset **) malloc(sizeof(THD_3dim_dataset *)*file_count);
   for( ii=0 ; ii < file_count ; ii++ ) dsetar[ii] = NULL ;
   dsetar[0] = dset ;

   /* make an empty copy of the first dataset, then modify it */

   new_dset = EDIT_empty_copy( dset ) ;
   ival     = DSET_PRINCIPAL_VALUE(dset) ;

   if( MRG_datum >= 0 ) output_datum = MRG_datum ;
   else                 output_datum = DSET_BRICK_TYPE(dset,ival) ;

   EDIT_dset_items( new_dset ,
                       ADN_prefix , MRG_output_prefix ,
                       ADN_label1 , MRG_output_prefix ,
                       ADN_directory_name , MRG_output_session ,
                    ADN_none ) ;
   strcat( new_dset->self_name , "(MG)" ) ;

   /* 29 Aug 1996: change the dataset type, depending on the merger type */

  if (! MRG_doall)   /* 02 Feb 1998 */
   switch( MRG_cflag_gthr ){
       default:
          EDIT_dset_items( new_dset , ADN_nvals,1 , ADN_ntt,0 , ADN_none ) ;
          if( ISFUNC(dset) )
             EDIT_dset_items( new_dset , ADN_func_type,FUNC_FIM_TYPE , ADN_none ) ;
       break ;

       case THFLAG_FICO:  /* do nothing to the dataset now */
          num_fico = 0 ;
       break ;
   }

   if( THD_deathcon() && THD_is_file(new_dset->dblk->diskptr->header_name) ){
      fprintf(stderr,
              "*** Output file %s already exists -- cannot continue!\n",
              new_dset->dblk->diskptr->header_name ) ;
      exit(1) ;
   }

   if( ! MRG_be_quiet && MRG_keepthr )
      printf("-- ignoring -keepthr option\n") ;


  /* 02 Feb 1998 */
   if (MRG_doall)
     {  iv_bot = 0;  iv_top = DSET_NVALS(dset);  }
   else if( MRG_ivfim >= 0 )                       /* Nov 1998 */
     {  iv_bot = MRG_ivfim ; iv_top = iv_bot+1 ; }
   else
     {  iv_bot = DSET_PRINCIPAL_VALUE(dset);  iv_top = iv_bot + 1;  }

  /*----- Iterate over sub-bricks -----*/

  ivout = 0 ;  /* Nov 1998 */

  for (iv = iv_bot;  iv < iv_top;  iv++)
    {
      if ((!MRG_be_quiet) ) printf ("-- Editing sub-brick %d \n", iv);

      MRG_edopt.iv_fim = iv;

   /* make space for the merger computations */

   tfim = (float *) RwcMalloc( sizeof(float) * nxyz ) ;  /* dataset copy */
   gfim = (float *) RwcMalloc( sizeof(float) * nxyz ) ;  /* results */
   gnum = (short *) RwcMalloc( sizeof(short) * nxyz ) ;  /* counts */
   for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] = 0.0 ;      /* initialize */
   for( ii=0 ; ii < nxyz ; ii++ ) gnum[ii] = 0 ;

   /* 29 Aug 1996: make space for merger of thresholds, if desired */

   if( MRG_cflag_gthr != THFLAG_NONE ){
      ttfim = (float *) RwcMalloc( sizeof(float) * nxyz ) ;  /* thresh copy */
      ggfim = (float *) RwcMalloc( sizeof(float) * nxyz ) ;  /* thresh results */
      for( ii=0 ; ii < nxyz ; ii++ ) ggfim[ii] = 0.0 ;      /* initialize */

      for( ii=0 ; ii < MAX_STAT_AUX ; ii++ ) thr_stataux[ii] = 0 ;
   }

   if( ! MRG_be_quiet ){
      float nbytes = (2.0*sizeof(float)+sizeof(short))*nxyz ;
      if( MRG_cflag_gthr != THFLAG_NONE ) nbytes += 2.0*sizeof(float)*nxyz ;
      printf("-- allocated %.1f MB scratch memory\n", nbytes/MEGA ) ;
   }

   /***--- read datasets, edit them, add them into gfim and gnum ---***/

   num_dset = 0 ;
   for( file_num=first_file; file_num < argc ; file_num++ ){

      /** read dataset header if not already input **/

      dset = dsetar[ file_num - first_file ] ;  /* Nov 1998 */
      if( dset == NULL ){
         dset = dsetar[ file_num - first_file ] = DSET_OPEN( argv[file_num] ) ;
         if( ! ISVALID_3DIM_DATASET(dset) ){
           fprintf(stderr,"*** cannot open dataset %s\n",argv[file_num]); exit(1);
         }
      }

      /* check for dimensional mismatch */

      if( dset->daxes->nxx!=nx || dset->daxes->nyy!=ny || dset->daxes->nzz!=nz )
        ERROR_exit("dataset brick size mismatch at file %s\n",argv[file_num] );
      else if( !EQUIV_GRIDS(dset,new_dset) )
        WARNING_message("dataset grid mismatch at file %s",argv[file_num]);

      /* 02 Feb 1998 */
      if ( (MRG_doall) && (DSET_NVALS(dset) != iv_top) )
         fprintf (stderr, "*** -doall dataset nvals mismatch at file %s\n",
                  argv[file_num]);


      if( DSET_NUM_TIMES(dset) > 1 && (!MRG_doall && MRG_ivfim < 0) ){ /* no time     */
         fprintf(stderr,                                               /* dependence! */
                 "*** cannot use time-dependent dataset %s"
                 " without -doall or -1dindex\n",argv[file_num]) ;
         exit(1) ;
      }

      /* check for dataset type, if needed for the merging operations ordered */

      if( MRG_cflag_gthr == THFLAG_FICO ){

         if( !ISFUNC(dset) ){
            fprintf(stderr,
                  "*** dataset from file %s is anatomical using '-tgfisher'!\n",
                argv[file_num] ) ;
            exit(1) ;
         }

         switch( dset->func_type ){
            default:
               fprintf(stderr,
                "*** dataset from file %s is illegal type using '-tgfisher'!\n",
                   argv[file_num] ) ;
            exit(1) ;

            case FUNC_COR_TYPE:   /* add up degrees-of-freedom */
               num_fico ++ ;
               for( ii=0 ; ii < FUNC_need_stat_aux[FUNC_COR_TYPE] ; ii++ )
                 thr_stataux[ii] += dset->stat_aux[ii] ;
            break ;

            case FUNC_THR_TYPE:  /* do nothing */
            break ;
         }
      }

      /* get the control information about this dataset */

#if 1
      ival = iv ;                  /* Nov 1998 */
#else
      if (MRG_doall)  ival = iv;   /* 02 Feb 1998 */
      else            ival = DSET_PRINCIPAL_VALUE(dset) ;
#endif
      datum = DSET_BRICK_TYPE(dset,ival) ;

      if( ! AFNI_GOOD_FUNC_DTYPE(datum) ){
         fprintf(stderr,"*** Illegal datum for 3dmerge: %s in file %s ***\n" ,
                 MRI_TYPE_name[datum] , argv[file_num] ) ;
         exit(1) ;
      }

      if( ! MRG_be_quiet ){
         printf("-- processing file %s" , argv[file_num] ) ;
         fflush(stdout) ;
      }

      /* mess with the input data, maybe */

      if( MRG_have_edopt )
         EDIT_one_dataset( dset , &MRG_edopt ) ;  /* some real work */
      else
         DSET_load(dset) ;
      CHECK_LOAD_ERROR(dset) ;

      /* 02 Nov 2001: check for data that doesn't exist */

      CHECK_LOAD_ERROR(dset) ;

      if( ! MRG_be_quiet ){ printf(".") ; fflush(stdout) ; }

      /* copy it into tfim , scaling if needed */

      fimfac = DSET_BRICK_FACTOR(dset,ival) ;          /* normal case */

      if( MRG_cflag_g == CFLAG_FISHER             &&   /* special case */
          DSET_BRICK_TYPE(dset,ival) == MRI_short &&
          (fimfac==0.0 || fimfac==1.0)              ){

         fimfac = 1.0 / FUNC_COR_SCALE_SHORT ;
      }

      if( num_dset == 0 ){
         first_fimfac = fimfac ;  /* save for later */
         first_datum  = datum  ;
      }

      /** 08 Jan 1998: check if all inputs are integer types **/

      is_int = is_int && (MRI_IS_INT_TYPE(datum) && fimfac == 0.0) ;

      /* the actual copy+scaling operation */

      EDIT_coerce_scale_type( nxyz , fimfac ,
                              DSET_BRICK_TYPE(dset,ival) , DSET_ARRAY(dset,ival) ,
                              MRI_float , tfim ) ;

      /* 29 Aug 1996: get the threshold data into ttfim , if needed */

      if( MRG_cflag_gthr != THFLAG_NONE && (tval=DSET_THRESH_VALUE(dset)) >= 0 ){

         int thdatum = DSET_BRICK_TYPE(dset,tval) ;

         if( ! AFNI_GOOD_FUNC_DTYPE(thdatum) ){
            fprintf(stderr,"*** Illegal threshold for 3dmerge: %s in file %s ***\n" ,
                    MRI_TYPE_name[thdatum] , argv[file_num] ) ;
            exit(1) ;
         }

         thrfac = DSET_BRICK_FACTOR(dset,tval) ;          /* normal case */

         if( MRG_cflag_gthr == THFLAG_FICO           &&   /* special case */
             DSET_BRICK_TYPE(dset,tval) == MRI_short &&
             (thrfac==0.0 || thrfac==1.0)              ){

            thrfac = 1.0 / FUNC_COR_SCALE_SHORT ;
         }

         EDIT_coerce_scale_type( nxyz , thrfac ,
                                 thdatum , DSET_ARRAY(dset,tval) ,
                                 MRI_float , ttfim ) ;
      }

      DSET_unload(dset) ;  /* Nov 1998: don't need data bricks in memory any more */

      /*** merge tfim into gfim and gnum ***/

      if( MRG_cflag_g == CFLAG_MMAX ){
         for( ii=0 ; ii < nxyz ; ii++ ){
            if( tfim[ii] != 0 ){
               gnum[ii]++ ;
               if( tfim[ii] > gfim[ii] ) gfim[ii] = tfim[ii] ;
            }
         }
      } else if( MRG_cflag_g == CFLAG_AMAX ){
         float dab ;
         for( ii=0 ; ii < nxyz ; ii++ ){
            if( tfim[ii] != 0 ){
               gnum[ii]++ ;
               dab = fabs(tfim[ii]) ;
               if( dab > gfim[ii] ) gfim[ii] = dab ;
            }
         }
      } else if( MRG_cflag_g == CFLAG_SMAX ){
         float dab ;
         for( ii=0 ; ii < nxyz ; ii++ ){
            if( tfim[ii] != 0 ){
               gnum[ii]++ ;
               dab = fabs(tfim[ii]) ;
               if( dab > fabs(gfim[ii]) ) gfim[ii] = tfim[ii] ;
            }
         }
      } else if( MRG_cflag_g == CFLAG_ORDER ){
         for( ii=0 ; ii < nxyz ; ii++ ){
            if( tfim[ii] != 0 ){
               gnum[ii]++ ;
               if( gfim[ii] == 0 ) gfim[ii] = tfim[ii] ;
            }
         }
      } else if( MRG_cflag_g == CFLAG_FISHER ){
         for( ii=0 ; ii < nxyz ; ii++ ){
            if( tfim[ii] != 0 ){
               gnum[ii]++ ; gfim[ii] += ATANH(tfim[ii]) ;
            }
         }
      } else {                                /* default = sum up */
         for( ii=0 ; ii < nxyz ; ii++ ){
            if( tfim[ii] != 0 ){
               gnum[ii]++ ; gfim[ii] += tfim[ii] ;
            }
         }
      }

      /* 29 Aug 1996: merge the threshold data, if any */

      if( MRG_cflag_gthr == THFLAG_FICO ){
         for( ii=0 ; ii < nxyz ; ii++ ) ggfim[ii] += ATANH(ttfim[ii]) ;
      }

      if( ! MRG_be_quiet ){ printf(".\n") ; fflush(stdout) ; }

      num_dset++ ;
   }  /* end of combiner loop over datasets */

   myRwcFree(tfim) ;                      /* not needed any more */
   if( ttfim != NULL ) myRwcFree(ttfim) ;

   if( MRG_edopt.nfmask > 0 ){
     free(MRG_edopt.fmask) ; MRG_edopt.fmask = NULL ; MRG_edopt.nfmask = 0 ;
   }

   /*** if only one dset encountered, some error! ***/

   if( num_dset <= 1 ){
      fprintf(stderr,"*** Only found 1 dataset -- computations aborted!\n") ;
      exit(1) ;
   }

   if( ! MRG_be_quiet ) printf("-- merging results into sub-brick %d\n",ivout) ;

   /*** now, edit the merged dataset:
        cast out voxels that weren't hit enough ***/

   if( MRG_hits_g > 0 ){
      for( ii=0 ; ii < nxyz ; ii++ )
         if( gnum[ii] < MRG_hits_g ) { gfim[ii] = 0 ; gnum[ii] = 0 ; }
   }

   /*** decide if the output is to be stored as integers ***/
   /*   [at this point, is_int is true if all]
         inputs were integers and unscaled.  ] 08 Jan 1998 */

   switch( MRG_cflag_g ){
      default:          is_int = 0 ; break ;   /* not allowed */

      case CFLAG_COUNT: is_int = 1 ; break ;   /* it IS an integer */

      case CFLAG_MMAX:                         /* if all inputs */
      case CFLAG_SMAX:                         /* are integers, */
      case CFLAG_AMAX:                         /* then one of these */
      case CFLAG_ORDER: break ;                /* will be integer also */
   }

   /*** do the averaging as ordered ***/

   switch( MRG_cflag_g ){
      default: break ;

      case CFLAG_COUNT:
         first_fimfac = 0.0 ;
         for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] = gnum[ii] ;
      break ;

      case CFLAG_MEAN:
         fac = 1.0 / num_dset ;
         for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] *= fac ;
      break ;

      case CFLAG_NZMEAN:
         for( ii=0 ; ii < nxyz ; ii++ )
            if( gnum[ii] > 0 ) gfim[ii] /= gnum[ii] ;
      break ;

      case CFLAG_FISHER:
         fac = 1.0 / num_dset ;
         for( ii=0 ; ii < nxyz ; ii++ ) gfim[ii] = TANH( fac * gfim[ii] ) ;
      break ;
   }

   /* 29 Aug 1996: clean up the merged threshold, too */

   switch( MRG_cflag_gthr ){
      default: break ;

      case THFLAG_FICO:
         fac = 1.0 / num_dset ;
         for( ii=0 ; ii < nxyz ; ii++ ) ggfim[ii] = TANH( fac * ggfim[ii] ) ;
      break ;
   }

   /**** at this point, don't need the count brick "gnum" anymore;
         "gfim" contains the results we'll eventually write to disk ****/

   myRwcFree( gnum ) ;

   /*** if desired, edit the result for cluster size ***/

   rmm   = MRG_clust_rmm_g ;
   vmul  = MRG_clust_vmul_g ;
   dxyz  = dx*dy*dz ;
   ptmin = vmul / dxyz + 0.99 ;

   if( (rmm >= dx || rmm >= dy || rmm >= dz) && ptmin > 1 ){
      if( ! MRG_be_quiet ) printf("-- editing merger for cluster size\n") ;

      clar  = MCW_find_clusters( nx,ny,nz , dx,dy,dz , MRI_float,gfim , rmm) ;
      nclu  = 0 ;
      if( clar != NULL ){
         for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
            if( clar->clar[iclu] != NULL && clar->clar[iclu]->num_pt < ptmin ){
               KILL_CLUSTER(clar->clar[iclu]) ;
            } else if( clar->clar[iclu] != NULL ){
               nclu++ ;
            }
         }
      }

      if( nclu > 0 ){
         for( iclu=0 ; iclu < clar->num_clu ; iclu++ ){
            if( clar->clar[iclu] != NULL && clar->clar[iclu]->num_pt > 0 )
               MCW_cluster_to_vol( nx,ny,nz , MRI_float,gfim , clar->clar[iclu] ) ;
         }
      }
   }

   /*** scan results for non-zero-ositifulness ***/

   if( !MRG_doall ){  /* Nov 1998 */
      for( ii=0 ; ii < nxyz ; ii++ ) if( gfim[ii] != 0 ) break ;
      if( ii == nxyz ){
         fprintf(stderr,
           "*** Merged dataset has no nonzero entries -- will not write\n" ) ;
         exit(0) ;
      }
   }

   /*** attach new data to output brick      ***/
   /*   [Nov 1998: changed from ival to ivout] */

   switch( output_datum ){

      default:
         fprintf(stderr,
                 "*** Fatal Error ***\n"
                 "*** Somehow ended up with output_datum = %d\n",output_datum) ;
      exit(1) ;

      case MRI_complex:{
         void * dfim ;
         dfim = (void *) RwcMalloc( sizeof(complex) * nxyz ) ;
         EDIT_coerce_type( nxyz , MRI_float,gfim , MRI_complex,dfim ) ;
         myRwcFree( gfim ) ;
         EDIT_substitute_brick( new_dset , ivout , MRI_complex , dfim ) ;
         DSET_BRICK_FACTOR(new_dset,ivout) = 0.0 ;
      }
      break ;

      case MRI_float:
         EDIT_substitute_brick( new_dset , ivout , MRI_float , gfim ) ;
         DSET_BRICK_FACTOR(new_dset,ivout) = 0.0 ;
      break ;

      case MRI_byte:
      case MRI_short:{
         void * dfim ;
         float gtop ;

         gtop = MCW_vol_amax( nx,ny,nz , MRI_float,gfim ) ;

         if( MRG_cflag_g == CFLAG_FISHER ){
            fimfac = FUNC_COR_SCALE_SHORT ;
         } else if( gtop == 0.0 ||           /* 08 Jan 1998 */
                    MRG_nscale  ||           /* 15 Sep 2000 */
                    (is_int && gtop <= MRI_TYPE_maxval[output_datum]) ){
            fimfac = 0.0 ;
         } else {
            fimfac = MRI_TYPE_maxval[output_datum] / gtop ;
         }

         dfim = (void *) RwcMalloc( mri_datum_size(output_datum) * nxyz ) ;
         EDIT_coerce_scale_type( nxyz,fimfac , MRI_float,gfim , output_datum,dfim ) ;
         myRwcFree( gfim ) ;
         EDIT_substitute_brick( new_dset , ivout , output_datum , dfim ) ;
         DSET_BRICK_FACTOR(new_dset,ivout) = (fimfac != 0.0) ? 1.0/fimfac : 0.0 ;
      }
      break ;
   }

   /** 29 Aug 1996: attach output threshold, if any **/

   if( MRG_cflag_gthr != THFLAG_NONE ){
      short * dfim ;

      dfim   = (short *) RwcMalloc( sizeof(short) * nxyz ) ;
      thrfac = FUNC_COR_SCALE_SHORT ;
      EDIT_coerce_scale_type( nxyz,thrfac , MRI_float,ggfim , MRI_short,dfim ) ;
      myRwcFree( ggfim ) ;
      EDIT_substitute_brick( new_dset , DSET_THRESH_VALUE(new_dset) ,
                             MRI_short , dfim ) ;
      DSET_BRICK_FACTOR(new_dset,1) = 1.0/thrfac ;

      /* if all datasets were fico, then output is fico,
         and needs to get the degrees-of-freedom parameters stataux */

      if( num_fico == num_dset ){
         (void) EDIT_dset_items( new_dset , ADN_stat_aux,thr_stataux , ADN_none ) ;

      /* some datasets were fith, so output is fith */

      } else {
          EDIT_dset_items( new_dset , ADN_func_type,FUNC_THR_TYPE , ADN_none ) ;
      }
   }

   ivout++ ;

 }  /* iv    End of iteration over sub-bricks */


   /*** write to disk!!! ***/

   if( MRG_nozero && THD_count_nonzero_bricks(new_dset) == 0 ){ /* 17 Jan 2017 */
     INFO_message("output dataset is all zero ==> not writing because of -allzero") ;
     exit(0) ;
   }

   tross_Make_History( "3dmerge" , argc , argv , new_dset ) ;
   THD_load_statistics( new_dset ) ;
   THD_write_3dim_dataset( NULL,NULL , new_dset , True ) ;
   if( ! MRG_be_quiet )
     fprintf(stderr,"-- Wrote merged dataset: %s\n" , DSET_BRIKNAME(new_dset) ) ;
   exit(0) ;
}
