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

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

MRI_vectim * THD_dset_list_censored_to_vectim( int nds, THD_3dim_dataset **ds,                                               byte *mask , int nkeep , int *keep ){   MRI_vectim *vout , **vim ;   int kk , jj ;   if( nds < 1 || ds == NULL ) return NULL ;   if( nds == 1 )   /* trivial case */     return THD_dset_censored_to_vectim( ds[0],mask,nkeep,keep );   for( kk=0 ; kk < nds ; kk++ ){     if( !ISVALID_DSET(ds[kk]) ) return NULL ;     if( DSET_NVALS(ds[kk]) != DSET_NVALS(ds[0]) ) return NULL ;   }#pragma omp critical (MALLOC)   vim = (MRI_vectim **)malloc(sizeof(MRI_vectim *)*nds) ;   for( kk=0 ; kk < nds ; kk++ ){     vim[kk] = THD_dset_censored_to_vectim( ds[kk] , mask , nkeep,keep ) ;     /** DSET_unload( ds[kk] ) ; **/     if( vim[kk] == NULL ){       for( jj=0 ; jj < kk ; jj++ ) VECTIM_destroy(vim[jj]) ;       free(vim) ; return NULL ;     }   }   vout = THD_tcat_vectims( nds , vim ) ;   for( jj=0 ; jj < nds ; jj++ ) VECTIM_destroy(vim[jj]) ;   free(vim) ; return vout ;}

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

/* Calculates the average value for each voxel in dset across all timepoints.   PRE:   dset is a valid 3D+T dataset with at least 1 timepoint;   ignore is the number of timepoints to ignore at the beginning of dset;   0 <= ignore < number of timepoints in dset;   POST:   return value is NULL on error,   else, return value is an array of floats with the same dimensions as the   subbricks of dset, containing the voxel value averages;   Note: The caller is responsible for free()ing the returned block of memory   when done.   Note2: The complex datatype is not supported, and any such bricks will   result in an error (NULL return value).*/double * RIC_CalcVoxelMeans(const THD_3dim_dataset * dset, int ignore) {    double * avg;       /* The voxel averages to be returned */    float scalefactor; /* Current dset brick scaling factor */    int ival, nvals;   /* Current, number of dset timepoints */    int ivox, nvoxs;   /* Current, number of dset brick voxels */    /* Quick check of arguments */    if (!ISVALID_3DIM_DATASET(dset) || DSET_NVALS(dset) < 1 ||	ignore < 0 || ignore >= DSET_NVALS(dset)) {	return NULL;    }    /* Initialize */    DSET_load(dset);    nvals = DSET_NVALS(dset);    nvoxs = dset->daxes->nxx * dset->daxes->nyy * dset->daxes->nzz;    avg = malloc(sizeof(double) * nvoxs);    if (avg == NULL) {	return NULL;    }    /* Calculate the voxel averages; treat matrices as 1-D arrays */    /* Zero the voxel sums */    for (ivox = 0; ivox < nvoxs; ivox += 1) {	avg[ivox] = 0.0;    }    /* Sum each voxel across time (and hope there are not too many points) */    for (ival = ignore; ival < nvals; ival += 1) {	scalefactor = DSET_BRICK_FACTOR(dset, ival);	switch (DSET_BRICK_TYPE(dset, ival)) {	case MRI_short:	    RIC_CALCVOXELMEANS__DO_VOXSUM(short);	    break;	case MRI_byte:	    RIC_CALCVOXELMEANS__DO_VOXSUM(byte);	    break;	case MRI_float:	    RIC_CALCVOXELMEANS__DO_VOXSUM(float);	    break;	default: /* Unsupported datatype */	    free(avg);	    return NULL;	}    }    /* Divide by number of timepoints to get average */    nvals -= ignore;  /* We do not average over the ignored timepoints */    for (ivox = 0; ivox < nvoxs; ivox += 1) {	avg[ivox] /= nvals;    }    return avg;}

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

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

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

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

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

void THD_vectim_indexed_to_dset( MRI_vectim *mrv, int nlist, int *ilist,                                 THD_3dim_dataset *dset ){   int nvals , nvec ,  jj,kk ;   float *tar , *var ;ENTRY("THD_vectim_indexed_to_dset") ;   if( mrv   == NULL || !ISVALID_DSET(dset) ||       nlist <= 0    || ilist == NULL       || nlist > DSET_NVALS(dset)  ){     ERROR_message("THD_vectim_indexed_to_dset: illegal inputs (nlist=%d)",nlist) ;     EXRETURN ;   }   nvec  = mrv->nvec ;   nvals = mrv->nvals ;   for( kk=0 ; kk < nlist ; kk++ ){     if( ilist[kk] < 0 || ilist[kk] >= nvals ){       ERROR_message("THD_vectim_indexed_to_dset: illegal ilist[%d]=%d",kk,ilist[kk]) ;       EXRETURN ;     }   }   tar = (float *)malloc(sizeof(float)*nlist) ;   for( kk=0 ; kk < nvec ; kk++ ){     var = VECTIM_PTR(mrv,kk) ;     for( jj=0 ; jj < nlist ; jj++ ) tar[jj] = var[ilist[jj]] ;     THD_insert_series( mrv->ivec[kk] , dset ,                        nlist , MRI_float , tar , 0 ) ;   }   free(tar) ; EXRETURN ;}

void THD_vectim_to_dset( MRI_vectim *mrv , THD_3dim_dataset *dset ){   int nvals , nvec ,  kk , ign ;ENTRY("THD_vectim_to_dset") ;   if( mrv == NULL || !ISVALID_DSET(dset)           ) EXRETURN ;   if( mrv->nvals + mrv->ignore != DSET_NVALS(dset) ) EXRETURN ;   nvec  = mrv->nvec ;   nvals = mrv->nvals ;   ign   = mrv->ignore ;   if( ign == 0 ){     for( kk=0 ; kk < nvec ; kk++ )       THD_insert_series( mrv->ivec[kk] , dset ,                          nvals , MRI_float , VECTIM_PTR(mrv,kk) , 0 ) ;   } else {     float *var ;#pragma omp critical (MALLOC)     var = (float *)malloc(sizeof(float)*(nvals+ign)) ;     for( kk=0 ; kk < nvec ; kk++ ){       (void)THD_extract_array( mrv->ivec[kk] , dset , 0 , var ) ;       AAmemcpy( var+ign , VECTIM_PTR(mrv,kk) , sizeof(float)*nvals ) ;       THD_insert_series( mrv->ivec[kk] , dset ,                          nvals , MRI_float , var , 0 ) ;     }   }   EXRETURN ;}

void THD_vectim_to_dset_indexed( MRI_vectim *mrv ,                                 THD_3dim_dataset *dset , int *tlist ){   int nvals , nvec ,  jj,kk , tmax=0 ;   float *tar , *var ;ENTRY("THD_vectim_to_dset_indexed") ;   if( mrv == NULL || !ISVALID_DSET(dset) || tlist == NULL ) EXRETURN ;   nvec  = mrv->nvec ;   nvals = mrv->nvals ;   for( kk=0 ; kk < nvals ; kk++ ){     if( tlist[kk] < 0    ) EXRETURN ;     if( tlist[kk] > tmax ) tmax = tlist[kk] ;   }   tmax++ ; if( DSET_NVALS(dset) < tmax ) EXRETURN ;   tar = (float *)malloc(sizeof(float)*tmax) ;   for( kk=0 ; kk < nvec ; kk++ ){     var = VECTIM_PTR(mrv,kk) ;     for( jj=0 ; jj < tmax  ; jj++ ) tar[jj]        = 0.0f    ;     for( jj=0 ; jj < nvals ; jj++ ) tar[tlist[jj]] = var[jj] ;     THD_insert_series( mrv->ivec[kk] , dset ,                        tmax , MRI_float , tar , 0 ) ;   }   free(tar) ; EXRETURN ;}

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

int main( int argc , char * argv[] ){   THD_3dim_dataset * dset ;   int iarg , all = 0, verbose = 0, cnt = 0;   if( argc < 2 || strncmp(argv[1],"-help",4) == 0 ) Syntax() ;   iarg = 1 ;   cnt = 1;   while (cnt < argc) {      if( strncmp(argv[iarg],"-all",4) == 0 ){ all = 1 ; iarg++ ; }      else if( strncmp(argv[iarg],"-verbose",5) == 0 ){ verbose = 1 ; iarg++ ; }      ++cnt;   }   for( ; iarg < argc ; iarg++ ){      dset = THD_open_dataset( argv[iarg] ) ;      if( dset == NULL ){         printf("-1/n") ;         continue ;      }      if (!all) {         if (verbose) {            printf("%s: %d/n",                     DSET_HEADNAME(dset),                     DSET_NVALS(dset)) ;         } else {            printf("%d/n",DSET_NVALS(dset)) ;         }      } else {         if (verbose) {            printf("%s: %d %d %d %d/n",                     DSET_HEADNAME(dset),                     DSET_NX(dset), DSET_NY(dset), DSET_NZ(dset),                     DSET_NVALS(dset)) ;         } else {            printf("%d %d %d %d/n",                     DSET_NX(dset), DSET_NY(dset), DSET_NZ(dset),                     DSET_NVALS(dset)) ;         }      }      THD_delete_3dim_dataset( dset , False ) ;   }   exit(0) ;}

int is_integral_dset ( THD_3dim_dataset *dset, int check_values){   int i=0;   if(   !ISVALID_DSET(dset)  ) return(0);   for (i=0; i<DSET_NVALS(dset); ++i) {      if (!is_integral_sub_brick(dset, i, check_values)) return(0);   }   return(1);}

// use Ndim to set number of dimensions to check: 3 or 4int CompareSetDims(THD_3dim_dataset *A, THD_3dim_dataset *B, int Ndim){   int DimA[4] = {0,0,0,0}, DimB[4] = {0,0,0,0};   int i;   if ( Ndim > 4)      ERROR_exit("Bad call to CompareSetDims-- overtime!");   DimA[0] = DSET_NX(A);   DimA[1] = DSET_NY(A);    DimA[2] = DSET_NZ(A);   DimA[3] = DSET_NVALS(A);   DimB[0] = DSET_NX(B);   DimB[1] = DSET_NY(B);    DimB[2] = DSET_NZ(B);   DimB[3] = DSET_NVALS(B);   for ( i=0 ; i<Ndim ; i++)      if ( DimA[i] != DimA[i] )         ERROR_exit("Bad dimensional matching of inputs: '%s' and '%s'!",                    DSET_PREFIXSTR(A), DSET_PREFIXSTR(B));   return 0;}

THD_3dim_dataset * THD_despike9_dataset( THD_3dim_dataset *inset , byte *mask ){   THD_3dim_dataset *outset ;   MRI_vectim *mrv ;   int ii ;ENTRY("THD_despike9_dataset") ;   if( !ISVALID_DSET(inset) || DSET_NVALS(inset) < 9 ) RETURN(NULL) ;   mrv = THD_dset_to_vectim(inset,mask,0) ;  DSET_unload(inset) ;   if( mrv == NULL ) RETURN(NULL) ;   (void)THD_vectim_despike9(mrv) ;   outset = EDIT_empty_copy(inset) ;   for( ii=0 ; ii < DSET_NVALS(outset) ; ii++ )     EDIT_substitute_brick(outset,ii,MRI_float,NULL) ;   THD_vectim_to_dset(mrv,outset) ; VECTIM_destroy(mrv) ;   RETURN(outset) ;}

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

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

/*   Return the sub-brick indices referenced in str   This is a new version of MCW_get_intlist, which allows   the use of sub-brick labels       See is_in_labels for details.         ZSS Dec 09*/int * MCW_get_thd_intlist( THD_3dim_dataset *dset , char *str ){   /* test for brick_lab is not needed, and breaks sub-brick selection      of NIfTI datasets                             6 Jan 2009 [rickr]      return( (dset && dset->dblk && dset->dblk->brick_lab) ?                MCW_get_labels_intlist (dset->dblk->brick_lab,                                        DSET_NVALS(dset), str) : NULL );   */   if( !dset || !dset->dblk ) return NULL;   return MCW_get_labels_intlist(dset->dblk->brick_lab, DSET_NVALS(dset), str);}

THD_3dim_dataset * process_dataset (){  THD_3dim_dataset * new_dset = NULL;     /* output bucket dataset */  int ibrick, nbricks;                    /* sub-brick indices */  int statcode;                           /* type of stat. sub-brick */ENTRY("process_dataset") ;  /*----- Make full copy of input dataset -----*/  new_dset = EDIT_full_copy(FDR_dset, FDR_output_prefix);  /*----- Record history of dataset -----*/  tross_Copy_History( FDR_dset , new_dset ) ;  if( commandline != NULL )    {      tross_Append_History ( new_dset, commandline);      free(commandline) ;    }  /*----- Deallocate memory for input dataset -----*/     THD_delete_3dim_dataset (FDR_dset , False );  FDR_dset = NULL ;  /*----- Loop over sub-bricks in the dataset -----*/  nbricks = DSET_NVALS(new_dset);  STATUS("start loop over bricks") ;  for (ibrick = 0;  ibrick < nbricks;  ibrick++)    {      statcode = DSET_BRICK_STATCODE(new_dset, ibrick);      if (FUNC_IS_STAT(statcode) || FDR_force )	{	  /*----- Process the statistical sub-bricks -----*/	  if (!FDR_quiet)  	    printf ("ibrick = %3d   statcode = %5s /n", 		    ibrick, FUNC_prefixstr[MAX(statcode,0)]);	  process_subbrick (new_dset, ibrick);	}    }  RETURN(new_dset);}

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

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

void THD_load_tcat( THD_datablock *dblk ){   int ivout , dd , iv ;   THD_3dim_dataset *dset_in , *dset_out ;   NI_str_array *sar ;ENTRY("THD_load_tcat") ;   if( !ISVALID_DBLK(dblk) ) EXRETURN ;   dset_out = (THD_3dim_dataset *)dblk->parent ;   if( !ISVALID_DSET(dset_out) ) EXRETURN ;   sar = NI_decode_string_list( dset_out->tcat_list , "~" ) ;   if( sar == NULL ) EXRETURN ;   if( sar->num != dset_out->tcat_num ){ NI_delete_str_array(sar); EXRETURN; }   ivout = 0 ;   for( dd=0 ; dd < sar->num ; dd++ ){     dset_in = THD_open_dataset( sar->str[dd] ) ;     if( dset_in == NULL ){       NI_delete_str_array(sar) ; DSET_unload(dset_out) ;       EXRETURN ;     }     DSET_mallocize(dset_in) ; DSET_load(dset_in) ;     if( !DSET_LOADED(dset_in) ){       NI_delete_str_array(sar) ; DSET_unload(dset_out) ; DSET_delete(dset_in) ;       EXRETURN ;     }     for( iv=0 ; iv < DSET_NVALS(dset_in) ; iv++ ){       EDIT_substitute_brick( dset_out , ivout ,                              DSET_BRICK_TYPE(dset_in,iv), DSET_ARRAY(dset_in,iv) );       mri_fix_data_pointer( NULL , DSET_BRICK(dset_in,iv) ) ;       EDIT_BRICK_FACTOR( dset_out , ivout , DSET_BRICK_FACTOR(dset_in,iv) ) ;       EDIT_BRICK_LABEL(dset_out, ivout,                         DSET_BRICK_LABEL(dset_in, iv)); /* ZSS Aug. 27 2012 */       ivout++ ;     }     DSET_delete(dset_in) ;   }   NI_delete_str_array(sar) ; EXRETURN ;}

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

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