/** \file volume.c
 * \brief MINC 2.0 Volume Functions
 * \author Leila Baghdadi, Bert Vincent
 *
 * Functions to create, open, and close MINC volume objects.
 ************************************************************************/
#include <stdlib.h>
#include <hdf5.h>
#include <limits.h>
#include <float.h>
#include <math.h>
#include "minc2.h"
#include "minc2_private.h"

/**
 * \defgroup mi2Vol MINC 2.0 Volume Functions
 */

/* Forward declarations */

static void miinit_default_range(mitype_t mitype, double *valid_max, 
                                 double *valid_min);
static void miread_valid_range(mihandle_t volume, double *valid_max, 
                               double *valid_min);

static int _miset_volume_class(mihandle_t volume, miclass_t volclass);
static int _miget_volume_class(mihandle_t volume, miclass_t *volclass);

/** Create the actual image for the volume.
    Note that the image dataset muct be created in the hierarchy 
    before the image data can be added.
    \ingroup mi2Vol
 */
int
micreate_volume_image(mihandle_t volume)
{
    char dimorder[MI2_CHAR_LENGTH];
    int i;
    hid_t dataspace_id;
    hid_t dset_id;
    hsize_t hdf_size[MI2_MAX_VAR_DIMS];

    /* Try creating IMAGE dataset i.e. /minc-2.0/image/0/image
     */

    dimorder[0] = '\0'; /* Set string to empty */
    
    for (i = 0; i < volume->number_of_dims; i++) {
        hdf_size[i] = volume->dim_handles[i]->length;

        /* Create the dimorder string, ordered comma-separated
           list of dimension names.
        */
        strcat(dimorder, volume->dim_handles[i]->name);
        if (i != volume->number_of_dims - 1) {
            strcat(dimorder, ",");
        }
    }
    
   
    /* Create a SIMPLE dataspace  */
    dataspace_id = H5Screate_simple(volume->number_of_dims, hdf_size, NULL);
    if (dataspace_id < 0) {
        return (MI_ERROR);
    }
    
    dset_id = H5Dcreate1(volume->hdf_id, "/minc-2.0/image/0/image", 
                        volume->ftype_id, 
                        dataspace_id, volume->plist_id);
    
    if (dset_id < 0) {  
        return (MI_ERROR);
    }

    volume->image_id = dset_id;
    
    hdf_var_declare(volume->hdf_id, "image", "/minc-2.0/image/0/image",
                    volume->number_of_dims, hdf_size);
    /* Create the dimorder attribute, ordered comma-separated
       list of dimension names.
    */
    miset_attr_at_loc(dset_id, "dimorder", MI_TYPE_STRING,
                      strlen(dimorder), dimorder);

    H5Sclose(dataspace_id);
    
    if (volume->volume_class == MI_CLASS_REAL) {
        int ndims;
        hid_t dcpl_id;
        double dtmp;

        dcpl_id = H5Pcreate(H5P_DATASET_CREATE);
        if (dcpl_id < 0) {
            return (MI_ERROR);
        }

        if (volume->has_slice_scaling) {
            /* TODO: Find the slowest-varying spatial dimension; that forms
             * the basis for the image-min and image-max variables.  Right
             * now this is an oversimplification!
             */
            ndims = volume->number_of_dims - 2;
            dataspace_id = H5Screate_simple(ndims, hdf_size, NULL);
        }
        else {
            ndims = 0;
            dataspace_id = H5Screate(H5S_SCALAR);
        }

        if (ndims != 0) {
            dimorder[0] = '\0'; /* Set string to empty */
            for (i = 0; i < ndims; i++) {
                /* Create the dimorder string, ordered comma-separated
                   list of dimension names.
                */
                strcat(dimorder, volume->dim_handles[i]->name);
                if (i != volume->number_of_dims - 1) {
                    strcat(dimorder, ",");
                }
            }
        }

	/* Create the image minimum dataset for FULL-RESOLUTION storage of data
	 */
        dtmp = 0.0;
        H5Pset_fill_value(dcpl_id, H5T_NATIVE_DOUBLE, &dtmp);

        dset_id = H5Dcreate1(volume->hdf_id, "/minc-2.0/image/0/image-min",
                            H5T_IEEE_F64LE, dataspace_id, dcpl_id);
        if (ndims != 0) {
            miset_attr_at_loc(dset_id, "dimorder", MI_TYPE_STRING,
                              strlen(dimorder), dimorder);
        }
        volume->imin_id = dset_id;
        hdf_var_declare(volume->hdf_id, "image-min", "/minc-2.0/image/0/image-min", ndims, hdf_size);

        
	/* Create the image maximum dataset for FULL-RESOLUTION storage of data
	 */
        dtmp = 1.0;
        H5Pset_fill_value(dcpl_id, H5T_NATIVE_DOUBLE, &dtmp);

        dset_id = H5Dcreate1(volume->hdf_id, "/minc-2.0/image/0/image-max",
                            H5T_IEEE_F64LE, dataspace_id, dcpl_id);
        if (ndims != 0) {
            miset_attr_at_loc(dset_id, "dimorder", MI_TYPE_STRING,
                              strlen(dimorder), dimorder);
        }
        volume->imax_id = dset_id;
        hdf_var_declare(volume->hdf_id, "image-max", "/minc-2.0/image/0/image-max", ndims, hdf_size);
        H5Sclose(dataspace_id);
        H5Pclose(dcpl_id);
    }

    return (MI_NOERROR);
}

/** Set up the array of conversions from voxel to world coordinate order.
 */
static int
miset_volume_world_indices(mihandle_t hvol)
{
    int i;

    for (i = 0; i < hvol->number_of_dims; i++) {
        midimhandle_t hdim = hvol->dim_handles[i];

        hdim->world_index = -1;
        if (hdim->class == MI_DIMCLASS_SPATIAL) {
            if (!strcmp(hdim->name, MIxspace)) {
                hdim->world_index = MI2_X;
            }
            else if (!strcmp(hdim->name, MIyspace)) {
                hdim->world_index = MI2_Y;
            }
            else if (!strcmp(hdim->name, MIzspace)) {
                hdim->world_index = MI2_Z;
            }
        }
        else if (hdim->class == MI_DIMCLASS_SFREQUENCY) {
            if (!strcmp(hdim->name, MIxfrequency)) {
                hdim->world_index = MI2_X;
            }
            else if (!strcmp(hdim->name, MIyfrequency)) {
                hdim->world_index = MI2_Y;
            }
            else if (!strcmp(hdim->name, MIzfrequency)) {
                hdim->world_index = MI2_Z;
            }
        }
    }
    return (MI_NOERROR);
}

/** Create and initialize a MINC 2.0 volume structure.
 */
static mihandle_t 
mialloc_volume_handle(void)
{
    mihandle_t handle = (mihandle_t) malloc(sizeof(struct mivolume));

    if (handle != NULL) {
        /* Clear the memory by default. */
        memset(handle, 0, sizeof(struct mivolume));

        /* Set the defaults for the data structure */
        handle->scale_min = 0.0;
        handle->scale_max = 1.0;
        handle->image_id = -1;
        handle->imax_id = -1;
        handle->imin_id = -1;
        handle->plist_id = -1;
        handle->has_slice_scaling = FALSE;
        handle->is_dirty = FALSE;
        handle->dim_indices = NULL;
        handle->selected_resolution = 0;
    }
    return (handle);
}

/** Create a volume with the specified name, dimensions, 
    type, class, volume properties and retrieve the volume handle.
    \ingroup mi2Vol
 */ 
int
micreate_volume(const char *filename, int number_of_dimensions,
		midimhandle_t dimensions[], mitype_t volume_type,
		miclass_t volume_class, mivolumeprops_t create_props,
		mihandle_t *volume)
{
  int i;
  int stat;
  hid_t file_id;
  hid_t hdf_type;
  hid_t hdf_plist;
  hid_t fspc_id;
  hsize_t dim[1];
  hid_t grp_id; 
  herr_t status;
  hid_t dataset_id = -1;
  hid_t dataset_width = -1;
  hid_t dataspace_id = -1;
  char *name;
  int size;
  hsize_t hdf_size[MI2_MAX_VAR_DIMS];
  mihandle_t handle;
  mivolumeprops_t props_handle;
  char ident_str[128];
  hid_t tmp_type;

  /* Initialization. 
     For the actual body of this function look at m2utils.c
   */
  miinit();

  /* Validate the parameters. 
   */
  if (filename == NULL) {
    return (MI_ERROR);
  }

  if (dimensions == NULL && number_of_dimensions != 0) {
      return (MI_ERROR);
  }

  /* Allocate space for the volume handle
   */
  handle = mialloc_volume_handle();
  if (handle == NULL) {
    return (MI_ERROR);
  }
  
  /* Initialize some of the variables associated with the volume handle.
   */
  handle->mode = MI2_OPEN_RDWR;
  handle->number_of_dims = number_of_dimensions;

  /* convert minc type to hdf type
   */
  hdf_type = mitype_to_hdftype(volume_type, FALSE);

  /* Setting up volume type_id 
   */
  switch (volume_class) {
  case MI_CLASS_REAL:
  case MI_CLASS_INT:
      handle->ftype_id = hdf_type;
      handle->mtype_id = H5Tget_native_type(handle->ftype_id, 
                                            H5T_DIR_ASCEND);
      break;

  case MI_CLASS_LABEL:
      /* A volume of class LABEL must have an integer type (positive).
       */
      switch (volume_type) {
      case MI_TYPE_UBYTE:
      case MI_TYPE_USHORT:
      case MI_TYPE_UINT:
      case MI_TYPE_BYTE:
      case MI_TYPE_SHORT:
      case MI_TYPE_INT:
          handle->ftype_id = H5Tenum_create(hdf_type);
          if (handle->ftype_id < 0) {
              return (MI_ERROR);
          }

          tmp_type = H5Tget_native_type(hdf_type, H5T_DIR_ASCEND);
          H5Tclose(hdf_type);
          hdf_type = tmp_type;

          /* Create an enumerated type with the native type as it's base.
           */
          handle->mtype_id = H5Tenum_create(hdf_type);
          if (handle->mtype_id < 0) {
              return (MI_ERROR);
          }
          H5Tclose(hdf_type);

          miinit_enum(handle->ftype_id);
          miinit_enum(handle->mtype_id);
          break;
      default:
          return (MI_ERROR);
      }
      break;

  case MI_CLASS_COMPLEX:
      switch (volume_type) {
      case MI_TYPE_SCOMPLEX:
      case MI_TYPE_ICOMPLEX:
      case MI_TYPE_FCOMPLEX:
      case MI_TYPE_DCOMPLEX:
          handle->ftype_id = hdf_type;
          handle->mtype_id = mitype_to_hdftype(volume_type, TRUE);
          break;
      default:
          return (MI_ERROR);
      }
      break;

  case MI_CLASS_UNIFORM_RECORD:
      handle->ftype_id = H5Tcreate(H5T_COMPOUND, H5Tget_size(hdf_type));
      handle->mtype_id = H5Tcreate(H5T_COMPOUND, H5Tget_size(hdf_type));
      H5Tclose(hdf_type);
      break;

  default:
      return (MI_ERROR);
  }

  handle->volume_class = volume_class;

  /* Create file in HDF5 with the given filename and
     H5F_ACC_TRUNC: Truncate file, if it already exists, 
                   erasing all data previously stored in the file.
     and create ID and ID access as default.
  */
  
  file_id = hdf_create(filename, H5F_ACC_TRUNC, NULL);
  if (file_id < 0) {
    return (MI_ERROR);
  }

  handle->hdf_id = file_id;

  micreate_ident(ident_str, sizeof(ident_str));
  miset_attribute(handle, MI_ROOT_PATH, "ident", MI_TYPE_STRING,
                  strlen(ident_str), ident_str);
  miset_attribute(handle, MI_ROOT_PATH, "minc_version", MI_TYPE_STRING,
                  strlen(VERSION), VERSION);

  _miset_volume_class(handle, handle->volume_class);

  /* Create a new property list for the volume
   */
  hdf_plist = H5Pcreate(H5P_DATASET_CREATE);
  if (hdf_plist < 0) {
    return (MI_ERROR);
  }

  handle->plist_id = hdf_plist;

  /* Set fill value to guarantee valid data on incomplete datasets.
   */
  if (volume_class != MI_CLASS_LABEL &&
      volume_class != MI_CLASS_UNIFORM_RECORD) {
    size_t siz = H5Tget_size(handle->ftype_id);
    char *tmp = calloc(1, siz);
    H5Pset_fill_value(hdf_plist, handle->ftype_id, tmp);
    free(tmp);
  }
  
  /* See if chunking and/or compression should be enabled 
     and if yes set the type of storage used to store the 
     raw data for a dataset.
   */
  
  if (create_props != NULL &&
      (create_props->compression_type == MI_COMPRESS_ZLIB ||
       create_props->edge_count != 0)) {
      /* Set the storage to CHUNKED */
      stat = H5Pset_layout(hdf_plist, H5D_CHUNKED);
      if (stat < 0) {
          return (MI_ERROR);
      }
      /* Create an array, hdf_size, containing the size of each chunk 
       */
      for (i=0; i < number_of_dimensions; i++) {
          hdf_size[i] = create_props->edge_lengths[i];
	   /* If the size of each chunk is greater than the size of
	     the corresponding dimension, set the chunk size to the
	     dimension size
	  */
          if (hdf_size[i] > dimensions[i]->length) {
              hdf_size[i] = dimensions[i]->length;
          }
      }
    
      /* Sets the size of the chunks used to store a chunked layout dataset */
      stat = H5Pset_chunk(hdf_plist, number_of_dimensions, hdf_size);
      if (stat < 0) {
          return (MI_ERROR);
      }
      /* Sets compression method and compression level */
      stat = H5Pset_deflate(hdf_plist, create_props->zlib_level);
      if (stat < 0) {
          return (MI_ERROR);
      }
      
  }
  else { /* No COMPRESSION or CHUNKING is enabled */
    stat = H5Pset_layout(hdf_plist, H5D_CONTIGUOUS); /*  CONTIGUOUS data */
    if (stat < 0) {
      return (MI_ERROR);
    }
  }
  
  /* See if Multi-res is set to a level above 0 and if yes create subgroups
     i.e., /minc-2.0/image/1/..
           /minc-2.0/image/2/..  etc
   */
  // must add some code to make sure that the res level is possible
  if (create_props != NULL && create_props->depth > 0) {
    for (i=0; i < create_props->depth ; i++) {
      if (minc_create_thumbnail(handle, i+1) < 0) {
	return (MI_ERROR);
      }
    }
  }

  /* Try creating DIMENSIONS GROUP i.e. /minc-2.0/dimensions
   */
  grp_id = H5Gopen1(file_id, MI_FULLDIMENSIONS_PATH);
  if (grp_id < 0) {
      return (MI_ERROR);
  }
  /* Once the DIMENSIONS GROUP is opened, create each dimension.
  */
  
  for (i=0; i < number_of_dimensions ; i++) {
    
    /* First create the dataspace required to create a 
       dimension variable (dataset)
     */
    if (dimensions[i]->attr & MI_DIMATTR_NOT_REGULARLY_SAMPLED) {
      dim[0] = dimensions[i]->length;
      dataspace_id = H5Screate_simple(1, dim, NULL);
    }
    else {
      dataspace_id = H5Screate(H5S_SCALAR);
    }

    if (dataspace_id < 0) {
      return (MI_ERROR);
    }
     
    /* Create a dataset(dimension variable name) in DIMENSIONS GROUP */
    dataset_id = H5Dcreate1(grp_id, dimensions[i]->name, 
                           H5T_IEEE_F64LE, dataspace_id, H5P_DEFAULT);
    
    /* Dimension variable for a regular dimension contains
       no meaningful data. Whereas, Dimension variable for 
       an irregular dimension contains a vector with the lengths 
       equal to the sampled points along the dimension.
       Also, create a variable named "<dimension>-width" and 
       write the dimension->widths.
     */
    
    /* Check for irregular dimension and make sure
       offset values are provided for this dimension 
     */
    
    if (dimensions[i]->attr & MI_DIMATTR_NOT_REGULARLY_SAMPLED) {
      if (dimensions[i]->offsets == NULL) {
	return (MI_ERROR);
      }
      else {
	
	/* If dimension is regularly sampled */
	fspc_id = H5Dget_space(dataset_id);
	if (fspc_id < 0) {
	  return (MI_ERROR);
	}
	/* Write the raw data from buffer (dimensions[i]->offsets)
	   to the dataset.
	 */
	status = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, dataspace_id, 
                          fspc_id, H5P_DEFAULT, dimensions[i]->offsets);
	if (status < 0) {
	  return (MI_ERROR);
	}
	/* Write the raw data from buffer (dimensions[i]->offsets)
	   to the dataset.
	 */
	size = strlen(dimensions[i]->name) + 6 + 1;
        name = malloc(size);
        strcpy(name, dimensions[i]->name);
	strcat(name, "-width");

	/* Create dataset dimension_name-width */
	dataset_width = H5Dcreate1(grp_id, name, H5T_IEEE_F64LE, 
                                  dataspace_id, H5P_DEFAULT);
	/* Return an Id for the dataspace of the dataset dataset_width */
	fspc_id = H5Dget_space(dataset_width);
	if (fspc_id < 0) {
	  return (MI_ERROR);
	}
	/* Write the raw data from buffer (dimensions[i]->widths)
	   to the dataset.
	 */
	status = H5Dwrite(dataset_width, H5T_NATIVE_DOUBLE, dataspace_id, fspc_id, H5P_DEFAULT, dimensions[i]->widths);
	if (status < 0) {
	  return (MI_ERROR);
	}
	/* Create new attribute "length", with appropriate 
	   type (to hdf5) conversion.
	   miset_attr_at_loc(..) is implemented at m2utils.c
	*/
        miset_attr_at_loc(dataset_width, "length", MI_TYPE_INT,
                          1, &dimensions[i]->length);
	/* Close the specified datatset */
	H5Dclose(dataset_width);
        free(name);
      }
    }
    
    if (dimensions[i]->attr & MI_DIMATTR_NOT_REGULARLY_SAMPLED) {
        name = "irregular";
    }
    else {
        name = "regular__";
    }
     /* Create attribute "spacing" and set its value to 
       "regular__" or "irregular"
     */
    miset_attr_at_loc(dataset_id, "spacing", MI_TYPE_STRING,
                      strlen(name), name);
    
    switch (dimensions[i]->class) {
    case MI_DIMCLASS_SPATIAL:
        name = "spatial";
        break;
    case MI_DIMCLASS_TIME:
        name = "time___";
        break;
    case MI_DIMCLASS_SFREQUENCY:
        name = "sfreq__";
        break;
    case MI_DIMCLASS_TFREQUENCY:
        name = "tfreq__";
        break;
    case MI_DIMCLASS_USER:
        name = "user___";
        break;
    case MI_DIMCLASS_RECORD:
        name = "record_";
        break;
    case MI_DIMCLASS_ANY:
    default:
        /* These should not be seen in this context!!!
         */
        return (MI_ERROR);
    }
     
    miset_attr_at_loc(dataset_id, "class", MI_TYPE_STRING, strlen(name),
                      name);

   /* Create Dimension attribute "direction_cosines"  */
    miset_attr_at_loc(dataset_id, "direction_cosines", MI_TYPE_DOUBLE,
                      3, dimensions[i]->direction_cosines);
    
    /* Save dimension length */
    miset_attr_at_loc(dataset_id, "length", MI_TYPE_INT,
                      1, &dimensions[i]->length);

    /* Save step value. */
    miset_attr_at_loc(dataset_id, "step", MI_TYPE_DOUBLE,
                      1, &dimensions[i]->step);

    /* Save start value. */
    miset_attr_at_loc(dataset_id, "start", MI_TYPE_DOUBLE,
                      1, &dimensions[i]->start);
    
    /* Save units. */
    miset_attr_at_loc(dataset_id, "units", MI_TYPE_STRING,
                      strlen(dimensions[i]->units), dimensions[i]->units);
    
    /* Save sample width. */
    miset_attr_at_loc(dataset_id, "width", MI_TYPE_DOUBLE,
                      1,  &dimensions[i]->width);

    /* Save comments. If user has not specified
       any comments, do not add this attribute
     */
 
    if (dimensions[i]->comments != NULL) {
      miset_attr_at_loc(dataset_id, "comments", MI_TYPE_STRING,
			strlen(dimensions[i]->comments),
                        dimensions[i]->comments);
    }
 
    /* Close the dataset with the specified Id
     */
    H5Dclose(dataset_id);
    
    
  } //for (i=0; i < number_of_dimensions ; i++) 
  
  /* Close the group with the specified Id
   */
  H5Gclose(grp_id);

  /* Allocate space for all the dimension handles
     Note, each volume handle is associated with an array of
     dimension handles in the order that they were create (i.e, file order)
   */
  handle->dim_handles = (midimhandle_t *)malloc(number_of_dimensions * 
                                                sizeof(midimhandle_t));  
   
  if (handle->dim_handles == NULL) {
    return (MI_ERROR);
  }

  /* Once the space for all dimension handles is created
     fill the dimension handle array with the appropriate 
     dimension handle.
   */
  for (i = 0; i < number_of_dimensions; i++) {
    handle->dim_handles[i] = dimensions[i];
    dimensions[i]->volume_handle = handle;
  }

  miset_volume_world_indices(handle);

  /* Verify the volume type.
   */
  switch (volume_type) {
  case MI_TYPE_BYTE:
  case MI_TYPE_SHORT:
  case MI_TYPE_INT:
  case MI_TYPE_FLOAT:
  case MI_TYPE_DOUBLE:
  case MI_TYPE_STRING:
  case MI_TYPE_UBYTE:
  case MI_TYPE_USHORT:
  case MI_TYPE_UINT:
  case MI_TYPE_SCOMPLEX:
  case MI_TYPE_ICOMPLEX:
  case MI_TYPE_FCOMPLEX:
  case MI_TYPE_DCOMPLEX:
  case MI_TYPE_UNKNOWN:
    break;
  default:
    return (MI_ERROR);
  }

  handle->volume_type = volume_type;

  /* Set the initial value of the valid-range 
   */
  miinit_default_range(handle->volume_type,
                       &handle->valid_max, 
                       &handle->valid_min);

  /* Get the voxel to world transform for the volume
   */
  miget_voxel_to_world(handle, handle->v2w_transform);

  /* Calculate the inverse transform */
  miinvert_transform(handle->v2w_transform, handle->w2v_transform);

  /* Allocated space for the volume properties */
  props_handle = (mivolumeprops_t)malloc(sizeof(struct mivolprops));
  /* Initialize volume properties with zero */
  memset(props_handle, 0, sizeof (struct mivolprops));
  /* If volume properties is specified by the user
     set all the properties of the volume handle
   */
  if (create_props != NULL) {
       /* Set the enable_flag for multi-resolution */
      props_handle->enable_flag = create_props->enable_flag;
      /* Set the depth of multi-resolution, i.e., how many
	 levels of resolution is specified maximum is 16.
       */
      props_handle->depth = create_props->depth;
      /* Set compression type, currently two values
	 either no compression or zlib is applicable.
       */
      switch (create_props->compression_type) {
      case MI_COMPRESS_NONE:
          props_handle->compression_type = MI_COMPRESS_NONE;
          break;
      case MI_COMPRESS_ZLIB:
          props_handle->compression_type = MI_COMPRESS_ZLIB;
          break;
      default:
          return (MI_ERROR);
      }
      /* Note that setting compression on (i.e., MI_COMPRESS_ZLIB)
	 turns chunking on by default. Need to set the number of chunks
	 (edge_count)
      */
      props_handle->zlib_level = create_props->zlib_level;
      props_handle->edge_count = create_props->edge_count;
      /* Allocate space for an array which holds the size of each chunk
	 and fill the array with the appropriiate chunk sizes.
       */
      props_handle->edge_lengths = (int *)malloc(create_props->max_lengths*sizeof(int));
      for (i=0;i<create_props->max_lengths; i++) {
          props_handle->edge_lengths[i] = create_props->edge_lengths[i];
      }
 
      props_handle->max_lengths = create_props->max_lengths;
      props_handle->record_length = create_props->record_length;

      /* Explicitly allocate storage for name
       */
      if (create_props->record_name != NULL) {
          props_handle->record_name =malloc(strlen(create_props->record_name) + 1 );
          strcpy(props_handle->record_name, create_props->record_name);
      }
      props_handle->template_flag = create_props->template_flag;
  }
  /* Set the handle to volume properties */
  handle->create_props = props_handle;
  /* Return volume handle */
  *volume = handle;
 
  return (MI_NOERROR);
}

/** Return the number of dimensions associated with this volume.
    \ingroup mi2Vol
 */
int
miget_volume_dimension_count(mihandle_t volume, midimclass_t cls,
			     midimattr_t attr, int *number_of_dimensions)
{
    int i, count=0;
    /* Validate the parameters */
    if (volume == NULL || number_of_dimensions == NULL) {
        return (MI_ERROR);
    }
    /* For each dimension check to make sure that dimension class and 
       attribute match with the specified parameters and if yes 
       increment the dimension count
     */
    for (i=0; i< volume->number_of_dims; i++) {
        if ((cls == MI_DIMCLASS_ANY || volume->dim_handles[i]->class == cls) &&
            (attr == MI_DIMATTR_ALL || volume->dim_handles[i]->attr == attr)) {
            count++;
        }
    }

    *number_of_dimensions = count;
    return (MI_NOERROR);
}

/** Returns the number of voxels in the volume.
    \ingroup mi2Vol
 */
int
miget_volume_voxel_count(mihandle_t volume, int *number_of_voxels)
{
    char path[MI2_MAX_PATH];
    hid_t dset_id;
    hid_t fspc_id;

    /* Validate parameters */
    if (volume == NULL || number_of_voxels == NULL) {
        return (MI_ERROR);
    }

    /* Quickest way to do this is with the dataspace identifier of the
     * volume. Use the volume's current resolution.
     */
    sprintf(path, "/minc-2.0/image/%d/image", volume->selected_resolution);
    /* Open the dataset with the specified path
     */
    dset_id = H5Dopen1(volume->hdf_id, path);
    if (dset_id < 0) {
        return (MI_ERROR);
    }
    /* Get an Id to the copy of the dataspace */ 
    fspc_id = H5Dget_space(dset_id);
    if (fspc_id < 0) {
        return (MI_ERROR);
    }
    /* Determines the number of elements in the dataspace and
       cast the result to an integer.
     */
    *number_of_voxels = (int) H5Sget_simple_extent_npoints(fspc_id);
    /* Close the dataspace */
    H5Sclose(fspc_id);
    /* Close the dataset */
    H5Dclose(dset_id);
    return (MI_NOERROR);
}

/* Get the number of dimensions in the file */
static int 
_miget_file_dimension_count(hid_t file_id)
{
    hid_t dset_id;
    hid_t space_id;
    int result = -1;
    /* hdf5 macro can temporarily disable the automatic error printing */
    H5E_BEGIN_TRY {
        dset_id = midescend_path(file_id, "/minc-2.0/image/0/image");
    } H5E_END_TRY;

    if (dset_id >= 0) {
      /* Get an Id to the copy of the dataspace */
        space_id = H5Dget_space(dset_id);
        if (space_id > 0) {
	    /* Determine the dimensionality of the dataspace */
            result = H5Sget_simple_extent_ndims(space_id);
	    /* Close the dataspace */
            H5Sclose(space_id);
        }
	/* Close the dataset */
        H5Dclose(dset_id);
    }
    return (result);
}

/* Get dimension variable attributes for the given dimension name */
static int
_miset_volume_class(mihandle_t volume, miclass_t volume_class)
{
    const char *class_ptr;

    switch (volume_class) {
    case MI_CLASS_REAL:
        class_ptr = "real___";
        break;
    case MI_CLASS_INT:
        class_ptr = "integer";
        break;
    case MI_CLASS_LABEL:
        class_ptr = "label__";
        break;
    case MI_CLASS_COMPLEX:
        class_ptr = "complex";
        break;
    case MI_CLASS_UNIFORM_RECORD:
        class_ptr = "array__";
        break;
    default:
        return (MI_ERROR);
    }
    miset_attribute(volume, MI_ROOT_PATH, "class", MI_TYPE_STRING,
                    strlen(class_ptr), class_ptr);
    return (MI_NOERROR);
}

static int
_miget_volume_class(mihandle_t volume, miclass_t *volume_class)
{
    char class_buf[MI2_CHAR_LENGTH];

    miget_attribute(volume, MI_ROOT_PATH, "class", MI_TYPE_STRING,
                    MI2_CHAR_LENGTH, class_buf);

    if (!strcmp(class_buf, "label__")) {
        *volume_class = MI_CLASS_LABEL;
    }
    else if (!strcmp(class_buf, "integer")) {
        *volume_class = MI_CLASS_INT;
    }
    else if (!strcmp(class_buf, "complex")) {
        *volume_class = MI_CLASS_COMPLEX;
    }
    else if (!strcmp(class_buf, "array__")) {
        *volume_class = MI_CLASS_UNIFORM_RECORD;
    }
    else {
        *volume_class = MI_CLASS_REAL;
    }
    return (MI_NOERROR);
}

/* Get dimension variable attributes for the given dimension name */
static int
_miget_file_dimension(mihandle_t volume, const char *dimname, 
                      midimhandle_t *hdim_ptr)
{
    char path[MI2_CHAR_LENGTH];
    char temp[MI2_CHAR_LENGTH];
    midimhandle_t hdim;

    /* Create a path with the dimension name */
    sprintf(path, "/minc-2.0/dimensions/%s", dimname);
    /* Allocate space for the dimension handle */
    hdim = (midimhandle_t) malloc(sizeof (*hdim));
    /* Initialize everything to zero */
    memset(hdim, 0, sizeof (*hdim));

    hdim->name = strdup(dimname);
    
    /* hdf5 macro can temporarily disable the automatic error printing */
    H5E_BEGIN_TRY {
        int r;
	/* Get the attribute (spacing) from a minc file */
        r = miget_attribute(volume, path, "spacing", MI_TYPE_STRING, 
                            MI2_CHAR_LENGTH, temp);
        if (!strcmp(temp, "irregular")) {
            hdim->attr |= MI_DIMATTR_NOT_REGULARLY_SAMPLED;
        }
        else {
            hdim->attr |= MI_DIMATTR_REGULARLY_SAMPLED;
        }

	/* Get the attribute (class) from a minc file */
        r = miget_attribute(volume, path, "class", MI_TYPE_STRING, 
                            MI2_CHAR_LENGTH, temp);
        if (r < 0) {
            /* Get the default class. */
            if (!strcmp(dimname, "time")) {
                hdim->class = MI_DIMCLASS_TIME;
            }
            else if (!strcmp(dimname, "vector_dimension")) {
	      hdim->class = MI_DIMCLASS_RECORD;
	      hdim->step = 0.0;
	    }
	    else {
                hdim->class =  MI_DIMCLASS_SPATIAL;
            }
        }
        else {
            if (!strcmp(temp, "spatial")) {
                hdim->class = MI_DIMCLASS_SPATIAL;
            }
            else if (!strcmp(temp, "time___")) {
                hdim->class = MI_DIMCLASS_TIME;
            }
            else if (!strcmp(temp, "sfreq__")) {
                hdim->class = MI_DIMCLASS_SFREQUENCY;
            }
            else if (!strcmp(temp, "tfreq__")) {
                hdim->class = MI_DIMCLASS_TFREQUENCY;
            }
            else if (!strcmp(temp, "user___")) {
                hdim->class = MI_DIMCLASS_USER;
            }
            else if (!strcmp(temp, "record_")) {
                hdim->class = MI_DIMCLASS_RECORD;
            }
            else {
                /* TODO: error message?? */
            }
        }
	/* Get the attribute (length) from a minc file */
        r = miget_attribute(volume, path, "length", MI_TYPE_UINT, 1, &hdim->length);
        if (r < 0) {
            fprintf(stderr, "Can't get length\n");
        }
	/* Get the attribute (start) from a minc file for NON vector_dimension only */
	if (strcmp(dimname, "vector_dimension")) {
	  r = miget_attribute(volume, path, "start", MI_TYPE_DOUBLE, 1, &hdim->start);
	  if (r < 0) {
            hdim->start = 0.0;
	  }
	  /* Get the attribute (step) from a minc file */
	  r = miget_attribute(volume, path, "step", MI_TYPE_DOUBLE, 1, &hdim->step);
	  if (r < 0) {
            hdim->step = 1.0;
	  }
	}
	/* Get the attribute (direction_cosines) from a minc file */
        r = miget_attribute(volume, path, "direction_cosines", MI_TYPE_DOUBLE, 3, 
                            hdim->direction_cosines);
        if (r < 0) {
            hdim->direction_cosines[MI2_X] = 0.0;
            hdim->direction_cosines[MI2_Y] = 0.0;
            hdim->direction_cosines[MI2_Z] = 0.0;
            if (!strcmp(dimname, MIxspace)) {
                hdim->direction_cosines[MI2_X] = 1.0;
            }
            else if (!strcmp(dimname, MIyspace)) {
                hdim->direction_cosines[MI2_Y] = 1.0;
            }
            else if (!strcmp(dimname, MIzspace)) {
                hdim->direction_cosines[MI2_Z] = 1.0;
            }
        }

        r = miget_attribute(volume, path, "units", MI_TYPE_STRING,
                            MI2_CHAR_LENGTH, temp);
        if (r < 0) {
            hdim->units = strdup("");
        }
        else {
            hdim->units = strdup(temp);
        }
            
    } H5E_END_TRY;
    /* Return the dimension handle */
    *hdim_ptr = hdim;
    hdim->volume_handle = volume;
    return (MI_NOERROR);
}


/** Opens an existing MINC volume for read-only access if mode argument is
    MI2_OPEN_READ, or read-write access if mode argument is MI2_OPEN_RDWR.
    \ingroup mi2Vol
 */
int
miopen_volume(const char *filename, int mode, mihandle_t *volume)
{
    hid_t file_id;
    hid_t dset_id;
    hid_t space_id;
    mihandle_t handle;
    int hdf_mode;
    char dimorder[MI2_CHAR_LENGTH];
    int i,r;
    char *p1, *p2;
    H5T_class_t class;
    size_t nbytes;
    int is_signed;

    /* Initialization. 
       For the actual body of this function look at m2utils.c
     */
    miinit();
     /* Convert the specified mode to hdf mode */
    if (mode == MI2_OPEN_READ) {
        hdf_mode = H5F_ACC_RDONLY;
    }
    else if (mode == MI2_OPEN_RDWR) {
        hdf_mode = H5F_ACC_RDWR;
    }
    else {
        return (MI_ERROR);
    }
    /* Open the hdf file using the given filename and mode */
    file_id = hdf_open(filename, hdf_mode);
    if (file_id < 0) {
	return (MI_ERROR);
    }
     /* Allocate space for the volume handle */
    handle = mialloc_volume_handle();
    if (handle == NULL) {
      return (MI_ERROR);
    }
    /* Set some varibales associated with the volume handle */
    handle->hdf_id = file_id;
    handle->mode = mode;

    /* Get the volume class.
     */
    _miget_volume_class(handle, &handle->volume_class);

    /* GET THE DIMENSION COUNT
     */
    handle->number_of_dims = _miget_file_dimension_count(file_id);
    
    /* READ EACH OF THE DIMENSIONS 
     */
    handle->dim_handles = (midimhandle_t *)malloc(handle->number_of_dims *
                                                  sizeof(midimhandle_t));
    /* Get the attribute (dimorder) from the image dataset */
    r =  miget_attribute(handle, "/minc-2.0/image/0/image", "dimorder", 
                    MI_TYPE_STRING, sizeof(dimorder), dimorder);
    
    if ( r < 0) {
      return (MI_ERROR);
    }
    p1 = dimorder;
    /* Break the ordered, comma-separated list of dimension names 
       to get each individual dimension name */
    for (i = 0; i < handle->number_of_dims; i++) {
        p2 = strchr(p1, ',');
        if (p2 != NULL) {
            *p2 = '\0';
        }
	/* Get dimension variable attributes for each dimension */
        _miget_file_dimension(handle, p1, &handle->dim_handles[i]);
        p1 = p2 + 1;
    }

    miset_volume_world_indices(handle);
    
    /* SEE IF SLICE SCALING IS ENABLED
     */
    handle->has_slice_scaling = FALSE;
    /* hdf5 macro can temporarily disable the automatic error printing */
    H5E_BEGIN_TRY {
         /* Open the dataset image-max at the specified path*/
        dset_id = H5Dopen1(file_id, "/minc-2.0/image/0/image-max");
    } H5E_END_TRY;
    if (dset_id >= 0) {
        /* Get the Id of the copy of the dataspace of the dataset */
	space_id = H5Dget_space(dset_id);
	if (space_id >= 0) {
	    /* If the dimensionality of the image-max variable is one or
	     * greater, we consider this volume to have slice-scaling enabled.
	     */
	    if (H5Sget_simple_extent_ndims(space_id) >= 1) {
		handle->has_slice_scaling = TRUE;
	    }
	    H5Sclose(space_id);	/* Close the dataspace handle */
	}
	H5Dclose(dset_id);	/* Close the dataset handle */
    }

    if (!handle->has_slice_scaling) {
        /* Read the minimum scalar of the given type at the specified path */
        miget_scalar(handle->hdf_id, H5T_NATIVE_DOUBLE,
                     "/minc-2.0/image/0/image-min", &handle->scale_min);
	/* Read the maximum scalar of the given type at the specified path */
        miget_scalar(handle->hdf_id, H5T_NATIVE_DOUBLE,
                     "/minc-2.0/image/0/image-max", &handle->scale_max);
    }
    
    /* Read the current voxel-to-world transform */
    miget_voxel_to_world(handle, handle->v2w_transform);

    /* Calculate the inverse transform */
    miinvert_transform(handle->v2w_transform, handle->w2v_transform);

    /* Open the image dataset */
    handle->image_id = H5Dopen1(file_id, "/minc-2.0/image/0/image");
    if (handle->image_id < 0) {
	return (MI_ERROR);
    }
    /* Get the Id for the copy of the datatype for the dataset */
    handle->ftype_id = H5Dget_type(handle->image_id);
    if (handle->ftype_id < 0) {
	return (MI_ERROR);
    }

    switch (H5Tget_class(handle->ftype_id)) {
    case H5T_INTEGER:
    case H5T_FLOAT:
        handle->mtype_id = H5Tget_native_type(handle->ftype_id, 
                                              H5T_DIR_ASCEND);
        break;

    case H5T_COMPOUND:
        handle->mtype_id = H5Tcreate(H5T_COMPOUND, 
                                     H5Tget_size(handle->ftype_id));
        for (i = 0; i < H5Tget_nmembers(handle->ftype_id); i++) {
            hid_t tmp_id = H5Tget_member_type(handle->ftype_id, i);
            size_t tmp_off = H5Tget_member_offset(handle->ftype_id, i);
            char *tmp_nm = H5Tget_member_name(handle->ftype_id, i);
            hid_t tmp2_id = H5Tget_native_type(tmp_id, H5T_DIR_ASCEND);
            H5Tinsert(handle->mtype_id, tmp_nm, tmp_off, tmp2_id);

            free(tmp_nm);
            H5Tclose(tmp_id);
            H5Tclose(tmp2_id);
        }
        break;

    case H5T_ENUM:
        handle->mtype_id = H5Tcopy(handle->ftype_id);
	miinit_enum(handle->ftype_id);
	miinit_enum(handle->mtype_id);
        /* Set native order ---> is not allowed after order is set */
        //H5Tset_order(handle->mtype_id, H5Tget_order(H5T_NATIVE_INT));
        break;

    default:
        return (MI_ERROR);
    }

    /* hdf5 macro can temporarily disable the automatic error printing */
    H5E_BEGIN_TRY {
        /* Open both image-min and image-max datasets */
        handle->imax_id = H5Dopen1(file_id, "/minc-2.0/image/0/image-max");
        handle->imin_id = H5Dopen1(file_id, "/minc-2.0/image/0/image-min");
    } H5E_END_TRY;

    /* Convert the type to a MINC type.
     */
    /* Get the class Id for the datatype */
    class = H5Tget_class(handle->ftype_id);
    /* Get the size of the datatype */
    nbytes = H5Tget_size(handle->ftype_id);

    switch (class) {
    case H5T_INTEGER:
    case H5T_ENUM:              /* label images */
	is_signed = (H5Tget_sign(handle->ftype_id) == H5T_SGN_2);

	switch (nbytes) {
	case 1:
	    handle->volume_type = (is_signed ? MI_TYPE_BYTE : MI_TYPE_UBYTE);
	    break;
	case 2:
	    handle->volume_type = (is_signed ? MI_TYPE_SHORT : MI_TYPE_USHORT);
	    break;
	case 4:
	    handle->volume_type = (is_signed ? MI_TYPE_INT : MI_TYPE_UINT);
	    break;
	default:
            return (MI_ERROR);
	}
	break;
    case H5T_FLOAT:
	handle->volume_type = (nbytes == 4) ? MI_TYPE_FLOAT : MI_TYPE_DOUBLE;
	break;
    case H5T_STRING:
	handle->volume_type = MI_TYPE_STRING;
	break;
    case H5T_ARRAY:
	/* TODO: handle this case for uniform records (arrays)? */
	break;
    case H5T_COMPOUND:
	/* TODO: handle this case for non-uniform records? */
	break;
    default:
        return (MI_ERROR);
    }

    /* Read the current settings for valid-range */
    miread_valid_range(handle, &handle->valid_max, &handle->valid_min);

    *volume = handle;
    return (MI_NOERROR);
}

/** Writes any changes associated with the volume to disk.
    \ingroup mi2Vol
 */
int
miflush_volume(mihandle_t volume)
{
    if ((volume->mode & MI2_OPEN_RDWR) != 0) {
        H5Fflush(volume->hdf_id, H5F_SCOPE_GLOBAL);
        misave_valid_range(volume);
    }
    return (MI_NOERROR);
}

/** Close an existing MINC volume. If the volume was newly created,
 *  all changes will be written to disk. In all cases this function closes
 *  the open volume and frees memory associated with the volume handle.
 *  \ingroup mi2Vol
 */
int 
miclose_volume(mihandle_t volume)
{
    if (volume == NULL) {
        return (MI_ERROR);
    }

    if (volume->is_dirty) {
        minc_update_thumbnails(volume);
        volume->is_dirty = FALSE;
    }

    miflush_volume(volume);

    if (volume->image_id > 0) {
        H5Dclose(volume->image_id);
    }
    if (volume->imax_id > 0) {
        H5Dclose(volume->imax_id);
    }
    if (volume->imin_id > 0) {
        H5Dclose(volume->imin_id);
    }
    if (volume->ftype_id > 0) {
        H5Tclose(volume->ftype_id);
    }
    if (volume->mtype_id > 0) {
        H5Tclose(volume->mtype_id);
    }
    if (volume->plist_id > 0) {
        H5Pclose(volume->plist_id);
    }
    //if (H5Fclose(volume->hdf_id) < 0) {
    if (hdf_close(volume->hdf_id) < 0) {
      return (MI_ERROR);
    }
    if (volume->dim_handles != NULL) {
        free(volume->dim_handles);
    }
    if (volume->dim_indices != NULL) {
        free(volume->dim_indices);
    }
    if (volume->create_props != NULL) {
      mifree_volume_props(volume->create_props);
    }
    free(volume);

    return (MI_NOERROR);
}



/* Internal functions
 */
static void
miinit_default_range(mitype_t mitype, double *valid_max, double *valid_min)
{
    switch (mitype) {
    case MI_TYPE_BYTE:
        *valid_min = CHAR_MIN;
        *valid_max = CHAR_MAX;
        break;
    case MI_TYPE_SHORT:
        *valid_min = SHRT_MIN;
        *valid_max = SHRT_MAX;
        break;
    case MI_TYPE_INT:
        *valid_min = INT_MIN;
        *valid_max = INT_MAX;
        break;
    case MI_TYPE_UBYTE:
        *valid_min = 0;
        *valid_max = UCHAR_MAX;
        break;
    case MI_TYPE_USHORT:
        *valid_min = 0;
        *valid_max = USHRT_MAX;
        break;
    case MI_TYPE_UINT:
        *valid_min = 0;
        *valid_max = UINT_MAX;
        break;
    case MI_TYPE_FLOAT:
        *valid_min = -FLT_MAX;
        *valid_max = FLT_MAX;
        break;
    case MI_TYPE_DOUBLE:
        *valid_min = -DBL_MAX;
        *valid_max = DBL_MAX;
        break;
    default:
        *valid_min = 0;
        *valid_max = 1;
        break;
    }
}

static void
miread_valid_range(mihandle_t volume, double *valid_max, double *valid_min)
{
    int r;
    double range[2];

    H5E_BEGIN_TRY {
        r = miget_attribute(volume, "/minc-2.0/image/0/image", "valid_range",
                            MI_TYPE_DOUBLE, 2, range);
    } H5E_END_TRY;
    if (r == MI_NOERROR) {
        if (range[0] < range[1]) {
            *valid_min = range[0];
            *valid_max = range[1];
        }
        else {
            *valid_min = range[1];
            *valid_max = range[0];
        }
    }
    else {
        /* Didn't find the attribute, so assign default values. */
        miinit_default_range(volume->volume_type, valid_max, valid_min);
    }
}

/** \internal
 * This function saves the current valid range set for a MINC file.
 */
void
misave_valid_range(mihandle_t volume)
{
    double range[2];
    range[0] = volume->valid_min;
    range[1] = volume->valid_max;
    miset_attribute(volume, "/minc-2.0/image/0/image", "valid_range",
                    MI_TYPE_DOUBLE, 2, range);
}