/*
 * for each input dataset name
 *    open (check dims, etc.)
 *    dilate (zeropad, make binary, dilate, unpad, apply)
 *    fill list of bytemask datasets
 *
 * also, count total volumes
 */
int process_input_dsets(param_t * params)
{
   THD_3dim_dataset * dset, * dfirst=NULL;
   int                iset, nxyz;

   ENTRY("process_input_dsets");

   if( !params ) ERROR_exit("NULL inputs to PID");

   if( params->ndsets <= 0 ) {
      ERROR_message("process_input_dsets: no input datasets");
      RETURN(1);
   }

   /* allocate space for dsets array */
   params->dsets = (THD_3dim_dataset **)malloc(params->ndsets*
                                               sizeof(THD_3dim_dataset*));
   if( !params->dsets ) ERROR_exit("failed to allocate dset pointers");

   if( params->verb ) INFO_message("processing %d input datasets...",
                                   params->ndsets);
   
   /* warn user of dilations */
   if(params->verb && params->ndsets) {
      int pad = needed_padding(&params->IND);
      INFO_message("padding all datasets by %d (for dilations)", pad);
   }

   /* process the datasets */
   nxyz = 0;
   for( iset=0; iset < params->ndsets; iset++ ) {
      /* open and verify dataset */
      dset = THD_open_dataset(params->inputs[iset]);
      if( !dset ) ERROR_exit("failed to open mask dataset '%s'",
                             params->inputs[iset]);
      DSET_load(dset);  CHECK_LOAD_ERROR(dset);

      if( params->verb>1 ) INFO_message("loaded dset %s, with %d volumes",
                                        DSET_PREFIX(dset), DSET_NVALS(dset));

      if( nxyz == 0 ) { /* make an empty copy of the first dataset */
         nxyz = DSET_NVOX(dset);
         dfirst = EDIT_empty_copy(dset);
      }

      /* check for consistency in voxels and grid */
      if( DSET_NVOX(dset) != nxyz ) ERROR_exit("nvoxel mis-match");
      if( ! EQUIV_GRIDS(dset, dfirst) )
         WARNING_message("grid from dset %s does not match that of dset %s",
                         DSET_PREFIX(dset), DSET_PREFIX(dfirst));

      /* apply dilations to all volumes, returning bytemask datasets */
      params->dsets[iset] = apply_dilations(dset, &params->IND,1,params->verb);
      if( ! params->dsets[iset] ) RETURN(1);
   } 

   DSET_delete(dfirst); /* and nuke */

   RETURN(0);
}

THD_3dim_dataset *Seg_load_dset_eng( char *set_name, char *view ) 
{
   static char FuncName[]={"Seg_load_dset_eng"};
   THD_3dim_dataset *dset=NULL, *sdset=NULL;
   int i=0;
   byte make_cp=0;
   int verb=0;
   char sprefix[THD_MAX_PREFIX+10], *stmp=NULL;
   
   SUMA_ENTRY;
   
   dset = THD_open_dataset( set_name );
   if( !ISVALID_DSET(dset) ){
     fprintf(stderr,"**ERROR: can't open dataset %s\n",set_name) ;
     SUMA_RETURN(NULL);
   }
   
   DSET_mallocize(dset)   ; DSET_load(dset);
   
   for (i=0; i<DSET_NVALS(dset); ++i) {
      if (DSET_BRICK_TYPE(dset,i) != MRI_short) {
         if (verb) INFO_message("Sub-brick %d in %s not of type short.\n"
                       "Creating new short copy of dset ", 
                       i, DSET_PREFIX(dset));
         make_cp=1; break;
      }
   }
   
   if (make_cp) {
      if (!SUMA_ShortizeDset(&dset, -1.0)) {
         SUMA_S_Err("**ERROR: Failed to shortize");
         SUMA_RETURN(NULL);
      }
   }
   
   if (DSET_IS_MASTERED(dset)) {
      if (verb) INFO_message("Dset is mastered, making copy...");
      stmp = SUMA_ModifyName(set_name, "append", ".cp", NULL);
      sdset = dset;
      dset = EDIT_full_copy(sdset, stmp);
      free(stmp); DSET_delete(sdset); sdset = NULL;  
   }
      
   
   if (view) {
      if (view) {
               if (!strstr(view,"orig")) 
            EDIT_dset_items(dset,ADN_view_type, VIEW_ORIGINAL_TYPE, ADN_none); 
         else  if (!strstr(view,"acpc")) 
            EDIT_dset_items(dset,ADN_view_type, VIEW_ACPCALIGNED_TYPE, ADN_none);
         else  if (!strstr(view,"tlrc")) 
            EDIT_dset_items(dset ,ADN_view_type, VIEW_TALAIRACH_TYPE, ADN_none);
         else SUMA_S_Errv("View of %s is rubbish", view);
      }
   }
   
   SUMA_RETURN(dset);
}

SEXP R_SUMA_HistString (SEXP SCallingFunc, SEXP Sarg, SEXP Shold) {
   char *CallingFunc=NULL;
   SEXP Rname = R_NilValue;
   char *fname = NULL, *hold=NULL, **arg=NULL, *res=NULL;
   int debug=0, nprot=0, narg=0, i= 0;
   
   if (!debug) debug = get_odebug();
   if (isNull(SCallingFunc) || isNull(Sarg)) {
      ERROR_message("Null input to R_SUMA_HistString");
      return(Rname);
   }
   /* get the executable name */
   PROTECT(SCallingFunc = AS_CHARACTER(SCallingFunc)); ++nprot;
   fname = R_alloc(strlen(CHAR(STRING_ELT(SCallingFunc,0)))+1, sizeof(char));
   strcpy(fname, CHAR(STRING_ELT(SCallingFunc,0)));
   if (debug) INFO_message("filename %s\n", fname);
  
   /* get the arg */
   PROTECT(Sarg = AS_CHARACTER(Sarg)); ++nprot;
   narg = (LENGTH(Sarg));
   arg = (char **)calloc(narg+1, sizeof(char *));
   if (fname) arg[0] = strdup(fname);
   else arg[0] = strdup("UnChevalSansNom");
   for (i=1; i<=narg; ++i) {
      arg[i] = (char *)calloc(strlen(CHAR(STRING_ELT(Sarg,i-1)))+1, 
                              sizeof(char));
      strcpy(arg[i], CHAR(STRING_ELT(Sarg,i-1)));
      if (debug) INFO_message("arg %d/%d %s\t", i, narg, arg[i]);
   }
   
   /* any old history ? */
   if (!isNull(Shold)) {
      PROTECT(Shold = AS_CHARACTER(Shold)); ++nprot;
      hold = R_alloc(strlen(CHAR(STRING_ELT(Shold,0)))+1, sizeof(char));
      strcpy(hold, CHAR(STRING_ELT(Shold,0)));
      if (debug) INFO_message("hold %s\n", hold);
   }
    
   if (( res = SUMA_HistString (fname, narg+1, arg, hold))) {
      PROTECT(Rname = allocVector(STRSXP, 1)); ++nprot;
      SET_STRING_ELT(Rname, 0, mkChar(res)); 
      if (debug) INFO_message("hist is %s\n", res);
      SUMA_free(res); res=NULL;
   } else {
      ERROR_message("Call to SUMA_HistString %s failed", fname);
   }
   
   for (i=0; i<=narg; ++i) { if (arg[i]) free(arg[i]); } free(arg); arg=NULL;
   
   UNPROTECT(nprot); 
   return(Rname);
}

SEXP R_SUMA_ParseModifyName(SEXP Sfname, SEXP Swhat, SEXP Sval, SEXP Scwd) 
{
   SEXP Rname = R_NilValue;
   char *fname = NULL, *what=NULL, *val=NULL, *cwd=NULL, *res=NULL;
   int debug=0, nprot=0;
   
   if (!debug) debug = get_odebug();
   if (isNull(Sfname) || isNull(Swhat) || isNull(Sval)) {
      ERROR_message("Null input to R_SUMA_ModifyName");
      return(Rname);
   }
   /* get the filename */
   PROTECT(Sfname = AS_CHARACTER(Sfname)); ++nprot;
   fname = R_alloc(strlen(CHAR(STRING_ELT(Sfname,0)))+1, sizeof(char));
   strcpy(fname, CHAR(STRING_ELT(Sfname,0)));
   if (debug) INFO_message("filename %s\n", fname);
  
   /* get the what */
   PROTECT(Swhat = AS_CHARACTER(Swhat)); ++nprot;
   what = R_alloc(strlen(CHAR(STRING_ELT(Swhat,0)))+1, sizeof(char));
   strcpy(what, CHAR(STRING_ELT(Swhat,0)));
   if (debug) INFO_message("what %s\n", what);

   /* get the val */
   PROTECT(Sval = AS_CHARACTER(Sval)); ++nprot;
   val = R_alloc(strlen(CHAR(STRING_ELT(Sval,0)))+1, sizeof(char));
   strcpy(val, CHAR(STRING_ELT(Sval,0)));
   if (debug) INFO_message("val %s\n", val);

   /* get the cwd */
   if (!isNull(Scwd)) {
      PROTECT(Scwd = AS_CHARACTER(Scwd)); ++nprot;
      cwd = R_alloc(strlen(CHAR(STRING_ELT(Scwd,0)))+1, sizeof(char));
      strcpy(cwd, CHAR(STRING_ELT(Scwd,0)));
      if (debug) INFO_message("cwd %s\n", cwd);
   } 
   
   if (debug) INFO_message("Modifying %s\n", fname);
   if ((res = SUMA_ModifyName(fname, what, val, cwd))) {
      PROTECT(Rname = allocVector(STRSXP, 1)); ++nprot;
      SET_STRING_ELT(Rname, 0, mkChar(res)); 
      SUMA_free(res); res=NULL;
   } else {
      ERROR_message("Call to SUMA_ModifyName %s %s %s failed", fname, what, val);
   }
   UNPROTECT(nprot); 
   
   return(Rname);
}

int fill_mask(options_t * opts)
{
   THD_3dim_dataset * mset;
   int nvox;

ENTRY("fill_mask");

   if( opts->automask ) {
      if( opts->verb ) INFO_message("creating automask...");

      opts->mask = THD_automask(opts->inset);
      if( ! opts->mask ) {
         ERROR_message("failed to apply -automask");
         RETURN(1);
      }

      RETURN(0);
   }

   if( opts->mask_name ) {
      if( opts->verb )
         INFO_message("reading mask dset from %s...", opts->mask_name);

      mset = THD_open_dataset( opts->mask_name );
      if( ! mset ) ERROR_exit("cannot open mask dset '%s'", opts->mask_name);
      nvox = DSET_NVOX(opts->inset);
      if( DSET_NVOX(mset) != nvox ) {
         ERROR_message("mask does not have the same voxel count as input");
         RETURN(1);
      }

      /* fill mask array and mask_nxyz, remove mask dset */
      DSET_load(mset); CHECK_LOAD_ERROR(mset);

      opts->mask = THD_makemask(mset, 0, 1, 0);
      DSET_delete(mset);

      if( ! opts->mask ) {
         ERROR_message("cannot make mask from '%s'", opts->mask_name);
         RETURN(1);
      }

      if( opts->verb > 1 )
         INFO_message("have mask with %d voxels", nvox);
   }

   RETURN(0);
}

/*
 * A hole is defined as a connected set of zero voxels that does
 * not reach an edge.
 *
 * The core functionality was added to libmri.a in THD_mask_fill_holes.
 */
int fill_holes(THD_3dim_dataset * dset, int verb)
{
   short * sptr;     /* to for filling holes */
   byte  * bmask;    /* computed result */
   int     nfilled;
   int     nx, ny, nz, nvox, index, fill=0;

   ENTRY("fill_holes");

   bmask = THD_makemask(dset, 0, 1, 0); /* copy input as byte mask */
   nx = DSET_NX(dset);  ny = DSET_NY(dset);  nz = DSET_NZ(dset);
   nvox = DSET_NVOX(dset);

   /* created filled mask */
   nfilled = THD_mask_fill_holes(nx,ny,nz, bmask, verb);
   if( nfilled < 0 ) { ERROR_message("failed to fill holes");  RETURN(1); }

   /* apply to short volume */
   sptr = DBLK_ARRAY(dset->dblk, 0);
   for( index = 0; index < nvox; index++ )
      if( !sptr[index] && bmask[index] ) { fill++;  sptr[index] = 1; }

   if(verb>2) INFO_message("final check: fill=%d, nfilled=%d", fill, nfilled);

   RETURN(0);
}

int main( int argc , char *argv[] )
{
   THD_3dim_dataset *dset ; int aa,ll ; char *cpt ; float val ;

   if( argc < 2 ){
     printf("Usage: 3dSatCheck dataset [...]\n"
            "\n"
            "Prints the 'raw' initial transient (saturation) check\n"
            "value for each dataset on the command line.  Round this\n"
            "number to the nearest integer to get an estimate of\n"
            "how many non-saturated time points start a dataset.\n"
           ) ;
     exit(0) ;
   }
   for( aa=1 ; aa < argc ; aa++ ){
     dset = THD_open_dataset( argv[aa] ) ; if( !ISVALID_DSET(dset) ) continue ;
     if( DSET_NVALS(dset) < 9 ) continue ;
     DSET_load(dset) ; if( !DSET_LOADED(dset) ) continue ;
     val = THD_saturation_check( dset , NULL , 0,0 ) ;
     ll = strlen(argv[aa]) ;
     cpt = (ll <= 50) ? argv[aa] : argv[aa]+(ll-50) ;
     INFO_message("%-50.50s = %.3f",cpt,val) ;
     DSET_delete(dset) ;
   }
   exit(0) ;
}

void AFNI_coord_filer_setup( Three_D_View *im3d )
{
   char ename[32] , *eval ; int ic ;

ENTRY("AFNI_coord_filer_setup") ;

   if( !IM3D_OPEN(im3d) ) EXRETURN ;
   ic = AFNI_controller_index(im3d) ;
   if( ic < 0 || ic >= MAX_CONTROLLERS || fpc[ic] != NULL ) EXRETURN ;

   sprintf(ename,"AFNI_FILE_COORDS_%c",abet[ic]) ;
   eval = my_getenv(ename) ;
   if( eval == NULL || *eval == '\0' ){ fpc[ic] = NULL ; EXRETURN ; }

   if( strcmp(eval,"-") == 0 || strncmp(eval,"stdout",6) == 0 )
     fpc[ic] = stdout ;
   else {
     fpc[ic] = fopen( eval , "w" ) ;
     if( fpc[ic] == NULL ){
       ERROR_message("Unable to open file %s from %s",eval,ename) ;
       EXRETURN ;
     }
   }

   AFNI_receive_init( im3d , RECEIVE_VIEWPOINT_MASK ,
                             AFNI_filer_viewpoint_CB ,
                             im3d , "AFNI_filer_viewpoint_CB" ) ;

   INFO_message("Logging [%c] viewpoint changes to '%s'",abet[ic],eval) ;
   EXRETURN ;
}

int main( int argc , char *argv[] )
{
   int ndset=0 , ii ;
   THD_3dim_dataset **dset ;

   if( argc < 3 || strcasecmp(argv[1],"-help") == 0 ){
     printf(
       "** Program 3dConformist reads in a collection of datasets and\n"
       "   zero pads them to the same size.\n"
       "** The output volume size is the smallest region that includes\n"
       "   all datasets (i.e., the minimal covering box).\n"
       "** If the datasets cannot be processed (e.g., different grid\n"
       "   spacings), then nothing will happen except for error messages.\n"
       "** The purpose of this program is to be used in scripts that\n"
       "   process lots of datasets and needs to make them all conform\n"
       "   to the same size for collective voxel-wise analyses.\n"
       "** The input datasets ARE ALTERED (embiggened)! <<<<<<------******\n"
       "   Therefore, don't use this program casually.\n"
     ) ;
     exit(0) ;
   }

   mainENTRY("3dConformist") ; machdep() ; PRINT_VERSION("3dConformist") ;

   ndset = argc-1 ;
   dset  = (THD_3dim_dataset **)malloc(sizeof(THD_3dim_dataset *)*ndset) ;
   for( ii=0 ; ii < ndset ; ii++ ){
     dset[ii] = THD_open_dataset(argv[ii+1]) ;
     CHECK_OPEN_ERROR(dset[ii],argv[ii+1]) ;
   }

   ii = THD_conformist(ndset,dset,CONFORM_REWRITE,NULL) ;

   switch(ii){
     default: INFO_message ("3dConformist: Re-wrote %d datasets",ii) ; break ;
     case  0: INFO_message ("3dConformist: all datasets matched OK") ; break ;
     case -1: ERROR_message("3dConformist: bad input")               ; break ;
     case -2: ERROR_message("3dConformist: bad inputs")              ; break ;
     case -3: ERROR_message("3dConformist: can't match grids")       ; break ;
   }

   exit(0) ;
}

int write_results(options_t * opts)
{
   FILE * fp;

   int c;

   ENTRY("write_results");

   if( ! opts->result ) {
      ERROR_message("no results to write!");
      RETURN(1);
   }

   if( ! opts->prefix )                         fp = stdout;
   else if ( !strcmp(opts->prefix, "-") )       fp = stdout;
   else if ( !strcmp(opts->prefix, "stdout") )  fp = stdout;
   else if ( !strcmp(opts->prefix, "stderr") )  fp = stderr;
   else {
      fp = fopen(opts->prefix, "w");
      if( ! fp ) {
         ERROR_message("failed to open '%s' for writing", opts->prefix);
         RETURN(1);
      }
   }

   if( opts->verb ) {
      if     ( fp == stdout ) INFO_message("writing to stdout...");
      else if( fp == stdout ) INFO_message("writing to stderr...");
      else if( opts->prefix ) INFO_message("writing to '%s'...", opts->prefix);
      else                    INFO_message("writing to unknown?!?");
   }

   /* actually write results */
   for( c=0; c < opts->nt; c++ )
      fprintf(fp, "%f\n", opts->result[c]);

   if( fp != stdout && fp != stderr )
      fclose(fp);

   RETURN(0);
}

/*--------------- main routine ---------------*/
int main( int argc, char *argv[] )
{
   THD_3dim_dataset * countset=NULL;
   param_t          * params = &g_params;
   int                rv, limit;

   if( argc < 1 ) { show_help();  return 0; }

   /* general stuff */
   mainENTRY("3dmask_tool"); machdep(); AFNI_logger("3dmask_tool",argc,argv);
   enable_mcw_malloc();

   /* process options: a negative return is considered an error */
   rv = process_opts(params, argc, argv);
   if( rv ) RETURN(rv < 0);

   /* open, convert to byte, zeropad, dilate, unzeropad */
   if( process_input_dsets(params) ) RETURN(1);

   /* create mask count dataset and return num volumes (delete old dsets) */
   if( count_masks(params->dsets, params->ndsets, params->verb,
                   &countset, &params->nvols) ) RETURN(1);

   /* limit to frac of nvols (if not counting, convert to 0/1 mask) */
   limit = ceil((params->frac>1) ? params->frac : params->nvols*params->frac );
   if( params->verb )
      INFO_message("frac %g over %d volumes gives min count %d\n",
                   params->frac, params->nvols, limit);
   if( limit <= 0 ) limit = 1;

   /* if not counting, result is binary 0/1 */
   if( limit_to_frac(countset, limit, params->count, params->verb) )
      RETURN(1);

   /* maybe apply dilations to output */
   if( params->RESD.num > 0 ) {
      countset = apply_dilations(countset, &params->RESD, 0, params->verb);
      if( !countset ) RETURN(1);
   }

   /* maybe fill any remaining holes */
   if( params->fill )
      if ( fill_holes(countset, params->verb) ) RETURN(1);

   /* create output */
   if( write_result(params, countset, argc, argv) ) RETURN(1);

   /* clean up a little memory */
   DSET_delete(countset);
   free(params->dsets);

   RETURN(0);
}

/*
 * check count against limit
 *    - clear small values
 *    - if not count, set large values to 1
 */
int limit_to_frac(THD_3dim_dataset * cset, int limit, int count, int verb)
{
   short * dptr;
   int     index, nsub, nsuper;

   ENTRY("limit_to_frac");

   if( ! ISVALID_3DIM_DATASET(cset) ) {
      ERROR_message("invalid count dataset");
      RETURN(1);
   } else if( DSET_BRICK_TYPE(cset, 0) != MRI_short ) {
      ERROR_message("count dataset not of type short");
      RETURN(1);
   }

   if(verb > 1) INFO_message("limiting to %d (count = %d)\n",limit,count);

   /* note how many voxels are affected, just for kicks */
   dptr = DBLK_ARRAY(cset->dblk, 0);
   nsub = nsuper = 0;
   for(index = 0; index < DSET_NVOX(cset); index++, dptr++) {
      if( ! *dptr ) continue;           /* 0, so skip */
      else if( *dptr < limit ) {        /* small, so clear */
         *dptr = 0;
         nsub++;
      }
      else {                            /* big enough */
         if ( ! count ) *dptr = 1;
         nsuper++;
      }
   }

   /* entertain the user */
   if( verb )
      INFO_message("voxel limits: %d clipped, %d survived, %d were zero\n",
                   nsub, nsuper, DSET_NVOX(cset)-nsub-nsuper);

   RETURN(0);
}

static void toy_tsfunc( double tzero, double tdelta ,
                        int npts, float ts[],
                        double ts_mean , double ts_slope ,
                        void *ud, int nbriks, float *val          )
{
   static complex *comp_array=NULL;
   float mag=0, pow=0, pow6=0;
   int  jj;
   TOY_UD *rpud = (TOY_UD *)ud; 
   
   if( val == NULL ){
      INFO_message("toy_tsfunc: %s notification call, npts=%d\n", 
                   npts?"Start":"End", npts);
      if( npts > 0 ){  /* the "start notification" */
         /* This is when you perform any setup you don't want to repeat
            each time the function is called for a new voxel.
            Such a setup typically involves allocating for temporary
            arrays */
         
         /* allocate for fft array */
         comp_array = (complex *) calloc( sizeof(complex) , rpud->nfft);
      } else {  /* the "end notification" */
         /* Last call, meant for cleanup */
         if (comp_array) free(comp_array); comp_array = NULL;
      }
      return ;
   }
   
   /* if we get here then we have data to process */
   
   /* Load time series */
   for( jj=0 ; jj < npts ; jj++ ) {
      comp_array[jj].r = ts[jj]; comp_array[jj].i = 0.0f ;
   }
   /* zero pad */
   for( jj=npts ; jj < rpud->nfft ; jj++ )
       comp_array[jj].r = comp_array[jj].i = 0.0f ;  
   
   csfft_cox( -1 , rpud->nfft, comp_array ) ;   /* DFT */
    
   for( jj=0 ; jj < rpud->nfft ; jj++ ) {
      mag = CABS(comp_array[jj]) ;
      pow += mag*mag;
      if (jj<rpud->nfft/6) pow6 += mag*mag ; 
   }    
   
   val[0] =  pow;
   val[1] =  pow6;
   
   return ;
}

void AFNI_pbar_lock_carryout( Three_D_View *im3d )
{
   Three_D_View *qq3d ;
   static int busy = 0 ;  /* !=0 if this routine is "busy" */
   int glock , cc,ii ;

ENTRY("AFNI_pbar_lock_carryout") ;

   /* first, determine if there is anything to do */

   glock = GLOBAL_library.controller_lock ;     /* not a handgun */

   if( busy )                         EXRETURN;  /* routine already busy */
   if( glock == 0 )                   EXRETURN;  /* nothing to do */
   if( !IM3D_OPEN(im3d) )             EXRETURN;  /* bad input */
   if( GLOBAL_library.ignore_lock )   EXRETURN;  /* ordered not to do anything */
   if( !AFNI_check_pbar_lock() )      EXRETURN;  /* not locked? */

   ii = AFNI_controller_index(im3d) ;           /* which one am I? */

   if( ii < 0 ) EXRETURN ;                      /* nobody? bad input! */
   if( ((1<<ii) & glock) == 0 ) EXRETURN ;      /* input not locked */

#if 0
INFO_message("AFNI_pbar_lock_carryout( %d ) ******************* ",ii) ;
TRACEBACK ;
#endif

   /* something to do? */

   busy = 1 ;  /* don't let this routine be called recursively */

   /* loop through other controllers:
        for those that ARE open, ARE NOT the current one, and ARE locked */

   for( cc=0 ; cc < MAX_CONTROLLERS ; cc++ ){

     qq3d = GLOBAL_library.controllers[cc] ; /* controller */

     if( IM3D_OPEN(qq3d) && qq3d != im3d && ((1<<cc) & glock) != 0 ){
#if 0
ININFO_message(" equate_pbars( %d )",cc) ;
#endif
       AFNI_equate_pbars( qq3d , im3d ) ;
     }
   }

   busy = 0 ;  /* OK, let this routine be activated again */
   EXRETURN ;
}

/* convert by hand, since no scaling will be done
 * (byte seems inappropriate and float does not need it)  */
int write_result(param_t * params, THD_3dim_dataset * oset,
                 int argc, char * argv[])
{
   short * sptr;
   int     nvox = DSET_NVOX(oset), ind;

   ENTRY("write_results");

   EDIT_dset_items(oset, ADN_prefix, params->prefix, ADN_none);

   if( params->verb )
      INFO_message("writing result %s...\n", DSET_PREFIX(oset));

   switch( params->datum ) {
      default: ERROR_exit("invalid datum for result: %d", params->datum);
      case MRI_short: break;     /* nothing to do */
      case MRI_float: {
         float * data = (float *)malloc(nvox*sizeof(float));
         sptr = DBLK_ARRAY(oset->dblk, 0);
         if( ! data ) ERROR_exit("failed to alloc %d output floats\n", nvox);
         for( ind = 0; ind < nvox; ind++ ) data[ind] = (float)sptr[ind];
         EDIT_substitute_brick(oset, 0, params->datum, data);
      }
      break;
      case MRI_byte: {
         byte * data = (byte *)malloc(nvox*sizeof(byte));
         int errs = 0;
         sptr = DBLK_ARRAY(oset->dblk, 0);
         if( ! data ) ERROR_exit("failed to alloc %d output bytes\n", nvox);
         for( ind = 0; ind < nvox; ind++ ) {
            if( sptr[ind] > 255 ) {     /* watch for overflow */
               data[ind] = (byte)255;
               errs++;
            } else data[ind] = (byte)sptr[ind];
         }
         EDIT_substitute_brick(oset, 0, params->datum, data);
         if(errs) WARNING_message("convert to byte: %d truncated voxels",errs);
      }
      break;
   }

   tross_Make_History( "3dmask_tool", argc, argv, oset );

   DSET_write(oset);
   WROTE_DSET(oset);

   RETURN(0);
}

SEXP getListElement(SEXP list, const char *str)
{
   int debug=0;
   R_len_t i;
   SEXP elmt = R_NilValue, names = getAttrib(list, R_NamesSymbol);
   
   if (!debug) debug = get_odebug();
   
   for ( i = 0; i < length(list); i++) {
      if (debug) INFO_message("Element %d/%d: named %s\n", 
                     i, length(list), CHAR(STRING_ELT(names, i)));
      if(strcmp(CHAR(STRING_ELT(names, i)), str) == 0) {
         elmt = VECTOR_ELT(list, i);
         break;
      }
   }
   return elmt;
}

void AFNI_all_locks_carryout( Three_D_View *im3d )
{
   static int busy = 0 ;  /* !=0 if this routine is "busy" */

ENTRY("AFNI_all_locks_carryout") ;

   if( busy || !IM3D_VALID(im3d)  ) EXRETURN ;
   if( GLOBAL_library.ignore_lock ) EXRETURN ;

#if 0
INFO_message("AFNI_all_locks_carryout: im3d index = %d",AFNI_controller_index(im3d)) ;
#endif

   AFNI_lock_carryout       ( im3d ) ;
   AFNI_time_lock_carryout  ( im3d ) ;
   AFNI_thresh_lock_carryout( im3d ) ;
   AFNI_pbar_lock_carryout  ( im3d ) ;
   AFNI_range_lock_carryout ( im3d ) ;

   AFNI_sleep(1) ; busy = 0 ; EXRETURN ;
}

void EDIT_misfit_report( char *name, int ib,
                         int nxyz, float fac, short *sar, float *far )
{
    float mf ;
    int im ;
    static char *msg[5] = { "* Caution"  , "** Take Care"     ,
                            "*** Beware" , "**** Red Alert ****" ,
                            "***** Purple Alert! *****"
                          } ;
    static int first=1 ;

    mf = 100.0f * EDIT_scale_misfit( nxyz , fac , sar , far ) ;
    if( mf <=  9.0f ) return ;  /* OK */
    if( mf <= 13.0f ) im = 0 ;
    else if( mf <= 19.0f ) im = 1 ;
    else if( mf <= 27.0f ) im = 2 ;
    else if( mf <= 39.0f ) im = 3 ;
    else                   im = 4 ;
    if( first )
        WARNING_message("+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");

    WARNING_message("%s[%d] scale to shorts mean misfit error = %.1f%% -- %s",
                    name , ib , mf , msg[im] ) ;
    if( first ) {
        ININFO_message("a) Numerical precision has been lost when truncating results\n"
                       "       from 32-bit floating point to 16-bit integers (shorts)." );
        ININFO_message("b) Consider writing datasets out in float format.\n"
                       "       In most AFNI programs, use the '-float' option.");
        ININFO_message("c) This warning is a new message, but is an old issue\n"
                       "       that arises when storing results in an integer format."  );
        ININFO_message("d) Don't panic! These messages likely originate in peripheral\n"
                       "       or unimportant voxels. They mean that you must examine your output.\n"
                       "       \"Assess the situation and keep a calm head about you,\n"
                       "        because it doesn't do anybody any good to panic.\"\n" ) ;
        INFO_message("++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++");
        first = 0 ;
    }
    return ;
}

void RT_test_callback(void *junk)
{
   RT_status *rts = GLOBAL_library.realtime_status ;
   int cc , nval,nbr ;

   if( rts == NULL ){ ERROR_message("bad call to RT_test_callback"); return; }

   INFO_message("RT_test_callback: numchan=%d status=%d numdset=%d",
                rts->numchan , rts->status , rts->numdset ) ;

   for( cc=0 ; cc < rts->numdset ; cc++ ){     /* print out some dataset info */
     if( !ISVALID_DSET(rts->dset[cc]) ){
       ININFO_message(" dset[%d] invalid!",cc) ;       /* should never happen */
     } else {
       nval = DSET_NVALS(rts->dset[cc]) ;                 /* number of bricks */
       nbr  = THD_count_databricks(rts->dset[cc]->dblk) ; /* number with data */
       ININFO_message(" dset[%d] '%s': nvals=%d  nbr=%d",
                      cc , DSET_HEADNAME(rts->dset[cc]) , nval,nbr ) ;
     }
   }
   return ;
}

/*
 * create empty count dataset
 * for each input dataset and each sub-volume
 *    for each voxel, if set: increment
 * close datasets as they are processed
 */
int count_masks(THD_3dim_dataset * dsets[], int ndsets, int verb, /* inputs */
                THD_3dim_dataset ** cset, int * nvol)             /* outputs */
{
   THD_3dim_dataset * dset;
   short * counts = NULL;             /* will become data for returned cset */
   byte  * bptr;                      /* always points to mask volumes      */
   int     nxyz, iset, ivol, ixyz;

   ENTRY("count_masks");

   if( !dsets || !cset || !nvol )
      ERROR_exit("NULL inputs to count_masks");

   if( ndsets <= 0 ) {
      ERROR_message("count_masks: no input datasets");
      RETURN(1);
   }

   *nvol = 0;
   nxyz = DSET_NVOX(dsets[0]);
   
   /* allocate memory for the counts */
   counts = (short *)calloc(nxyz, sizeof(short));
   if( !counts ) ERROR_exit("failed to malloc %d shorts", nxyz);

   /* for each volume of each dataset, count set voxels */
   for( iset=0; iset < ndsets; iset++ ) {
      dset = dsets[iset];
      *nvol += DSET_NVALS(dset);        /* accumulate num volumes */

      /* for each volume in this dataset, count set voxels */
      for( ivol=0; ivol < DSET_NVALS(dset); ivol++ ) {
         if( DSET_BRICK_TYPE(dset, ivol) != MRI_byte )
            ERROR_exit("in count_masks with non-byte data (set %d, vol %d)",
                       iset, ivol);

         bptr = DBLK_ARRAY(dset->dblk, ivol);
         for( ixyz = 0; ixyz < nxyz; ixyz++ ) 
            if( bptr[ixyz] ) counts[ixyz]++;
      }

      if( iset > 0 ) DSET_delete(dset); /* close the first one at end */
   }  /* dataset */

   if( verb > 1 ) {
      int maxval;
      for( maxval=counts[0], ixyz=1; ixyz < nxyz; ixyz++ )
         if( counts[ixyz] > maxval ) maxval = counts[ixyz];

      INFO_message("counted %d mask volumes in %d datasets (%d voxels)\n",
                   *nvol, ndsets, nxyz);
      INFO_message("   (maximum overlap = %d)\n", maxval);
   }

   if( *nvol >= (1<<15) )
      WARNING_message("too many volumes to count as shorts: %d", *nvol);

   /* create output dataset */
   *cset = EDIT_empty_copy(dsets[0]);
   EDIT_dset_items(*cset, ADN_nvals, 1,  ADN_ntt, 0, ADN_none);
   EDIT_substitute_brick(*cset, 0, MRI_short, counts);

   DSET_delete(dsets[0]);  /* now finished with first dataset */

   RETURN(0);
}