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