自学教程：C++ DSET_NVOX函数代码示例

/*    Get copy contents of sub-brick iv into an double array.   if iv == -1, get the entire dset */double *THD_extract_to_double( int iv , THD_3dim_dataset *dset ){   MRI_IMAGE *im ;   double *var=NULL, *vv=NULL;   register int ii , nvox ;   ENTRY("THD_extract_to_double") ;   if (!dset) RETURN(var);   if (iv >= 0) {      if (!(im = THD_extract_double_brick(iv, dset))) RETURN(var);      var = MRI_DOUBLE_PTR(im);mri_fix_data_pointer(NULL, im);                               mri_free(im);im=NULL;   } else if (iv == -1) {      if (!(var = (double *)calloc(DSET_NVOX(dset)*DSET_NVALS(dset),                                  sizeof(double)))){         ERROR_message("Failed to allocate");         RETURN(NULL);      }      for (ii=0; ii<DSET_NVALS(dset); ++ii) {         if (!(im = THD_extract_double_brick(ii, dset))) {            ERROR_message("Failed toextract sb %d from dset", ii);            if (var) free(var);            RETURN(NULL);         }         vv = MRI_DOUBLE_PTR(im);         memcpy(var+ii*DSET_NVOX(dset),vv, sizeof(double)*DSET_NVOX(dset));          mri_free(im);im=NULL;      }   } else {      ERROR_message("Bad value of %d/n", iv);   }      RETURN(var);}

/* * 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);}

int SUMA_ShortizeDset(THD_3dim_dataset **dsetp, float thisfac) {   static char FuncName[]={"SUMA_ShortizeDset"};   char sprefix[THD_MAX_PREFIX+10];   int i, j;   byte *bb=NULL;   short *sb=NULL;   float bbf=0.0;      THD_3dim_dataset *cpset=NULL, *dset=*dsetp;      SUMA_ENTRY;      if (!dset) {      SUMA_S_Err("NULL *dsetp at input!");      SUMA_RETURN(0);   }      sprintf(sprefix, "%s.s", dset->dblk->diskptr->prefix);   NEW_SHORTY(dset, DSET_NVALS(dset), "ss.cp", cpset);         for (i=0; i<DSET_NVALS(dset); ++i) {      if (DSET_BRICK_TYPE(dset,i) == MRI_byte) {         bb = (byte *)DSET_ARRAY(dset,i);         sb = (short *)DSET_ARRAY(cpset,i);         if (thisfac <= 0.0) {            for (j=0; j<DSET_NVOX(dset); ++j) {               sb[j] = (short)bb[j];            }            thisfac = DSET_BRICK_FACTOR(dset,i);         } else {            bbf = DSET_BRICK_FACTOR(dset,i); if (bbf == 0.0f) bbf = 1.0;            bbf = bbf/thisfac;            for (j=0; j<DSET_NVOX(dset); ++j) {               sb[j] = SHORTIZE((((float)bb[j])*bbf));            }         }         EDIT_BRICK_FACTOR( cpset,i,thisfac ) ;      } else {         EDIT_substscale_brick(cpset, i, DSET_BRICK_TYPE(dset,i),                             DSET_ARRAY(dset,i), MRI_short, thisfac);         if (DSET_BRICK_TYPE(dset,i) != MRI_short) {            DSET_FREE_ARRAY(dset, i);         } else {            DSET_NULL_ARRAY(dset, i);         }      }   }   /* preserve tables, if any */   THD_copy_labeltable_atr( cpset->dblk,  dset->dblk);    DSET_delete(dset); dset = NULL;    *dsetp=cpset;   SUMA_RETURN(1);}

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);}

/*---------------------------------------------------------------------  23 Feb 2012: Return the absolute value of the difference between                two volumes, divided by the number of voxels               and the number of sub-bricks. Voxels that are zero               in both sets are not counted.               Comparisons are done after conversion of data to double    return = -1.0 ERROR           =  0.0 Exactly the same-----------------------------------------------------------------------*/double THD_diff_vol_vals(THD_3dim_dataset *d1, THD_3dim_dataset *d2, int scl) {   double dd=0.0, denom=0.0;   int i=0, k=0;   double *a1=NULL, *a2=NULL;   MRI_IMAGE *b1 = NULL , *b2 = NULL;      ENTRY("THD_diff_vol_vals");      if (!d1 && !d2) RETURN(dd);   if (!d1 || !d2) RETURN(-1.0);   if (!EQUIV_GRIDS(d1,d2)) RETURN(-1.0);   if (DSET_NVALS(d1) != DSET_NVALS(d2)) RETURN(-1.0);     DSET_mallocize(d1) ; DSET_load(d1) ;   DSET_mallocize(d2) ; DSET_load(d2) ;   dd = 0.0; denom = 0;   for (i=0; i<DSET_NVALS(d1); ++i) {      b1 = THD_extract_double_brick(i, d1);      b2 = THD_extract_double_brick(i, d2);      a1 = MRI_DOUBLE_PTR(b1);      a2 = MRI_DOUBLE_PTR(b2);      for (k=0; k<DSET_NVOX(d1); ++k) {         dd += ABS(a1[k]-a2[k]);         if (a1[k]!=0.0 || a2[k]!=0.0) ++denom;      }      mri_clear_data_pointer(b1); mri_free(b1) ;      mri_clear_data_pointer(b2); mri_free(b2) ;   }   if (scl && denom>0.0) dd /= denom;      RETURN(dd);   }  

/* * 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);}

MRI_IMAGE * THD_median_brick( THD_3dim_dataset *dset ){   int nvox , nvals , ii ;   MRI_IMAGE *tsim , *medim ;   float *medar ;   float *tsar ;  /* 05 Nov 2001 */ENTRY("THD_median_brick") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ;   if( !DSET_LOADED(dset) ) RETURN(NULL) ;   nvals = DSET_NVALS(dset) ;   tsim  = DSET_BRICK(dset,0) ;   if( nvals == 1 ){     medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;     RETURN(medim) ;   }   medim = mri_new_conforming( tsim , MRI_float ) ;   medar = MRI_FLOAT_PTR(medim) ;   nvox  = DSET_NVOX(dset) ;   tsar = (float *) calloc( sizeof(float),nvals+1 ) ; /* 05 Nov 2001 */   for( ii=0 ; ii < nvox ; ii++ ){     THD_extract_array( ii , dset , 0 , tsar ) ;     /* 05 Nov 2001 */     medar[ii] = qmed_float( nvals , tsar ) ;   }   free(tsar) ; RETURN(medim) ;}

/* just make sure we have sufficient data for computations */int check_dims(options_t * opts){   int nt, nvox, nmask;   ENTRY("check_dims");   nt = DSET_NVALS(opts->inset);   nvox = DSET_NVOX(opts->inset);   if( opts->mask ) nmask = THD_countmask( nvox, opts->mask );   else             nmask = nvox;   /* make sure we have something to compute */   if( nvox < 1 ) {      ERROR_message("input dataset must have at least 1 voxel");      RETURN(1);   } else if( nmask < 1 ) {      ERROR_message("input mask must have at least 1 voxel");      RETURN(1);   } else if( nt < 2 ) {      ERROR_message("input dataset must have at least 2 time points");      RETURN(1);   }   RETURN(0);}

MRI_IMARR * THD_medmad_bricks( THD_3dim_dataset *dset ){   int nvox , nvals , ii ;   MRI_IMAGE *tsim , *madim, *medim ;   float             *madar, *medar ;   MRI_IMARR *imar ;   float *tsar ;ENTRY("THD_medmad_bricks") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   nvals = DSET_NVALS(dset) ; if( nvals == 1 ) RETURN(NULL) ;   DSET_load(dset) ;  if( !DSET_LOADED(dset) ) RETURN(NULL) ;   tsim  = DSET_BRICK(dset,0) ;   madim = mri_new_conforming( tsim , MRI_float ) ;   madar = MRI_FLOAT_PTR(madim) ;   medim = mri_new_conforming( tsim , MRI_float ) ;   medar = MRI_FLOAT_PTR(medim) ;   nvox  = DSET_NVOX(dset) ;   tsar = (float *) calloc( sizeof(float),nvals+1 ) ;   for( ii=0 ; ii < nvox ; ii++ ){     THD_extract_array( ii , dset , 0 , tsar ) ;     qmedmad_float( nvals , tsar , medar+ii , madar+ii ) ;   }   free(tsar) ;   INIT_IMARR(imar) ; ADDTO_IMARR(imar,medim) ; ADDTO_IMARR(imar,madim) ;   RETURN(imar) ;}

MRI_IMAGE * THD_rms_brick( THD_3dim_dataset *dset ){   int nvox , nvals , ii , jj ;   MRI_IMAGE *tsim , *medim ;   float *medar , sum,fac ;   float *tsar ;ENTRY("THD_rms_brick") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ;   if( !DSET_LOADED(dset) ) RETURN(NULL) ;   nvals = DSET_NVALS(dset)   ; fac = 1.0 / nvals ;   tsim  = DSET_BRICK(dset,0) ;   if( nvals == 1 ){     medim = mri_scale_to_float( DSET_BRICK_FACTOR(dset,0), tsim ) ;     RETURN(medim) ;   }   medim = mri_new_conforming( tsim , MRI_float ) ;   medar = MRI_FLOAT_PTR(medim) ;   nvox  = DSET_NVOX(dset) ;   tsar = (float *) calloc( sizeof(float),nvals+1 ) ;   for( ii=0 ; ii < nvox ; ii++ ){     THD_extract_array( ii , dset , 0 , tsar ) ;     for( sum=0.0,jj=0 ; jj < nvals ; jj++ ) sum += tsar[jj]*tsar[jj] ;     medar[ii] = sqrtf(fac * sum) ;   }   free(tsar) ; RETURN(medim) ;}

int64_t THD_vectim_size( THD_3dim_dataset *dset , byte *mask ){   int nvals , nvox , nmask ;   int64_t sz ;   ENTRY("THD_vectim_size") ;   if( !ISVALID_DSET(dset) ) RETURN(0) ;   nvals = DSET_NVALS(dset) ;   nvox  = DSET_NVOX(dset) ;   if( mask != NULL ) nmask = THD_countmask( nvox , mask ) ;   else               nmask = DSET_NVOX(dset) ;   sz = ((int64_t)nmask) * ( ((int64_t)nvals) * sizeof(float) + sizeof(int) ) ;   RETURN(sz) ;}

bytevec * THD_create_mask_from_string( char *str )  /* Jul 2010 */{   bytevec *bvec=NULL ; int nstr ; char *buf=NULL ;ENTRY("THD_create_mask") ;   if( str == NULL || *str == '/0' ) RETURN(NULL) ;   nstr = strlen(str) ;   bvec = (bytevec *)malloc(sizeof(bytevec)) ;   /* try to read it as a dataset */   if( nstr < THD_MAX_NAME ){     THD_3dim_dataset *dset = THD_open_one_dataset(str) ;     if( dset != NULL ){       bvec->nar = DSET_NVOX(dset) ;       bvec->ar  = THD_makemask( dset , 0 , 1.0f,0.0f ) ;       DSET_delete(dset) ;       if( bvec->ar == NULL ){         ERROR_message("Can't make mask from dataset '%s'",str) ;         free(bvec) ; bvec = NULL ;       }       RETURN(bvec) ;     }   }   /* if str is a filename, read that file;      otherwise, use the string itself to find the mask */   if( THD_is_file(str) ){     buf = AFNI_suck_file(str) ;     if( buf != NULL ) nstr = strlen(buf) ;   } else {     buf = str ;   }   /* try to read buf as a Base64 mask string */   if( strrchr(buf,'=') != NULL ){     int nvox ;     bvec->ar = mask_from_b64string( buf , &nvox ) ;     if( bvec->ar != NULL ){       bvec->nar = nvox ;     } else {       ERROR_message("Can't make mask from string '%.16s' %s",buf,(nstr<=16)?" ":"...") ;       free(bvec) ; bvec = NULL ;     }   } else {     ERROR_message("Don't understand mask string '%.16s'",buf,(nstr<=16)?" ":"...") ;     free(bvec) ; bvec = NULL ;   }   if( buf != str && buf != NULL ) free(buf) ;   RETURN(bvec) ;}

/* * 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);}

int THD_voxel_is_constant( int ind , THD_3dim_dataset *dset ){   float *far ; int ii,nvox,nvals ;   if( !ISVALID_DSET(dset) ) return 1 ;   if( ind < 0 || ind >= DSET_NVOX(dset) ) return 1 ;   nvals = DSET_NVALS(dset) ; if( nvals == 1 ) return 1 ;   far = (float *)malloc(sizeof(float)*nvals) ; NULL_CHECK(far) ;   ii = THD_extract_array( ind , dset , 0 , far ) ;   if( ii < 0 ){ free(far); return 1; }   for( ii=1 ; ii < nvals && far[ii]==far[0]; ii++ ) ; /*nada*/   free(far) ; return (ii==nvals) ;}

float THD_get_float_value( int ind , int ival , THD_3dim_dataset *dset ){   MRI_TYPE typ ; float val=0.0f ;   if( ind < 0 || ival < 0 || !ISVALID_DSET(dset) ||       ival >= DSET_NVALS(dset) || ind >= DSET_NVOX(dset) ) return val ;   typ = DSET_BRICK_TYPE(dset,ival) ;  /* raw data type */   switch( typ ){      default:           /* don't know what to do --> return nada */         return(-1);      break ;      case MRI_byte:{         byte *bar ;         bar = (byte *) DSET_ARRAY(dset,ival) ;         if( bar != NULL ) val = (float)bar[ind] ;      }      break ;      case MRI_short:{         short *bar ;         bar = (short *) DSET_ARRAY(dset,ival) ;         if( bar != NULL ) val = (float)bar[ind] ;      }      break ;      case MRI_float:{         float *bar ;         bar = (float *) DSET_ARRAY(dset,ival) ;         if( bar != NULL ) val = bar[ind] ;      }      break ;      case MRI_complex:{         complex *bar ;         bar = (complex *) DSET_ARRAY(dset,ival) ;         if( bar != NULL ) val = CABS(bar[ind]) ;      }      break ;   }   if( DSET_BRICK_FACTOR(dset,ival) > 0.0f )     val *= DSET_BRICK_FACTOR(dset,ival) ;   return val ;}

/* 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);}

int is_integral_sub_brick ( THD_3dim_dataset *dset, int isb, int check_values){   float mfac = 0.0;   void *vv=NULL;   if(   !ISVALID_DSET(dset)    ||            isb < 0                     ||            isb >= DSET_NVALS(dset)  ) {      fprintf(stderr,"** Bad dset or sub-brick index./n");      return (0) ;   }   if( !DSET_LOADED(dset) ) DSET_load(dset);   switch( DSET_BRICK_TYPE(dset,isb) ){      case MRI_short:      case MRI_byte:         if (check_values) {            mfac = DSET_BRICK_FACTOR(dset,isb) ;            if (mfac != 0.0f && mfac != 1.0f) return(0);         }         break;      case MRI_double:      case MRI_complex:      case MRI_float:         vv = (void *)DSET_ARRAY(dset,isb);         mfac = DSET_BRICK_FACTOR(dset,isb) ;         if (mfac != 0.0f && mfac != 1.0f) return(0);         if (!vv) {            fprintf(stderr,"** NULL array!/n");            return(0);         }         return(is_integral_data(DSET_NVOX(dset),                                 DSET_BRICK_TYPE(dset,isb),                                 DSET_ARRAY(dset,isb) ) );         break;      default:         return(0);   }   return(1);}

MRI_vectim * THD_dset_to_vectim_byslice( THD_3dim_dataset *dset, byte *mask ,                                         int ignore , int kzbot , int kztop  ){   byte *mmm ;   MRI_vectim *mrv=NULL ;   int kk,iv , nvals , nvox , nmask , nxy , nz ;ENTRY("THD_dset_to_vectim_byslice") ;                     if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ; if( !DSET_LOADED(dset)  ) RETURN(NULL) ;   nvals = DSET_NVALS(dset) ; if( nvals <= 0 ) RETURN(NULL) ;   nvox  = DSET_NVOX(dset) ;   nxy = DSET_NX(dset) * DSET_NY(dset) ; nz = DSET_NZ(dset) ;   if( kzbot <  0  ) kzbot = 0 ;   if( kztop >= nz ) kztop = nz-1 ;   if( kztop < kzbot ) RETURN(NULL) ;   if( kzbot == 0 && kztop == nz-1 ){     mrv = THD_dset_to_vectim( dset , mask, ignore ) ; RETURN(mrv) ;   }   /* make a mask that includes cutting out un-desirable slices */   { int ibot , itop , ii ;#pragma omp critical (MALLOC)     mmm = (byte *)malloc(sizeof(byte)*nvox) ;     if( mask == NULL ) AAmemset( mmm ,    1 , sizeof(byte)*nvox ) ;     else               AAmemcpy( mmm , mask , sizeof(byte)*nvox ) ;     if( kzbot > 0 )       AAmemset( mmm               , 0 , sizeof(byte)*kzbot       *nxy ) ;     if( kztop < nz-1 )       AAmemset( mmm+(kztop+1)*nxy , 0 , sizeof(byte)*(nz-1-kztop)*nxy ) ;   }   /* and make the vectim using the standard function */   mrv = THD_dset_to_vectim( dset , mmm , ignore ) ;   free(mmm) ;   RETURN(mrv) ;}

MRI_IMAGE * THD_dset_to_1Dmri( THD_3dim_dataset *dset ){   MRI_IMAGE *im ; float *far ;   int nx , ny , ii ;ENTRY("THD_dset_to_1D") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   DSET_load(dset) ;   if( !DSET_LOADED(dset) ) RETURN(NULL) ;   nx = DSET_NVALS(dset) ;   ny = DSET_NVOX(dset) ;   im = mri_new( nx , ny , MRI_float ) ; far = MRI_FLOAT_PTR(im) ;   for( ii=0 ; ii < ny ; ii++ )     THD_extract_array( ii , dset , 0 , far + ii*nx ) ;   RETURN(im) ;}

int main( int argc , char * argv[] ){   int kk , nvox , ii ;   THD_3dim_dataset * oset ;   float * far ;   /*-- read command line arguments --*/   if( argc < 2 || strncmp(argv[1],"-help",5) == 0 ) UC_syntax(NULL) ;   (void) my_getenv("junk") ;   UC_read_opts( argc , argv ) ;   set_unusuality_tail( UC_ptail ) ;   oset = EDIT_empty_copy( UC_dset ) ;   EDIT_dset_items( oset ,                       ADN_prefix      , UC_prefix ,                       ADN_ntt         , 0 ,                       ADN_nvals       , 1 ,                       ADN_datum_all   , MRI_float ,                       ADN_malloc_type , DATABLOCK_MEM_MALLOC ,                    ADN_none ) ;   nvox = DSET_NVOX(oset) ;   far = (float *) malloc( sizeof(float) * nvox ) ;   for( kk=0 ; kk < nvox ; kk++ ) far[kk] = 0.0 ;   EDIT_substitute_brick( oset , 0 , MRI_float , far ) ;   if( !UC_be_quiet ){ printf("--- computing u") ; fflush(stdout) ; }   for( kk=0 ; kk < UC_nvec ; kk++ ){      ii = (UC_iv == NULL) ? kk : UC_iv[kk] ;      far[ii] = UC_unusuality( UC_vdim, UC_vec[kk] , UC_nvec, UC_vec ) ;      if( !UC_be_quiet && kk%1000==999 ){ printf(".");fflush(stdout); }   }   if( !UC_be_quiet ) printf("/n--- writing output/n") ;   DSET_write(oset) ;   exit(0) ;}

void THD_extract_detrended_array( THD_3dim_dataset *dset ,                                  int nref, float **ref, MRI_IMARR *imar,                                  int ii, int scl, float *far ){   int tt , nval , qq ;   float val , **fitar , *var ;   MRI_IMAGE *qim ;ENTRY("THD_extract_detrended_array") ;   if( !ISVALID_DSET(dset) ||       nref < 1            || ref == NULL                ||       imar == NULL        || IMARR_COUNT(imar) < nref+1 ||       ii < 0              || ii >= DSET_NVOX(dset)      || far == NULL ) EXRETURN ;   qq = THD_extract_array( ii , dset , 0 , far ) ;  /* get data */   if( qq < 0 ) EXRETURN ;   nval = DSET_NVALS(dset) ;   fitar = (float **)malloc(sizeof(float *)*nref) ;   for( qq=0 ; qq < nref ; qq++ ){     qim = IMARR_SUBIM(imar,qq) ; fitar[qq] = MRI_FLOAT_PTR(qim) ;   }   qim = IMARR_SUBIM(imar,nref) ; var = MRI_FLOAT_PTR(qim) ;   for( tt=0 ; tt < nval ; tt++ ){  /* get residuals */     val = far[tt] ;     for( qq=0 ; qq < nref ; qq++ ) val -= ref[qq][tt] * fitar[qq][ii] ;     far[tt] = val ;   }   if( scl && var[ii] > 0.0f ){     val = 1.0f / var[ii] ;     for( tt=0 ; tt < nval ; tt++ ) far[tt] *= val ;   }      /* ZSS: Need to free fitar */   free(fitar); fitar=NULL;   EXRETURN ;}

THD_3dim_dataset * THD_detrend_dataset( THD_3dim_dataset *dset ,                                        int nref , float **ref ,                                        int meth , int scl ,                                        byte *mask , MRI_IMARR **imar ){   MRI_IMARR *qmar ;   int ii,jj,kk , nvals,nvox , iv ;   float *var ;   THD_3dim_dataset *newset ;ENTRY("THD_detrend_dataset") ;   if( !ISVALID_DSET(dset) ) RETURN(NULL) ;   nvals = DSET_NVALS(dset) ; nvox = DSET_NVOX(dset) ;   qmar = THD_time_fit_dataset( dset , nref,ref , meth , mask ) ;   if( qmar == NULL ) RETURN(NULL) ;   newset = EDIT_empty_copy(dset) ;   for( iv=0 ; iv < nvals ; iv++ ){     EDIT_substitute_brick( newset , iv , MRI_float , NULL ) ;     EDIT_BRICK_FACTOR( newset , iv , 0.0f ) ;  /* 04 Jun 2007 */   }   var = (float *)malloc(sizeof(float)*nvals) ;   for( ii=0 ; ii < nvox ; ii++ ){     if( mask == NULL || mask[ii] )       THD_extract_detrended_array( dset , nref,ref , qmar , ii,scl , var ) ;     else       memset(var,0,sizeof(float)*nvals) ;     THD_insert_series( ii , newset , nvals , MRI_float , var , 0 ) ;   }   free(var) ;   if( imar != NULL ) *imar = qmar ;   else               DESTROY_IMARR(qmar) ;   RETURN(newset) ;}

int THD_retrend_dataset( THD_3dim_dataset *dset ,                         int nref , float **ref ,                         int scl , byte *mask , MRI_IMARR *imar ){   int ii,qq,tt , nvals,nvox ;   MRI_IMAGE *qim ; float **fitar , *far , fac , *var , val ;ENTRY("THD_retrend_dataset") ;   if( !ISVALID_DSET(dset) ||       nref < 1            || ref == NULL ||       imar == NULL        || IMARR_COUNT(imar) <= nref ) RETURN(0) ;   nvals = DSET_NVALS(dset) ; nvox = DSET_NVOX(dset) ;   fitar = (float **)malloc(sizeof(float *)*nref) ;   for( qq=0 ; qq < nref ; qq++ ){     qim = IMARR_SUBIM(imar,qq) ; fitar[qq] = MRI_FLOAT_PTR(qim) ;   }   qim = IMARR_SUBIM(imar,nref) ; var = MRI_FLOAT_PTR(qim) ;   far  = (float *)malloc(sizeof(float)*nvals) ;   for( ii=0 ; ii < nvox ; ii++ ){     if( mask != NULL && !mask[ii] ) continue ;     qq = THD_extract_array( ii , dset , 0 , far ) ;  /* get data */     if( qq < 0 ) continue ;     fac = (scl) ? var[ii] : 1.0f ;     for( tt=0 ; tt < nvals ; tt++ ){  /* add fit back in */       val = far[tt] * fac ;       for( qq=0 ; qq < nref ; qq++ ) val += ref[qq][tt] * fitar[qq][ii] ;       far[tt] = val ;     }     THD_insert_series( ii , dset , nvals , MRI_float , far , 0 ) ;   }   free(far) ; free(fitar) ; RETURN(1) ;}

//.........这里部分代码省略.........		kk = THD_bandpass_set_nfft(nfft) ;		if( kk != nfft && verb )			INFO_message("Data length = %d  FFT length = %d",ntime,kk) ;   }   if( dt <= 0.0f ){		dt = DSET_TR(inset) ;		if( dt <= 0.0f ){			WARNING_message("Setting dt=1.0 since input dataset lacks a time axis!") ;			dt = 1.0f ;		}   }   ftopALL = 1./dt ;// Aug,2016: should solve problem of a too-large                    // value for THD_bandpass_vectors(), while still                    // being >f_{Nyquist}   if( !THD_bandpass_OK(ntime,dt,fbot,ftop,1) ) ERROR_exit("Can't continue!") ;   nx = DSET_NX(inset); ny = DSET_NY(inset); nz = DSET_NZ(inset); nvox = nx*ny*nz;   /* check mask, or create it */   if( verb ) INFO_message("Loading input dataset time series" ) ;   DSET_load(inset) ;   if( mask != NULL ){		if( mask_nx != nx || mask_ny != ny || mask_nz != nz )			ERROR_exit("-mask dataset grid doesn't match input dataset") ;   } else if( do_automask ){		mask = THD_automask( inset ) ;		if( mask == NULL )			ERROR_message("Can't create -automask from input dataset?") ;		nmask = THD_countmask( DSET_NVOX(inset) , mask ) ;		if( verb ) INFO_message("Number of voxels in automask = %d",nmask);		if( nmask < 1 ) ERROR_exit("Automask is too small to process") ;   } else {		mask = (byte *)malloc(sizeof(byte)*nvox) ; nmask = nvox ;		memset(mask,1,sizeof(byte)*nvox) ;		// if( verb ) // @@ alert if aaaalllllll vox are going to be analyzed!		INFO_message("No mask ==> processing all %d voxels",nvox);   }   /* A simple check of dataset quality [08 Feb 2010] */   if( !nosat ){		float val ;		INFO_message(						 "Checking dataset for initial transients [use '-notrans' to skip this test]") ;		val = THD_saturation_check(inset,mask,0,0) ; kk = (int)(val+0.54321f) ;		if( kk > 0 )			ININFO_message(								"Looks like there %s %d non-steady-state initial time point%s :-(" ,								((kk==1) ? "is" : "are") , kk , ((kk==1) ? " " : "s") ) ;		else if( val > 0.3210f )  /* don't ask where this threshold comes from! */			ININFO_message(								"MAYBE there's an initial positive transient of 1 point, but it's hard to tell/n") ;		else			ININFO_message("No widespread initial positive transient detected :-)") ;   }   /* check -dsort inputs for match to inset */   for( kk=0 ; kk < nortset ; kk++ ){		if( DSET_NX(ortset[kk])    != nx ||

//.........这里部分代码省略.........         }         polort = val ; nopt++ ; continue ;      }      if( strcmp(argv[nopt],"-mem_stat") == 0 ){         MEM_STAT = 1 ; nopt++ ; continue ;      }      if( strncmp(argv[nopt],"-mem_profile",8) == 0 ){         MEM_PROF = 1 ; nopt++ ; continue ;      }      /* check for 1d argument */      if ( strcmp(argv[nopt],"-out1D") == 0 ){          if (!(fout1D = fopen(argv[++nopt], "w"))) {             ERROR_message("Failed to open %s for writing", argv[nopt]);             exit(1);          }          nopt++ ; continue ;      }      ERROR_exit("Illegal option: %s",argv[nopt]) ;   }   /*-- open dataset, check for legality --*/   if( nopt >= argc ) ERROR_exit("Need a dataset on command line!?") ;   xset = THD_open_dataset(argv[nopt]); CHECK_OPEN_ERROR(xset,argv[nopt]);   if( DSET_NVALS(xset) < 3 )     ERROR_exit("Input dataset %s does not have 3 or more sub-bricks!",argv[nopt]) ;   DSET_load(xset) ; CHECK_LOAD_ERROR(xset) ;   /*-- compute mask array, if desired --*/   nvox = DSET_NVOX(xset) ; nvals = DSET_NVALS(xset) ;   INC_MEM_STATS((nvox * nvals * sizeof(double)), "input dset");   PRINT_MEM_STATS("inset");   /* if a mask was specified make sure it is appropriate */   if( mset ){      if( DSET_NVOX(mset) != nvox )         ERROR_exit("Input and mask dataset differ in number of voxels!") ;      mask  = THD_makemask(mset, 0, 1.0, 0.0) ;      /* update running memory statistics to reflect loading the image */      INC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );      PRINT_MEM_STATS( "mset load" );      nmask = THD_countmask( nvox , mask ) ;      INC_MEM_STATS( nmask * sizeof(byte), "mask array" );      PRINT_MEM_STATS( "mask" );      INFO_message("%d voxels in -mask dataset",nmask) ;      if( nmask < 2 ) ERROR_exit("Only %d voxels in -mask, exiting...",nmask);      /* update running memory statistics to reflect loading the image */      DEC_MEM_STATS( mset->dblk->total_bytes, "mask dset" );      DSET_unload(mset) ;      PRINT_MEM_STATS( "mset unload" );   }    /* if automasking is requested, handle that now */   else if( do_autoclip ){      mask  = THD_automask( xset ) ;      nmask = THD_countmask( nvox , mask ) ;      INFO_message("%d voxels survive -autoclip",nmask) ;      if( nmask < 2 ) ERROR_exit("Only %d voxels in -automask!",nmask);

/* * 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);}

THD_3dim_dataset * THD_localhistog( int nsar , THD_3dim_dataset **insar ,                                    int numval , int *rlist , MCW_cluster *nbhd ,                                    int do_prob , int verb ){   THD_3dim_dataset *outset=NULL , *inset ;   int nvox=DSET_NVOX(insar[0]) ;   int ids, iv, bb, nnpt=nbhd->num_pt ;   MRI_IMAGE *bbim ; int btyp ;   float **outar , **listar ;ENTRY("THD_localhistog") ;   /*---- create output dataset ----*/   outset = EDIT_empty_copy(insar[0]) ;   EDIT_dset_items( outset ,                      ADN_nvals     , numval    ,                      ADN_datum_all , MRI_float ,                      ADN_nsl       , 0         ,                      ADN_brick_fac , NULL      ,                    ADN_none ) ;   outar = (float **)malloc(sizeof(float *)*numval) ;   for( bb=0 ; bb < numval ; bb++ ){     EDIT_substitute_brick( outset , bb , MRI_float , NULL ) ;     outar[bb] = DSET_BRICK_ARRAY(outset,bb) ;   }   /*---- make mapping between values and arrays to get those values ----*/   listar = (float **)malloc(sizeof(float *)*TWO16) ;   for( bb=0 ; bb < TWO16 ; bb++ ) listar[bb] = outar[0] ;   for( bb=1 ; bb < numval ; bb++ ){     listar[ rlist[bb] + TWO15 ] = outar[bb] ;   }   /*----------- loop over datasets, add in counts for all voxels -----------*/   for( ids=0 ; ids < nsar ; ids++ ){              /* dataset loop */     inset = insar[ids] ; DSET_load(inset) ;     for( iv=0 ; iv < DSET_NVALS(inset) ; iv++ ){  /* sub-brick loop */       if( verb ) fprintf(stderr,".") ;       bbim = DSET_BRICK(inset,iv) ; btyp = bbim->kind ;       if( nnpt == 1 ){                            /* only 1 voxel in nbhd */         int qq,ii,jj,kk,ib,nb ;         switch( bbim->kind ){           case MRI_short:{             short *sar = MRI_SHORT_PTR(bbim) ;             for( qq=0 ; qq < nvox ; qq++ ) listar[sar[qq]+TWO15][qq]++ ;           }           break ;           case MRI_byte:{             byte *bar = MRI_BYTE_PTR(bbim) ;             for( qq=0 ; qq < nvox ; qq++ ) listar[bar[qq]+TWO15][qq]++ ;           }           break ;           case MRI_float:{             float *far = MRI_FLOAT_PTR(bbim) ; short ss ;             for( qq=0 ; qq < nvox ; qq++ ){ ss = SHORTIZE(far[qq]); listar[ss+TWO15][qq]++; }           }           break ;         }       } else {                                    /* multiple voxels in nbhd */ AFNI_OMP_START ;#pragma omp parallel { int qq,ii,jj,kk,ib,nb ; void *nar ; short *sar,ss ; byte *bar ; float *far ;   nar = malloc(sizeof(float)*nnpt) ;   sar = (short *)nar ; bar = (byte *)nar ; far = (float *)nar ;#pragma omp for         for( qq=0 ; qq < nvox ; qq++ ){           /* qq=voxel index */           ii = DSET_index_to_ix(inset,qq) ;           jj = DSET_index_to_jy(inset,qq) ;           kk = DSET_index_to_kz(inset,qq) ;           nb = mri_get_nbhd_array( bbim , NULL , ii,jj,kk , nbhd , nar ) ;           if( nb == 0 ) continue ;           switch( btyp ){             case MRI_short:               for( ib=0 ; ib < nb ; ib++ ) listar[sar[ib]+TWO15][qq]++ ;             break ;             case MRI_byte:               for( ib=0 ; ib < nb ; ib++ ) listar[bar[ib]+TWO15][qq]++ ;             break ;             case MRI_float:               for( ib=0 ; ib < nb ; ib++ ){ ss = SHORTIZE(far[ib]); listar[ss+TWO15][qq]++; }             break ;           }         } /* end of voxel loop */   free(nar) ; } /* end of OpenMP */ AFNI_OMP_END ;       }     } /* end of sub-brick loop */     DSET_unload(inset) ;   } /* end of dataset loop */   if( verb ) fprintf(stderr,"/n") ;   free(listar) ;   /*---- post-process output ---*///.........这里部分代码省略.........

//.........这里部分代码省略.........   }   /* number of bytes needed:        sizeof(datum) * number of vectors per dataset                      * number of datasets                      * sum of per dataset vector lengths */   nbytes_needed = 0 ;   for( ids=0 ; ids < ndset ; ids++ ) nbytes_needed += nvals[ids] ;   nbytes_needed *= ((long long)nvec) * datum_size ;   if( nbytes_needed >= twogig &&       ( sizeof(void *) < 8 || sizeof(size_t) < 8 ) ) /* too much for 32-bit */     GQUIT("datafile size exceeds 2 GB -- you need a 64-bit computer!") ;   /* scale factor for each dataset */   atr = NI_get_attribute(nel,"fac") ;   if( atr == NULL ) GQUIT("fac attribute missing") ;   facar = NI_decode_float_list(atr,",") ;   if( facar == NULL || facar->num < ndset )     GQUIT("can't decode fac attribute") ;   fac = facar->ar ; facar->ar = NULL ; NI_delete_float_array(facar) ;   for( ids=0 ; ids < ndset ; ids++ ) if( fac[ids] <= 0.0f ) fac[ids] = 1.0f ;   /* grid definition */   atr = NI_get_attribute(nel,"geometry") ;   if( atr == NULL ) GQUIT("geometry attribute missing") ;   geometry_string = strdup(atr) ;   tdset = EDIT_geometry_constructor( geometry_string , "GrpInCorr" ) ;   if( tdset == NULL ) GQUIT("can't decode geometry attribute") ;   nvox = DSET_NVOX(tdset) ;   if(  no_ivec && nvox != nvec )     GQUIT("geometry attribute doesn't match nvec attribute") ;   if( !no_ivec && nvox <  nvec )     GQUIT("geometry attribute specifies too few voxels") ;   /* name of data file: check its size against what's needed */#if 0   atr = NI_get_attribute(nel,"datafile") ;   if( atr != NULL ){     dfname = strdup(atr) ; nbytes_dfname = THD_filesize(dfname) ;     if( nbytes_dfname <= 0 && strstr(dfname,"/") != NULL ){       char *tnam = THD_trailname(atr,0) ;       nbytes_dfname = THD_filesize(tnam) ;       if( nbytes_dfname > 0 ){ free(dfname); dfname = strdup(tnam); }     }   }#endif   if( nbytes_dfname <= 0 && strstr(fname,".niml") != NULL ){     if( dfname != NULL ) free(dfname) ;     dfname = strdup(fname) ; strcpy(dfname+strlen(dfname)-5,".data") ;     nbytes_dfname = THD_filesize(dfname) ;   }   if( nbytes_dfname <= 0 ){     char mess[THD_MAX_NAME+256] ;     sprintf(mess,"datafile is missing (%s)",dfname) ; GQUIT(mess) ;   } else if( nbytes_dfname < nbytes_needed ){     char mess[THD_MAX_NAME+1024] ;     sprintf(mess,"datafile %s has %s bytes but needs at least %s",              dfname ,              commaized_integer_string(nbytes_dfname) ,              commaized_integer_string(nbytes_needed) ) ;

//.........这里部分代码省略.........   /* get the maximum integer in the unique array */   imax = 0;   for (iunq=0; iunq<N_final_unq; ++iunq) {      if (final_unq[iunq] > imax) imax = final_unq[iunq];       if (fout) fprintf(fout, "%d   %d/n", iunq, final_unq[iunq]);      hd = (INT_HASH_DATUM*)calloc(1,sizeof(INT_HASH_DATUM));      hd->id = final_unq[iunq];      hd->index = iunq;      HASH_ADD_INT(rmap, id, hd);    }      fclose(fout); fout=NULL;   /* now cycle over all dsets and replace their voxel values with rank */   for (ib = 0; ib<nbriks; ++ib) {      for (isb=0; isb<DSET_NVALS(dsets_in[ib]); ++isb) {         EDIT_BRICK_LABEL(  dsets_in[ib],isb, "rank" ) ;         EDIT_BRICK_TO_NOSTAT(  dsets_in[ib],isb ) ;         EDIT_BRICK_FACTOR( dsets_in[ib],isb, 0.0);/* no factors for rank*/         switch (DSET_BRICK_TYPE(dsets_in[ib],isb) ){            default:               fprintf(stderr,                        "** Bad dset type for unique operation./n"                        "Only Byte, Short, and float dsets are allowed./n");               break ; /* this should not happen here,                         so don't bother returning*/            case MRI_short:{               short *mar = (short *) DSET_ARRAY(dsets_in[ib],isb) ;               if (imax >  MRI_TYPE_maxval[MRI_short]) {                  WARNING_message("Maximum rank value of %d is/n"                                  "than maximum value for dset datatype of %d/n",                                  imax, MRI_TYPE_maxval[MRI_short]);               }               for( ii=0 ; ii < DSET_NVOX(dsets_in[ib]) ; ii++ )                  if (!cmask || cmask[ii]) {                     id = (int)mar[ii];                     HASH_FIND_INT(rmap,&id ,hd);                     if (hd) mar[ii] = (short)(hd->index);                      else                        ERROR_exit("** Failed to find key %d in hash table/n",id);                  } else mar[ii] = 0;            }            break ;            case MRI_byte:{               byte *mar ;               if (imax >  MRI_TYPE_maxval[MRI_short]) {                  WARNING_message("Maximum rank value of %d is/n"                                  "than maximum value for dset datatype of %d/n",                                  imax, MRI_TYPE_maxval[MRI_byte]);               }               mar = (byte *) DSET_ARRAY(dsets_in[ib],isb) ;               for( ii=0 ; ii < DSET_NVOX(dsets_in[ib]) ; ii++ )                  if (!cmask || cmask[ii]) {                     id = (int)mar[ii];                     HASH_FIND_INT(rmap,&id ,hd);                     if (hd) mar[ii] = (byte)(hd->index);                      else                        ERROR_exit("** Failed to find key %d in hash table/n",id);                  } else mar[ii] = 0;            }            break ;            case MRI_float:{               float *mar = (float *) DSET_ARRAY(dsets_in[ib],isb) ;               for( ii=0 ; ii < DSET_NVOX(dsets_in[ib]) ; ii++ )                  if (!cmask || cmask[ii]) {                     id = (int)mar[ii]; /* Assuming float is integral valued */