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