/*****************************************************************************
   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.
******************************************************************************/

#include "mrilib.h"

/*-- these macros stolen from file thd.h --*/

#define ORCODE(aa) \
  ( (aa)=='R' ? ORI_R2L_TYPE : (aa)=='L' ? ORI_L2R_TYPE : \
    (aa)=='P' ? ORI_P2A_TYPE : (aa)=='A' ? ORI_A2P_TYPE : \
    (aa)=='I' ? ORI_I2S_TYPE : (aa)=='S' ? ORI_S2I_TYPE : ILLEGAL_TYPE )

#define OR3OK(x,y,z) ( ((x)&6) + ((y)&6) + ((z)&6) == 6 )

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

void Syntax(char * msg)
{
   if( msg != NULL ){
      fprintf(stderr,"*** %s\n",msg) ; exit(1) ;
   }

   printf(
    "Usage: 3dUndump [options] infile ...\n"
    "Assembles a 3D dataset from an ASCII list of coordinates and\n"
    "(optionally) values.\n"
    "\n"
    "Options:\n"
    "  -prefix ppp  = 'ppp' is the prefix for the output dataset\n"
    "                   [default = undump].\n"
    "  -master mmm  = 'mmm' is the master dataset, whose geometry\n"
    "    *OR*           will determine the geometry of the output.\n"
    "  -dimen I J K = Sets the dimensions of the output dataset to\n"
    "                   be I by J by K voxels.  (Each I, J, and K\n"
    "                   must be >= 2.)  This option can be used to\n"
    "                   create a dataset of a specific size for test\n"
    "                   purposes, when no suitable master exists.\n"
    "          ** N.B.: Exactly one of -master or -dimen must be given.\n"
    "  -mask kkk    = This option specifies a mask dataset 'kkk', which\n"
    "                   will control which voxels are allowed to get\n"
    "                   values set.  If the mask is present, only\n"
    "                   voxels that are nonzero in the mask can be\n"
    "                   set in the new dataset.\n"
    "                   * A mask can be created with program 3dAutomask.\n"
    "                   * Combining a mask with sphere insertion makes\n"
    "                     a lot of sense (to me, at least).\n"
    "  -datum type  = 'type' determines the voxel data type of the\n"
    "                   output, which may be byte, short, or float\n"
    "                   [default = short].\n"
    "  -dval vvv    = 'vvv' is the default value stored in each\n"
    "                   input voxel that does not have a value\n"
    "                   supplied in the input file [default = 1].\n"
    "  -fval fff    = 'fff' is the fill value, used for each voxel\n"
    "                   in the output dataset that is NOT listed\n"
    "                   in the input file [default = 0].\n"
    "  -ijk         = Coordinates in the input file are (i,j,k) index\n"
    "       *OR*        triples, as might be output by 3dmaskdump.\n"
    "  -xyz         = Coordinates in the input file are (x,y,z)\n"
    "                   spatial coordinates, in mm.  If neither\n"
    "                   -ijk or -xyz is given, the default is -ijk.\n"
    "          ** N.B.: -xyz can only be used with -master. If -dimen\n"
    "                   is used to specify the size of the output dataset,\n"
    "                   (x,y,z) coordinates are not defined (until you\n"
    "                   use 3drefit to define the spatial structure).\n"
    "  -srad rrr    = Specifies that a sphere of radius 'rrr' will be\n"
    "                   filled about each input (x,y,z) or (i,j,k) voxel.\n"
    "                   If the radius is not given, or is 0, then each\n"
    "                   input data line sets the value in only one voxel.\n"
    "                   * If '-master' is used, then 'rrr' is in mm.\n"
    "                   * If '-dimen' is used, then 'rrr' is in voxels.\n"
    "  -orient code = Specifies the coordinate order used by -xyz.\n"
    "                   The code must be 3 letters, one each from the pairs\n"
    "                   {R,L} {A,P} {I,S}.  The first letter gives the\n"
    "                   orientation of the x-axis, the second the orientation\n"
    "                   of the y-axis, the third the z-axis:\n"
    "                     R = right-to-left         L = left-to-right\n"
    "                     A = anterior-to-posterior P = posterior-to-anterior\n"
    "                     I = inferior-to-superior  S = superior-to-inferior\n"
    "                   If -orient isn't used, then the coordinate order of the\n"
    "                   -master dataset is used to interpret (x,y,z) inputs.\n"
    "          ** N.B.: If -dimen is used (which implies -ijk), then the\n"
    "                   only use of -orient is to specify the axes ordering\n"
    "                   of the output dataset.  If -master is used instead,\n"
    "                   the output dataset's axes ordering is the same as the\n"
    "                   -master dataset's, regardless of -orient.\n"
    "  -head_only   =  A secret option for creating only the .HEAD file which\n"
    "                  gets exploited by the AFNI matlab library function\n"
    "                  New_HEAD.m\n"
    "\n"
    "Input File Format:\n"
    " The input file(s) are ASCII files, with one voxel specification per\n"
    " line.  A voxel specification is 3 numbers (-ijk or -xyz coordinates),\n"
    " with an optional 4th number giving the voxel value.  For example:\n"
    "\n"
    "   1 2 3 \n"
    "   3 2 1 5\n"
    "   5.3 6.2 3.7\n"
    "   // this line illustrates a comment\n"
    "\n"
    " The first line puts a voxel (with value given by -dval) at point\n"
    " (1,2,3).  The second line puts a voxel (with value 5) at point (3,2,1).\n"
    " The third line puts a voxel (with value given by -dval) at point\n"
    " (5.3,6.2,3.7).  If -ijk is in effect, and fractional coordinates\n"
    " are given, they will be rounded to the nearest integers; for example,\n"
    " the third line would be equivalent to (i,j,k) = (5,6,4).\n"
    "\n"
    "Notes:\n"
    "* This program creates a 1 sub-brick file.  You can 'glue' multiple\n"
    "   files together using 3dbucket or 3dTcat to make multi-brick datasets.\n"
    "* If an input filename is '-', then stdin is used for input.\n"
    "* By default, the output dataset is of type '-fim', unless the -master\n"
    "   dataset is an anat type. You can change the output type using 3drefit.\n"
    "* You could use program 1dcat to extract specific columns from a\n"
    "   multi-column rectangular file (e.g., to get a specific sub-brick\n"
    "   from the output of 3dmaskdump), and use the output of 1dcat as input\n"
    "   to this program.\n"
    "* [19 Feb 2004] The -mask and -srad options were added this day.\n"
    "   Also, a fifth value on an input line, if present, is taken as a\n"
    "   sphere radius to be used for that input point only.  Thus, input\n"
    "      3.3 4.4 5.5 6.6 7.7\n"
    "   means to put the value 6.6 into a sphere of radius 7.7 mm centered\n"
    "   about (x,y,z)=(3.3,4.4,5.5).\n"
    "\n"
    "-- RWCox -- October 2000\n"
   ) ;

   exit(0) ;
}

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

#define NBUF 1024  /* line buffer size */

int main( int argc , char * argv[] )
{
   int do_ijk=1 , dimen_ii=0 , dimen_jj=0 , dimen_kk=0 , datum=MRI_short ;
   THD_3dim_dataset *mset=NULL ;
   char *prefix="undump" , *orcode=NULL ;
   THD_coorder cord ;
   float dval_float=1.0 , fval_float=0.0 , *fbr=NULL ;
   short dval_short=1   , fval_short=0   , *sbr=NULL ;
   byte  dval_byte =1   , fval_byte =0   , *bbr=NULL , *mmask=NULL ;

   FILE *fp ;
   THD_3dim_dataset *dset , *maskset=NULL ;
   int iarg , ii,jj,kk,ll,ijk , nx,ny,nz , nxyz , nn , do_head_only=0;
   float      xx,yy,zz,vv=0.0 ;
   short               sv=0   ;
   byte                bv=0   ;
   char linbuf[NBUF] , *cp ;

   float xxdown,xxup , yydown,yyup , zzdown,zzup ;

   float srad=0.0 , vrad,rii,rjj,rkk,qii,qjj,qkk , dx,dy,dz ;  /* 19 Feb 2004 */

   /*-- help? --*/

   if( argc < 3 || strcmp(argv[1],"-help") == 0 ) Syntax(NULL) ;

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

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

   machdep() ;
   { 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("3dUndump",argc,argv) ;

   /*-- command line options --*/

   do_head_only = 0 ;
   iarg = 1 ;
   while( iarg < argc && argv[iarg][0] == '-' ){

      if( strcmp(argv[iarg],"-") == 0 ){    /* a single - is an input filename */
         break ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-prefix") == 0 ){
         if( iarg+1 >= argc )
            Syntax("-prefix: no argument follows!?") ;
         else if( !THD_filename_ok(argv[++iarg]) )
            Syntax("-prefix: Illegal prefix given!") ;
         prefix = argv[iarg] ;
         iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-master") == 0 ){
         if( iarg+1 >= argc )
            Syntax("-master: no argument follows!?") ;
         else if( mset != NULL )
            Syntax("-master: can't have two -master options!") ;
         else if( dimen_ii > 0 )
            Syntax("-master: conflicts with previous -dimen!") ;

         mset = THD_open_dataset( argv[++iarg] ) ;
         if( mset == NULL )
            Syntax("-master: can't open dataset" ) ;

         iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-mask") == 0 ){
        if( iarg+1 >= argc )
          Syntax("-mask: no argument follows!?") ;
        else if( maskset != NULL )
          Syntax("-mask: can't have two -mask options!") ;

        maskset = THD_open_dataset( argv[++iarg] ) ;
        if( maskset == NULL )
          Syntax("-mask: can't open dataset" ) ;

        iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-dimen") == 0 ){
         if( iarg+3 >= argc )
            Syntax("-dimen: don't have 3 arguments following!?") ;
         else if( mset != NULL )
            Syntax("-dimen: conflicts with previous -master!") ;
         else if( dimen_ii > 0 )
            Syntax("-dimen: can't have two -dimen options!") ;
         dimen_ii = strtol(argv[++iarg],NULL,10) ;
         dimen_jj = strtol(argv[++iarg],NULL,10) ;
         dimen_kk = strtol(argv[++iarg],NULL,10) ;
         if( dimen_ii < 2 || dimen_jj < 2 || dimen_kk < 2 )
            Syntax("-dimen: values following are not all >= 2!") ;

         iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-datum") == 0 ){
         if( ++iarg >= argc )
            Syntax("-datum: no argument follows?!") ;

         if( strcmp(argv[iarg],"short") == 0 )
            datum = MRI_short ;
         else if( strcmp(argv[iarg],"float") == 0 )
            datum = MRI_float ;
         else if( strcmp(argv[iarg],"byte") == 0 )
            datum = MRI_byte ;
         else
            Syntax("-datum: illegal type given!") ;

         iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-srad") == 0 ){   /* 19 Feb 2004 */
        if( iarg+1 >= argc )
          Syntax("-srad: no argument follows!?") ;

        srad = strtod( argv[++iarg] , NULL ) ;
        if( srad <= 0.0 ){
          fprintf(stderr,"++ WARNING: -srad value of %g is ignored!\n",srad);
          srad = 0.0 ;
        }
        iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-dval") == 0 ){
         if( iarg+1 >= argc )
            Syntax("-dval: no argument follows!?") ;

         dval_float = strtod( argv[++iarg] , NULL ) ;
         dval_short = (short) rint(dval_float) ;
         dval_byte  = (byte)  dval_short ;
         iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-fval") == 0 ){
         if( iarg+1 >= argc )
            Syntax("-fval: no argument follows!?") ;

         fval_float = strtod( argv[++iarg] , NULL ) ;
         fval_short = (short) rint(fval_float) ;
         fval_byte  = (byte)  fval_short ;
         iarg++ ; continue ;
      }

      if( strcmp(argv[iarg],"-ijk") == 0 ){
         do_ijk = 1 ;
         iarg++ ; continue ;
      }

      /*-----*/

      if( strcmp(argv[iarg],"-xyz") == 0 ){
         do_ijk = 0 ;
         iarg++ ; continue ;
      }
      
      /*-----*/

      if( strcmp(argv[iarg],"-head_only") == 0 ){
         do_head_only = 1 ;
         iarg++ ; continue ;
      }

      /*-----*/
      if( strcmp(argv[iarg],"-orient") == 0 ){
         int xx,yy,zz ;
         if( iarg+1 >= argc )
            Syntax("-orient: no argument follows!?") ;

         orcode = argv[++iarg] ;
         if( strlen(orcode) != 3 )
            Syntax("-orient: illegal argument follows") ;

         xx = ORCODE(orcode[0]) ; yy = ORCODE(orcode[1]) ; zz = ORCODE(orcode[2]) ;
         if( xx < 0 || yy < 0 || zz < 0 || !OR3OK(xx,yy,zz) )
            Syntax("-orient: illegal argument follows") ;

         iarg++ ; continue ;
      }

      /*-----*/

      fprintf(stderr,"*** Unknown option: %s\n",argv[iarg]) ; exit(1) ;

   } /* end of loop over command line options */

   /*-- check for inconsistencies --*/

   if( iarg >= argc && !do_head_only)
      Syntax("No input files on command line!?") ;

   if( do_ijk == 0 && mset == NULL )
      Syntax("Can't use -xyz without -master also!") ;

   if( mset == NULL && dimen_ii < 2 )
      Syntax("Must use exactly one of -master or -dimen options on command line");

   if( (datum == MRI_short && dval_short == fval_short) ||
       (datum == MRI_float && dval_float == fval_float) ||
       (datum == MRI_byte  && dval_byte  == fval_byte )   ){

      fprintf(stderr,"+++ Warning: -dval and -fval are the same!\n") ;
   }

   /*-- set orcode to value from -master, if this is needed --*/

   if( mset != NULL && do_ijk == 0 && orcode == NULL ){
      orcode = malloc(4) ;
      orcode[0] = ORIENT_typestr[mset->daxes->xxorient][0] ;
      orcode[1] = ORIENT_typestr[mset->daxes->yyorient][0] ;
      orcode[2] = ORIENT_typestr[mset->daxes->zzorient][0] ;
      orcode[3] = '\0' ;
   }

   THD_coorder_fill( orcode , &cord ) ;  /* setup coordinate order */

   /*-- make empty dataset --*/

   if( mset != NULL ){                 /* from -master */

      dset = EDIT_empty_copy( mset ) ;
      EDIT_dset_items( dset ,
                          ADN_prefix    , prefix ,
                          ADN_datum_all , datum ,
                          ADN_nvals     , 1 ,
                          ADN_ntt       , 0 ,
                          ADN_func_type , ISANAT(mset) ? mset->func_type
                                                       : FUNC_FIM_TYPE ,

                          ADN_directory_name , "./" ,
                       ADN_none ) ;

   } else {                            /* from nothing */

     THD_ivec3 iv_nxyz   , iv_xyzorient ;
     THD_fvec3 fv_xyzorg , fv_xyzdel ;

     LOAD_IVEC3( iv_nxyz , dimen_ii , dimen_jj , dimen_kk ) ;
     LOAD_IVEC3( iv_xyzorient , cord.xxor , cord.yyor , cord.zzor ) ;
     LOAD_FVEC3( fv_xyzdel ,
                 ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? 1.0 : -1.0 ,
                 ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? 1.0 : -1.0 ,
                 ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? 1.0 : -1.0  ) ;
     LOAD_FVEC3( fv_xyzorg ,
                 ORIENT_sign[iv_xyzorient.ijk[0]]=='+' ? -0.5*dimen_ii : 0.5*dimen_ii,
                 ORIENT_sign[iv_xyzorient.ijk[1]]=='+' ? -0.5*dimen_jj : 0.5*dimen_jj,
                 ORIENT_sign[iv_xyzorient.ijk[2]]=='+' ? -0.5*dimen_kk : 0.5*dimen_kk ) ;

     dset = EDIT_empty_copy( NULL ) ;

     EDIT_dset_items( dset ,
                       ADN_nxyz      , iv_nxyz ,
                       ADN_xyzdel    , fv_xyzdel ,
                       ADN_xyzorg    , fv_xyzorg ,
                       ADN_xyzorient , iv_xyzorient ,
                       ADN_prefix    , prefix ,
                       ADN_datum_all , datum ,
                       ADN_nvals     , 1 ,
                       ADN_ntt       , 0 ,
                       ADN_type      , HEAD_FUNC_TYPE ,
                       ADN_func_type , FUNC_FIM_TYPE ,
                    ADN_none ) ;
   }

   if( THD_deathcon() && THD_is_file(DSET_HEADNAME(dset)) )
      Syntax("Output dataset already exists -- can't overwrite") ;

   if (do_head_only) {
      DSET_write_header(dset);
      exit(0);
   }
   
   /*-- make empty brick array for dataset --*/

   EDIT_substitute_brick( dset , 0 , datum , NULL ) ;  /* will make array */

   nx = DSET_NX(dset); ny = DSET_NY(dset); nz = DSET_NZ(dset); nxyz = nx*ny*nz;

   /* 19 Feb 2004: check and make mask if desired */

   if( maskset != NULL &&
       ( DSET_NX(maskset) != nx ||
         DSET_NY(maskset) != ny ||
         DSET_NZ(maskset) != nz   ) )
     Syntax("-mask dataset doesn't match dimension of output dataset") ;

   if( maskset != NULL ){
     mmask = THD_makemask( maskset , 0 , 1.0,-1.0 ) ;
     if( mmask == NULL ){
       fprintf(stderr,"++ WARNING: can't create mask for some reason!\n") ;
     } else {
       int nmask = THD_countmask( nxyz , mmask ) ;
       if( nmask == 0 ){
         fprintf(stderr,"++ WARNING: 0 voxels in mask -- ignoring it!\n") ;
         free((void *)mmask) ; mmask = NULL ;
       } else {
         fprintf(stderr,"++ %d voxels found in mask\n",nmask) ;
       }
     }
     DSET_delete(maskset) ;
   }

   /*-- fill new dataset brick with the -fval value --*/

   switch( datum ){
      case MRI_short:
         sbr = (short *) DSET_BRICK_ARRAY(dset,0) ;
         for( ii=0 ; ii < nxyz ; ii++ ) sbr[ii] = fval_short ;
      break ;

      case MRI_float:
         fbr = (float *) DSET_BRICK_ARRAY(dset,0) ;
         for( ii=0 ; ii < nxyz ; ii++ ) fbr[ii] = fval_float ;
      break ;

      case MRI_byte:
         bbr = (byte *) DSET_BRICK_ARRAY(dset,0) ;
         for( ii=0 ; ii < nxyz ; ii++ ) bbr[ii] = fval_byte ;
      break ;
   }

   /* 24 Nov 2000: get the bounding box for the dataset */

   dx = fabs(dset->daxes->xxdel) ; if( dx <= 0.0 ) dx = 1.0 ;
   dy = fabs(dset->daxes->yydel) ; if( dy <= 0.0 ) dy = 1.0 ;
   dz = fabs(dset->daxes->zzdel) ; if( dz <= 0.0 ) dz = 1.0 ;

   if( !do_ijk ){
#ifndef EXTEND_BBOX
      xxdown = dset->daxes->xxmin - 0.501 * dx ;
      xxup   = dset->daxes->xxmax + 0.501 * dx ;
      yydown = dset->daxes->yymin - 0.501 * dy ;
      yyup   = dset->daxes->yymax + 0.501 * dy ;
      zzdown = dset->daxes->zzmin - 0.501 * dz ;
      zzup   = dset->daxes->zzmax + 0.501 * dz ;
#else
      xxdown = dset->daxes->xxmin ;
      xxup   = dset->daxes->xxmax ;
      yydown = dset->daxes->yymin ;
      yyup   = dset->daxes->yymax ;
      zzdown = dset->daxes->zzmin ;
      zzup   = dset->daxes->zzmax ;
#endif
   }

   /*-- loop over input files and read them line by line --*/

   for( ; iarg < argc ; iarg++ ){  /* iarg is already set at start of this loop */

      /* get input file ready to read */

      if( strcmp(argv[iarg],"-") == 0 ){  /* stdin */
         fp = stdin ;
      } else {                            /* OK, open the damn file */
         fp = fopen( argv[iarg] , "r" ) ;
         if( fp == NULL ){
            fprintf(stderr,
                    "+++ Warning: can't open input file %s -- skipping it\n",
                    argv[iarg]) ;
            continue ;                    /* skip to end of iarg loop */
         }
      }

      /* read lines, process and store */

      ll = 0 ;
      while(1){
         ll++ ;                               /* line count */
         cp = fgets( linbuf , NBUF , fp ) ;
         if( cp == NULL ) break ;             /* end of file => end of loop */
         kk = strlen(linbuf) ;
         if( kk == 0 ) continue ;             /* empty line => get next line */

         /* find 1st nonblank */

         for( ii=0 ; ii < kk && isspace(linbuf[ii]) ; ii++ ) ; /* nada */
         if( ii == kk ) continue ;                                 /* all blanks */
         if( linbuf[ii] == '/' && linbuf[ii+1] == '/' ) continue ; /* comment */
         if( linbuf[ii] == '#'                        ) continue ; /* comment */

         /* scan line for data */

         vv   = dval_float ;   /* if not scanned in below, use the default value */
         vrad = srad ;         /* 19 Feb 2004: default sphere radius */
         nn   = sscanf(linbuf+ii , "%f%f%f%f%f" , &xx,&yy,&zz,&vv,&vrad ) ;
         if( nn < 3 ){
           fprintf(stderr,"+++ Warning: file %s line %d: incomplete\n",argv[iarg],ll) ;
           continue ;
         }

         /* get voxel index into (ii,jj,kk) */

         if( do_ijk ){   /* inputs are (ii,jj,kk) themselves */

            ii = (int) rint(xx) ; jj = (int) rint(yy) ; kk = (int) rint(zz) ;
            if( ii < 0 || ii >= nx ){
               fprintf(stderr,
                       "+++ Warning: file %s line %d: i index=%d is invalid\n",
                       argv[iarg],ll,ii) ;
               continue ;
            }
            if( jj < 0 || jj >= ny ){
               fprintf(stderr,
                       "+++ Warning: file %s line %d: j index=%d is invalid\n",
                       argv[iarg],ll,jj) ;
               continue ;
            }
            if( kk < 0 || kk >= nz ){
               fprintf(stderr,
                       "+++ Warning: file %s line %d: k index=%d is invalid\n",
                       argv[iarg],ll,kk) ;
               continue ;
            }

         } else {  /* inputs are coordinates => must convert to index */

            THD_fvec3 mv , dv ;                              /* temp vectors */
            THD_ivec3 iv ;

            THD_coorder_to_dicom( &cord , &xx,&yy,&zz ) ;    /* to Dicom order */
            LOAD_FVEC3( dv , xx,yy,zz ) ;
            mv = THD_dicomm_to_3dmm( dset , dv ) ;           /* to Dataset order */

            /* 24 Nov 2000: check (xx,yy,zz) for being inside the box */

            if( mv.xyz[0] < xxdown || mv.xyz[0] > xxup ){
               fprintf(stderr,"+++ Warning: file %s line %d: x coord=%g is outside %g .. %g\n" ,
                       argv[iarg],ll,mv.xyz[0] , xxdown,xxup ) ;
               continue ;
            }
            if( mv.xyz[1] < yydown || mv.xyz[1] > yyup ){
               fprintf(stderr,"+++ Warning: file %s line %d: y coord=%g is outside %g .. %g\n" ,
                       argv[iarg],ll,mv.xyz[1] , yydown , yyup ) ;
               continue ;
            }
            if( mv.xyz[2] < zzdown || mv.xyz[2] > zzup ){
               fprintf(stderr,"+++ Warning: file %s line %d: z coord=%g is outside %g .. %g\n" ,
                       argv[iarg],ll,mv.xyz[2] , zzdown , zzup ) ;
               continue ;
            }

            iv = THD_3dmm_to_3dind( dset , mv ) ;            /* to Dataset index */
            ii = iv.ijk[0]; jj = iv.ijk[1]; kk = iv.ijk[2];  /* save */
         }

         /* now load individual voxel (ii,jj,kk) */

         ijk = ii + jj*nx + kk*nx*ny ;
         if( mmask == NULL || mmask[ijk] ){
          switch( datum ){
            case MRI_float:{
              if( fbr[ijk] != fval_float && fbr[ijk] != vv )
                fprintf(stderr,"Overwrite voxel %d %d %d\n",ii,jj,kk) ;
              fbr[ijk] = vv ;
            }
            break ;
            case MRI_short:{
              sv = SHORTIZE(vv) ;
              if( sbr[ijk] != fval_short && sbr[ijk] != sv )
                fprintf(stderr,"Overwrite voxel %d %d %d\n",ii,jj,kk) ;
              sbr[ijk] = sv ;
            }
            break ;
            case MRI_byte:{
              bv = BYTEIZE(vv) ;
              if( bbr[ijk] != fval_byte && bbr[ijk] != bv )
                fprintf(stderr,"Overwrite voxel %d %d %d\n",ii,jj,kk) ;
              bbr[ijk] = bv ;
            }
            break ;
          }
         }

         /* 19 Feb 2004:
             Make up radius of ellipsoid in voxel indexes
             Will put value vv into all voxels (aa,bb,cc) with
               (aa-ii)**2/qii + (bb-jj)**2/qjj + (cc-kk)**2/qkk <= 1.0 */

         vrad *= 1.00001 ;                               /* expand a little */
         rii = vrad/dx ; rjj = vrad/dy ; rkk = vrad/dz ;
         qii = rii*rii ; qjj = rjj*rjj ; qkk = rkk*rkk ;

         if( rii >= 1.0 || rjj >= 1.0 || rkk >= 1.0 ){
           int aa,bb,cc , abot,atop,bbot,btop,cbot,ctop; float rr;
           abot = ii-(int)rint(rii) ; atop = ii+(int)rint(rii) ;
           if( abot < 0 ) abot = 0 ; if( atop >= nx ) atop = nx-1 ;
           bbot = jj-(int)rint(rjj) ; btop = jj+(int)rint(rjj) ;
           if( bbot < 0 ) bbot = 0 ; if( btop >= ny ) btop = ny-1 ;
           cbot = kk-(int)rint(rkk) ; ctop = kk+(int)rint(rkk) ;
           if( cbot < 0 ) cbot = 0 ; if( ctop >= nz ) ctop = nz-1 ;
           for( cc=cbot ; cc <= ctop ; cc++ ){
             for( bb=bbot ; bb <= btop ; bb++ ){
               for( aa=abot ; aa <= atop ; aa++ ){
                 rr =  (aa-ii)*(aa-ii)/qii
                     + (bb-jj)*(bb-jj)/qjj + (cc-kk)*(cc-kk)/qkk ;
                 if( rr <= 1.00001 ){
                   ijk = aa + bb*nx + cc*nx*ny ;    /* (aa,bb,cc) in dataset */
                   if( mmask == NULL || mmask[ijk] ){
                     switch( datum ){
                       case MRI_float: fbr[ijk] = vv ; break ;
                       case MRI_short: sbr[ijk] = sv ; break ;
                       case MRI_byte:  bbr[ijk] = bv ; break ;
                     }
                   }
                 }
               }
             }
           }
         }  /* 19 Feb 2004: end of inserting a sphere */

      } /* end of loop over input lines */

      /* close input file */

      if( fp != stdin ) fclose( fp ) ;

   } /* end of loop over input files */

   tross_Make_History( "3dUndump" , argc,argv , dset ) ;
   DSET_write(dset) ;
   fprintf(stderr,"+++ Wrote results to dataset %s\n",DSET_BRIKNAME(dset)) ;
   exit(0) ;
}