/** \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 = H5Dcreate(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 = H5Dcreate(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 = H5Dcreate(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.
*/
switch (volume_type) {
case MI_TYPE_BYTE:
case MI_TYPE_SHORT:
case MI_TYPE_INT:
case MI_TYPE_UBYTE:
case MI_TYPE_USHORT:
case MI_TYPE_UINT:
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_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 = H5Gopen(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 = H5Dcreate(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;
name = malloc(size);
strcpy(name, dimensions[i]->name);
strcat(name, "-width");
/* Create dataset dimension_name-width */
dataset_width = H5Dcreate(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 = H5Dopen(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, "xspace")) {
hdim->direction_cosines[MI2_X] = 1.0;
}
else if (!strcmp(dimname, "yspace")) {
hdim->direction_cosines[MI2_Y] = 1.0;
}
else if (!strcmp(dimname, "zspace")) {
hdim->direction_cosines[MI2_Z] = 1.0;
}
}
} 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 = H5Dopen(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 = H5Dopen(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);
/* Set native order */
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 = H5Dopen(file_id, "/minc-2.0/image/0/image-max");
handle->imin_id = H5Dopen(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:
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) {
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);
}
syntax highlighted by Code2HTML, v. 0.9.1