C++ DSET_load函数代码示例

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

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

/*----------------------------------------------------------------------
**
**  Main routine for this plugin (will be called from AFNI).
**
**----------------------------------------------------------------------
*/
char * HEMISUB_main( PLUGIN_interface * plint )
{
    THD_3dim_dataset * dset, * new_dset;
    MCW_idcode       * idc;
    hemi_s             hs = { 0 };
    char             * new_prefix;
    char             * ret_string = NULL;
    char             * tag;


    if ( plint == NULL )
	return  "------------------------\n"
		"HEMISUB_main: NULL input\n"
		"------------------------\n";

    PLUTO_next_option( plint );
    idc  = PLUTO_get_idcode( plint );
    dset = PLUTO_find_dset( idc );

    if( dset == NULL )
	return "-------------------------------\n"
	       "HEMISUB_main: bad input dataset\n"
	       "-------------------------------";

    DSET_load( dset );

    PLUTO_next_option( plint );
    new_prefix = PLUTO_get_string( plint );
    if ( ! PLUTO_prefix_ok( new_prefix ) )
	return  "------------------------\n"
		"HEMISUB_main: bad prefix\n"
		"------------------------\n";

    if ( ( new_dset = PLUTO_copy_dset( dset, new_prefix ) ) == NULL )
	return  "------------------------------------------\n"
		"HEMISUB_main: failed to copy input dataset\n"
		"------------------------------------------\n";

    tag = PLUTO_get_optiontag( plint );
    if ( tag && ! strcmp( tag, "Thresh Type" ) )
    {
	tag = PLUTO_get_string( plint );
	if ( tag != NULL )
	    hs.thresh_type = PLUTO_string_index( tag, NUM_T_OPTS, thresh_opts );
    }

    if ( ret_string = process_data( new_dset, &hs ) )
	return  ret_string;

    if ( PLUTO_add_dset( plint, new_dset, DSET_ACTION_MAKE_CURRENT ) )
    {
	THD_delete_3dim_dataset( new_dset, False );

	return  "---------------------------------------\n"
		"HEMISUB_main: failed to add new dataset\n"
		"---------------------------------------\n";
    }

    return NULL;
}

/* - moved from 3dhistog.c               18 Dec 2012 [rickr] */
int THD_slow_minmax_dset(THD_3dim_dataset *dset, float *dmin, float *dmax,
                         int iv_bot, int iv_top)
{
   int iv_fim;
   float fimfac;
   float vbot , vtop, temp_fbot=1.0, temp_ftop=0.0 ;

   DSET_load(dset);

   for( iv_fim=iv_bot ; iv_fim <= iv_top ; iv_fim++ ){
     /* minimum and maximum for sub-brick */
     vbot = mri_min( DSET_BRICK(dset,iv_fim) ) ;
     vtop = mri_max( DSET_BRICK(dset,iv_fim) ) ;
     fimfac = DSET_BRICK_FACTOR(dset,iv_fim) ;
     if (fimfac == 0.0)  fimfac = 1.0;
     vbot *= fimfac ; vtop *= fimfac ;

     /* update global min and max */
     if ( temp_fbot > temp_ftop ) { /* first time, just copy */
        temp_fbot = vbot;
        temp_ftop = vtop;
     } else {
        if( vbot < temp_fbot ) temp_fbot = vbot;
        if( vtop > temp_ftop ) temp_ftop = vtop;
     }
   }
   *dmin = temp_fbot;
   *dmax = temp_ftop;

   return(0);
}

int main( int argc , char *argv[] )
{
   THD_3dim_dataset *dset ;
   MRI_IMAGE *im ;
   int iarg=1 ;

   if( strcmp(argv[1],"-nper") == 0 ){
      nper = strtol( argv[2] , NULL , 10 ) ;
      iarg = 3 ;
   }

   for( ; iarg < argc ; iarg++ ){
     printf("======= dataset %s  nper=%d ========\n",argv[iarg],nper) ;
     dset = THD_open_dataset(argv[iarg]) ;
     if( dset == NULL ) continue ;
     DSET_load(dset) ;
     im = DSET_BRICK(dset,0) ;
     im->dx = DSET_DX(dset) ;
     im->dy = DSET_DY(dset) ;
     im->dz = DSET_DZ(dset) ;
     (void) THD_orient_guess( im ) ;

     printf( "Data Axes Orientation:\n"
             "  first  (x) = %s\n"
             "  second (y) = %s\n"
             "  third  (z) = %s\n" ,
        ORIENT_typestr[dset->daxes->xxorient] ,
        ORIENT_typestr[dset->daxes->yyorient] ,
        ORIENT_typestr[dset->daxes->zzorient]  ) ;
     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) ;
}

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

double ENTROPY_dataset( THD_3dim_dataset *dset )
{
   if( !ISVALID_DSET(dset) ) return(0.0) ;
   DSET_load(dset) ;
   if( !DSET_LOADED(dset) ) return(0.0) ;
   return ENTROPY_datablock( dset->dblk ) ;
}

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

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

void get_options
(
  int argc,                        /* number of input arguments */
  char ** argv                     /* array of input arguments */ 
)

{
  int nopt = 1;                     /* input option argument counter */
  int ival, index;                  /* integer input */
  float fval;                       /* float input */
  char message[MAX_STRING_LENGTH];  /* error message */


  /*----- does user request help menu? -----*/
  if (argc < 2 || strncmp(argv[1], "-help", 5) == 0)  display_help_menu();  
   

  /*----- main loop over input options -----*/
  while (nopt < argc )
    {

      /*-----   -anat filename   -----*/
      if (strncmp(argv[nopt], "-anat", 5) == 0)
	{
	  nopt++;
	  if (nopt >= argc)  estPDF_error ("need argument after -anat ");
	  anat_filename = malloc (sizeof(char) * MAX_STRING_LENGTH);
	  MTEST (anat_filename);
	  strcpy (anat_filename, argv[nopt]);

	  anat = THD_open_one_dataset (anat_filename);
	  if (!ISVALID_3DIM_DATASET (anat))
	    {
	      sprintf (message, "Can't open dataset: %s\n", anat_filename); 
	      estPDF_error (message); 
	    } 
	  DSET_load(anat); CHECK_LOAD_ERROR(anat);

	  nopt++;
	  continue;
	}
      

      /*----- unknown command -----*/
      sprintf(message,"Unrecognized command line option: %s\n", argv[nopt]);
      estPDF_error (message);
      
    }

  
}

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

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