/** \file m2util.c
* \brief MINC 2.0 Private utility functions
* \author Leila Baghdadi, Bert Vincent
*
************************************************************************/
#include <stdlib.h>
#include <math.h>
#include <hdf5.h>
#include <limits.h>
#include <float.h>
#include "minc2.h"
#include "minc2_private.h"
/* Uggh!!! The HDF5 team changed the definition of the H5Tconvert(),
* H5Tregister(), and H5T_conv_t functions, and the result is that we
* have to special-case these types. I am bummed.
*/
#if (H5_VERS_MAJOR > 1) || (H5_VERS_MINOR > 6) || (H5_VERS_RELEASE > 2)
#define H5_NELEMENTS_T size_t
#else
#define H5_NELEMENTS_T hsize_t
#endif
/* They also redefined the type of the 4th argument to H5Sselect_elements.
* This is harmless as long as sizeof(hssize_t) == sizeof(hsize_t).
*/
#if (H5_VERS_MAJOR > 1) || (H5_VERS_MINOR > 6) || (H5_VERS_RELEASE > 3)
#define H5_START_T hsize_t
#else
#define H5_START_T hssize_t
#endif
/*! Convert a MINC 2 datatype into a HDF5 datatype. Actually returns a copy
* of the datatype, so the returned value must be explicitly freed with a
* call to H5Tclose().
* \param mitype The MINC 2 data type to convert
* \param is_native Convert to native type if TRUE
*/
hid_t
mitype_to_hdftype(mitype_t mitype, int is_native)
{
hid_t type_id;
if (is_native) {
/* Native types are in the byte-ordering of the native system.
* They are appropriate for the "in-memory" types for data.
*/
switch (mitype) {
case MI_TYPE_BYTE:
type_id = H5Tcopy(H5T_NATIVE_SCHAR);
break;
case MI_TYPE_SHORT:
type_id = H5Tcopy(H5T_NATIVE_SHORT);
break;
case MI_TYPE_INT:
type_id = H5Tcopy(H5T_NATIVE_INT);
break;
case MI_TYPE_FLOAT:
type_id = H5Tcopy(H5T_NATIVE_FLOAT);
break;
case MI_TYPE_DOUBLE:
type_id = H5Tcopy(H5T_NATIVE_DOUBLE);
break;
case MI_TYPE_UBYTE:
type_id = H5Tcopy(H5T_NATIVE_UCHAR);
break;
case MI_TYPE_USHORT:
type_id = H5Tcopy(H5T_NATIVE_USHORT);
break;
case MI_TYPE_UINT:
type_id = H5Tcopy(H5T_NATIVE_UINT);
break;
case MI_TYPE_SCOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 4);
H5Tinsert(type_id, "real", 0, H5T_NATIVE_SHORT);
H5Tinsert(type_id, "imag", 2, H5T_NATIVE_SHORT);
break;
case MI_TYPE_ICOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 8);
H5Tinsert(type_id, "real", 0, H5T_NATIVE_INT);
H5Tinsert(type_id, "imag", 4, H5T_NATIVE_INT);
break;
case MI_TYPE_FCOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 8);
H5Tinsert(type_id, "real", 0, H5T_NATIVE_FLOAT);
H5Tinsert(type_id, "imag", 4, H5T_NATIVE_FLOAT);
break;
case MI_TYPE_DCOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 16);
H5Tinsert(type_id, "real", 0, H5T_NATIVE_DOUBLE);
H5Tinsert(type_id, "imag", 8, H5T_NATIVE_DOUBLE);
break;
default:
type_id = H5Tcopy(mitype); /* It is a standard HDF type handle? */
break;
}
}
else {
/* The non-native types are standardized to be in
* "little-endian" form. That's an arbitrary decision which
* could certainly be debated.
*/
switch (mitype) {
case MI_TYPE_BYTE:
type_id = H5Tcopy(H5T_STD_I8LE);
break;
case MI_TYPE_SHORT:
type_id = H5Tcopy(H5T_STD_I16LE);
break;
case MI_TYPE_INT:
type_id = H5Tcopy(H5T_STD_I32LE);
break;
case MI_TYPE_FLOAT:
type_id = H5Tcopy(H5T_IEEE_F32LE);
break;
case MI_TYPE_DOUBLE:
type_id = H5Tcopy(H5T_IEEE_F64LE);
break;
case MI_TYPE_UBYTE:
type_id = H5Tcopy(H5T_STD_U8LE);
break;
case MI_TYPE_USHORT:
type_id = H5Tcopy(H5T_STD_U16LE);
break;
case MI_TYPE_UINT:
type_id = H5Tcopy(H5T_STD_U32LE);
break;
case MI_TYPE_SCOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 4);
H5Tinsert(type_id, "real", 0, H5T_STD_I16LE);
H5Tinsert(type_id, "imag", 2, H5T_STD_I16LE);
break;
case MI_TYPE_ICOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 8);
H5Tinsert(type_id, "real", 0, H5T_STD_I32LE);
H5Tinsert(type_id, "imag", 4, H5T_STD_I32LE);
break;
case MI_TYPE_FCOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 8);
H5Tinsert(type_id, "real", 0, H5T_IEEE_F32LE);
H5Tinsert(type_id, "imag", 4, H5T_IEEE_F32LE);
break;
case MI_TYPE_DCOMPLEX:
type_id = H5Tcreate(H5T_COMPOUND, 16);
H5Tinsert(type_id, "real", 0, H5T_IEEE_F64LE);
H5Tinsert(type_id, "imag", 8, H5T_IEEE_F64LE);
break;
default:
type_id = H5Tcopy(mitype); /* It is a standard HDF type handle? */
break;
}
}
return (type_id);
}
/*! Convert a MINC 2 datatype into a NetCDF datatype.
* \param mitype The MINC 2 data type to convert
* \param is_signed Set to TRUE if the data type is signed, FALSE if unsigned.
*/
int
mitype_to_nctype(mitype_t mitype, int *is_signed)
{
int nctype;
*is_signed = 1; /* Assume signed by default. */
switch (mitype) {
case MI_TYPE_BYTE:
nctype = NC_BYTE;
break;
case MI_TYPE_SHORT:
nctype = NC_SHORT;
break;
case MI_TYPE_INT:
nctype = NC_INT;
break;
case MI_TYPE_FLOAT:
nctype = NC_FLOAT;
break;
case MI_TYPE_DOUBLE:
nctype = NC_DOUBLE;
break;
case MI_TYPE_UBYTE:
nctype = NC_BYTE;
*is_signed = 0;
break;
case MI_TYPE_USHORT:
nctype = NC_SHORT;
*is_signed = 0;
break;
case MI_TYPE_UINT:
nctype = NC_INT;
*is_signed = 1;
break;
default:
nctype = -1; /* ERROR!! */
break;
}
return (nctype);
}
/*! Return the group or dataset ID of the last item in a "path",
* specified like a UNIX pathname /black/white/red etc.
* \param file_id The HDF5 handle of the file (or group) at which to start
* the search.
* \param path A string consisting of a slash-separated list of
* HDF5 groupnames
*/
hid_t
midescend_path(hid_t file_id, const char *path)
{
hid_t tmp_id;
/* Put H5E_BEGIN_TRY/H5E_END_TRY around this to avoid the overzealous
* automatic error reporting of HDF5.
*/
H5E_BEGIN_TRY {
tmp_id = H5Dopen(file_id, path);
/* If the dataset open fails, try opening the object as a group.
*/
if (tmp_id < 0) {
tmp_id = H5Gopen(file_id, path);
}
} H5E_END_TRY;
return (tmp_id);
}
/** Get the number of voxels in the file - this is the total number,
* not just spatial dimensions
*/
mi_i64_t
miget_voxel_count(int mincid)
{
int imgid;
int dim[MI2_MAX_VAR_DIMS];
int idim;
int ndims;
mi_i64_t nvoxels;
long length;
/* Get the dimension ids for the image variable
*/
imgid = ncvarid(mincid, MIimage);
(void)ncvarinq(mincid, imgid, NULL, NULL, &ndims, dim, NULL);
/* Loop over them to get the total number of voxels
*/
nvoxels = 1;
for (idim = 0; idim < ndims; idim++) {
(void)ncdiminq(mincid, dim[idim], NULL, &length);
nvoxels *= length;
}
return (nvoxels);
}
int
miset_attr_at_loc(hid_t hdf_loc, const char *name, mitype_t data_type,
int length, const void *values)
{
hid_t ftyp_id;
hid_t mtyp_id;
hid_t spc_id;
hid_t hdf_attr;
hsize_t hdf_len;
H5E_BEGIN_TRY {
/* Delete attribute if it already exists. */
H5Adelete(hdf_loc, name);
} H5E_END_TRY;
switch (data_type) {
case MI_TYPE_INT:
ftyp_id = H5Tcopy(H5T_STD_I32LE);
mtyp_id = H5Tcopy(H5T_NATIVE_INT);
break;
case MI_TYPE_FLOAT:
ftyp_id = H5Tcopy(H5T_IEEE_F32LE);
mtyp_id = H5Tcopy(H5T_NATIVE_FLOAT);
break;
case MI_TYPE_DOUBLE:
ftyp_id = H5Tcopy(H5T_IEEE_F64LE);
mtyp_id = H5Tcopy(H5T_NATIVE_DOUBLE);
break;
case MI_TYPE_STRING:
ftyp_id = H5Tcopy(H5T_C_S1);
H5Tset_size(ftyp_id, length);
mtyp_id = H5Tcopy(ftyp_id);
length = 1; /* Apparent length is one. */
break;
default:
return (MI_ERROR);
}
if (length == 1) {
spc_id = H5Screate(H5S_SCALAR);
}
else {
hdf_len = (hsize_t) length;
spc_id = H5Screate_simple(1, &hdf_len, NULL);
}
if (spc_id < 0) {
return (MI_ERROR);
}
hdf_attr = H5Acreate(hdf_loc, name, ftyp_id, spc_id, H5P_DEFAULT);
if (hdf_attr < 0) {
return (MI_ERROR);
}
H5Awrite(hdf_attr, mtyp_id, values);
H5Aclose(hdf_attr);
H5Tclose(ftyp_id);
H5Tclose(mtyp_id);
H5Sclose(spc_id);
return (MI_NOERROR);
}
/** Set an attribute from a minc file */
int
miset_attribute(mihandle_t volume, const char *path, const char *name,
mitype_t data_type, int length, const void *values)
{
hid_t hdf_file;
hid_t hdf_loc;
/* Get a handle to the actual HDF file
*/
hdf_file = volume->hdf_id;
if (hdf_file < 0) {
return (MI_ERROR);
}
/* Search through the path, descending into each group encountered.
*/
hdf_loc = midescend_path(hdf_file, path);
if (hdf_loc < 0) {
return (MI_ERROR);
}
miset_attr_at_loc(hdf_loc, name, data_type, length, values);
/* The hdf_loc identifier could be a group or a dataset.
*/
if (H5Iget_type(hdf_loc) == H5I_GROUP) {
H5Gclose(hdf_loc);
}
else {
H5Dclose(hdf_loc);
}
return (MI_NOERROR);
}
/** Get a double attribute from a minc file */
int
miget_attribute(mihandle_t volume, const char *path, const char *name,
mitype_t data_type, int length, void *values)
{
hid_t hdf_file;
hid_t hdf_loc;
hid_t mtyp_id; /* Parameter type */
hid_t spc_id;
hid_t hdf_attr;
int status;
/* Get a handle to the actual HDF file
*/
hdf_file = volume->hdf_id;
if (hdf_file < 0) {
return (MI_ERROR);
}
/* Find the group or dataset referenced by the path.
*/
hdf_loc = midescend_path(hdf_file, path);
if (hdf_loc < 0) {
return (MI_ERROR);
}
H5E_BEGIN_TRY {
hdf_attr = H5Aopen_name(hdf_loc, name);
} H5E_END_TRY;
if (hdf_attr < 0) {
return (MI_ERROR);
}
switch (data_type) {
case MI_TYPE_INT:
mtyp_id = H5Tcopy(H5T_NATIVE_INT);
break;
case MI_TYPE_UINT:
mtyp_id = H5Tcopy(H5T_NATIVE_UINT);
break;
case MI_TYPE_FLOAT:
mtyp_id = H5Tcopy(H5T_NATIVE_FLOAT);
break;
case MI_TYPE_DOUBLE:
mtyp_id = H5Tcopy(H5T_NATIVE_DOUBLE);
break;
case MI_TYPE_STRING:
mtyp_id = H5Tcopy(H5T_C_S1);
H5Tset_size(mtyp_id, length);
break;
default:
return (MI_ERROR);
}
spc_id = H5Aget_space(hdf_attr);
if (spc_id < 0) {
return (MI_ERROR);
}
/* If we're retrieving a vector, make certain the length passed into this
* function is sufficient.
*/
if (H5Sget_simple_extent_ndims(spc_id) == 1) {
hsize_t hdf_dims[1];
H5Sget_simple_extent_dims(spc_id, hdf_dims, NULL);
if (length < hdf_dims[0]) {
return (MI_ERROR);
}
}
status = H5Aread(hdf_attr, mtyp_id, values);
if (status < 0) {
return (MI_ERROR);
}
/* Be certain that the string is null-terminated.
*/
if (data_type == MI_TYPE_STRING) {
hid_t atype; /* Attribute type */
int alength;
atype = H5Aget_type(hdf_attr);
alength = H5Tget_size(atype);
((char *)values)[alength] = '\0';
H5Tclose(atype);
}
H5Aclose(hdf_attr);
H5Tclose(mtyp_id);
H5Sclose(spc_id);
/* The hdf_loc identifier could be a group or a dataset.
*/
if (H5Iget_type(hdf_loc) == H5I_GROUP) {
H5Gclose(hdf_loc);
}
else {
H5Dclose(hdf_loc);
}
return (MI_NOERROR);
}
/* Get the mapping from spatial dimension - x, y, z - to file dimensions
* and vice-versa.
*/
void
mifind_spatial_dims(int mincid, int space_to_dim[], int dim_to_space[])
{
int imgid, dim[MI2_MAX_VAR_DIMS];
int idim, ndims, world_index;
char dimname[MAX_NC_NAME];
/* Set default values */
for (idim = 0; idim < 3; idim++)
space_to_dim[idim] = -1;
for (idim = 0; idim < MI2_MAX_VAR_DIMS; idim++)
dim_to_space[idim] = -1;
/* Get the dimension ids for the image variable
*/
imgid = ncvarid(mincid, MIimage);
(void)ncvarinq(mincid, imgid, NULL, NULL, &ndims, dim, NULL);
/* Loop over them to find the spatial ones
*/
for (idim = 0; idim < ndims; idim++) {
/* Get the name and check that this is a spatial dimension
*/
(void)ncdiminq(mincid, dim[idim], dimname, NULL);
if (!strcmp(dimname, MIxspace)) {
world_index = MI2_X;
}
else if (!strcmp(dimname, MIyspace)) {
world_index = MI2_Y;
}
else if (!strcmp(dimname, MIzspace)) {
world_index = MI2_Z;
}
else {
continue;
}
space_to_dim[world_index] = idim;
dim_to_space[idim] = world_index;
}
}
/** Get the voxel to world transform (for column vectors) */
void
miget_voxel_to_world(mihandle_t volume, mi_lin_xfm_t voxel_to_world)
{
int i;
int j;
int k;
double dircos[MI2_3D];
double step;
double start;
/* Zero the matrix */
for (i = 0; i < MI2_LIN_XFM_SIZE; i++) {
for (j = 0; j < MI2_LIN_XFM_SIZE; j++) {
voxel_to_world[i][j] = 0.0;
}
voxel_to_world[i][i] = 1.0;
}
for (j = 0; j < volume->number_of_dims; j++) {
midimhandle_t hdim = volume->dim_handles[j];
if (hdim->class == MI_DIMCLASS_SPATIAL ||
hdim->class == MI_DIMCLASS_SFREQUENCY) {
k = hdim->world_index;
}
else {
continue;
}
start = hdim->start;
step = hdim->step;
dircos[0] = hdim->direction_cosines[0];
dircos[1] = hdim->direction_cosines[1];
dircos[2] = hdim->direction_cosines[2];
minormalize_vector(dircos);
/* Put them in the matrix */
for (i = 0; i < MI2_3D; i++) {
voxel_to_world[i][k] = step * dircos[i];
voxel_to_world[i][MI2_3D] += start * dircos[i];
}
}
}
/** Normalize a 3 element vector */
void
minormalize_vector(double vector[MI2_3D])
{
int i;
double magnitude;
magnitude = 0.0;
for (i = 0; i < MI2_3D; i++) {
magnitude += (vector[i] * vector[i]);
}
magnitude = sqrt(magnitude);
if (magnitude > 0.0) {
for (i = 0; i < MI2_3D; i++) {
vector[i] /= magnitude;
}
}
}
/** Transforms a coordinate through a linear transform */
void
mitransform_coord(double out_coord[],
mi_lin_xfm_t transform,
const double in_coord[])
{
int i, j;
double in_homogeneous[MI2_3D + 1];
double out_homogeneous[MI2_3D + 1];
for (i = 0; i < MI2_3D; i++) {
in_homogeneous[i] = in_coord[i];
}
in_homogeneous[MI2_3D] = 1.0;
for (i = 0; i < MI2_3D + 1; i++) {
out_homogeneous[i] = 0.0;
for (j = 0; j < MI2_3D + 1; j++) {
out_homogeneous[i] += transform[i][j] * in_homogeneous[j];
}
}
#if 0
printf("W = %f\n", out_homogeneous[3]);
#endif /* 0 */
for (i = 0; i < MI2_3D; i++) {
out_coord[i] = out_homogeneous[i];
}
}
/** For conversions from double to integer, rounding may be performed
* by setting this variable to non-zero.
* However, at present, no API is available to control this.
*/
static int rounding_enabled = FALSE;
static void miswap8(unsigned char *tmp_ptr)
{
unsigned char x;
x = tmp_ptr[0];
tmp_ptr[0] = tmp_ptr[7];
tmp_ptr[7] = x;
x = tmp_ptr[1];
tmp_ptr[1] = tmp_ptr[6];
tmp_ptr[6] = x;
x = tmp_ptr[2];
tmp_ptr[2] = tmp_ptr[5];
tmp_ptr[5] = x;
x = tmp_ptr[3];
tmp_ptr[3] = tmp_ptr[4];
tmp_ptr[4] = x;
}
static void miswap4(unsigned char *tmp_ptr)
{
unsigned char x;
x = tmp_ptr[0];
tmp_ptr[0] = tmp_ptr[3];
tmp_ptr[3] = x;
x = tmp_ptr[1];
tmp_ptr[1] = tmp_ptr[2];
tmp_ptr[2] = x;
}
static void miswap2(unsigned char *tmp_ptr)
{
unsigned char x;
x = tmp_ptr[0];
tmp_ptr[0] = tmp_ptr[1];
tmp_ptr[1] = x;
}
/** Generic HDF5 integer-to-double converter.
*/
static herr_t
mi2_int_to_dbl(hid_t src_id,
hid_t dst_id,
H5T_cdata_t *cdata,
H5_NELEMENTS_T nelements,
size_t buf_stride,
size_t bkg_stride,
void *buf_ptr,
void *bkg_ptr,
hid_t dset_xfer_plist)
{
unsigned char *dst_ptr;
unsigned char *src_ptr;
int src_nb;
int dst_nb;
H5T_sign_t src_sg;
double t;
int dst_cnt;
int src_cnt;
int src_swap;
int dst_swap;
switch (cdata->command) {
case H5T_CONV_INIT:
cdata->need_bkg = H5T_BKG_NO;
src_nb = H5Tget_size(src_id);
if (src_nb != 1 && src_nb != 2 && src_nb != 4) {
return (-1);
}
dst_nb = H5Tget_size(dst_id);
if (dst_nb != 8) {
return (-1);
}
break;
case H5T_CONV_CONV:
src_nb = H5Tget_size(src_id);
src_sg = H5Tget_sign(src_id);
dst_nb = H5Tget_size(dst_id);
if (buf_stride == 0) {
dst_cnt = dst_nb;
src_cnt = src_nb;
}
else {
dst_cnt = buf_stride;
src_cnt = buf_stride;
}
/* Convert starting from "far side" of buffer... (Hope this works!)
*/
dst_ptr = ((unsigned char *) buf_ptr) + ((nelements - 1) * dst_nb);
src_ptr = ((unsigned char *) buf_ptr) + ((nelements - 1) * src_nb);
if (H5Tget_order(H5T_NATIVE_INT) != H5Tget_order(src_id)) {
src_swap = 1;
}
else {
src_swap = 0;
}
if (H5Tget_order(H5T_NATIVE_DOUBLE) != H5Tget_order(dst_id)) {
dst_swap = 1;
}
else {
dst_swap = 0;
}
if (src_sg == H5T_SGN_2) {
switch (src_nb) {
case 4:
while (nelements-- > 0) {
if (src_swap) {
miswap4(src_ptr);
}
t = *(int *)src_ptr;
*(double *)dst_ptr = t;
if (dst_swap) {
miswap8(dst_ptr);
}
src_ptr -= src_cnt;
dst_ptr -= dst_cnt;
}
break;
case 2:
while (nelements-- > 0) {
if (src_swap) {
miswap2(src_ptr);
}
t = *(short *)src_ptr;
*(double *)dst_ptr = t;
if (dst_swap) {
miswap8(dst_ptr);
}
src_ptr -= src_cnt;
dst_ptr -= dst_cnt;
}
break;
case 1:
while (nelements-- > 0) {
t = *(char *)src_ptr;
*(double *)dst_ptr = t;
if (dst_swap) {
miswap8(dst_ptr);
}
src_ptr -= src_cnt;
dst_ptr -= dst_cnt;
}
break;
default:
/* Should not happen! */
break;
}
}
else {
switch (src_nb) {
case 4:
while (nelements-- > 0) {
if (src_swap) {
miswap4(src_ptr);
}
t = *(unsigned int *)src_ptr;
*(double *)dst_ptr = t;
if (dst_swap) {
miswap8(dst_ptr);
}
src_ptr -= src_cnt;
dst_ptr -= dst_cnt;
}
break;
case 2:
while (nelements-- > 0) {
if (src_swap) {
miswap2(src_ptr);
}
t = *(unsigned short *)src_ptr;
*(double *)dst_ptr = t;
if (dst_swap) {
miswap8(dst_ptr);
}
src_ptr -= src_cnt;
dst_ptr -= dst_cnt;
}
break;
case 1:
while (nelements-- > 0) {
t = *(unsigned char *)src_ptr;
*(double *)dst_ptr = t;
if (dst_swap) {
miswap8(dst_ptr);
}
src_ptr -= src_cnt;
dst_ptr -= dst_cnt;
}
break;
default:
/* Should not happen! */
break;
}
}
break;
case H5T_CONV_FREE:
break;
default:
/* Unknown command */
return (-1);
}
return (0);
}
/** Generic HDF5 double-to-integer converter.
*/
static herr_t
mi2_dbl_to_int(hid_t src_id,
hid_t dst_id,
H5T_cdata_t *cdata,
H5_NELEMENTS_T nelements,
size_t buf_stride,
size_t bkg_stride,
void *buf_ptr,
void *bkg_ptr,
hid_t dset_xfer_plist)
{
unsigned char *dst_ptr;
unsigned char *src_ptr;
int src_nb;
int dst_nb;
H5T_sign_t dst_sg;
double t;
int dst_cnt;
int src_cnt;
int src_swap;
int dst_swap;
switch (cdata->command) {
case H5T_CONV_INIT:
cdata->need_bkg = H5T_BKG_NO;
/* Verify that we can handle this conversion.
*/
src_nb = H5Tget_size(src_id);
if (src_nb != 8) {
return (-1);
}
dst_nb = H5Tget_size(dst_id);
if (dst_nb != 4 && dst_nb != 2 && dst_nb != 1) {
return (-1);
}
break;
case H5T_CONV_CONV:
dst_nb = H5Tget_size(dst_id);
dst_sg = H5Tget_sign(dst_id);
src_nb = H5Tget_size(src_id);
dst_ptr = (unsigned char *) buf_ptr;
src_ptr = (unsigned char *) buf_ptr;
if (H5Tget_order(H5T_NATIVE_DOUBLE) != H5Tget_order(src_id)) {
src_swap = 1;
}
else {
src_swap = 0;
}
if (H5Tget_order(H5T_NATIVE_INT) != H5Tget_order(dst_id)) {
dst_swap = 1;
}
else {
dst_swap = 0;
}
/* The logic of HDF5 seems to be that if a stride is specified,
* both the source and destination pointers should advance by that
* amount. This seems wrong to me, but I've examined the HDF5 sources
* and that's what their own type converters do.
*/
if (buf_stride == 0) {
dst_cnt = dst_nb;
src_cnt = src_nb;
}
else {
dst_cnt = buf_stride;
src_cnt = buf_stride;
}
if (rounding_enabled) {
if (dst_sg == H5T_SGN_2) {
switch (dst_nb) {
case 4:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = rint(*(double *) src_ptr);
if (t > INT_MAX) {
t = INT_MAX;
}
else if (t < INT_MIN) {
t = INT_MIN;
}
*((int *)dst_ptr) = (int) t;
if (dst_swap) {
miswap4(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 2:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = rint(*(double *) src_ptr);
if (t > SHRT_MAX) {
t = SHRT_MAX;
}
else if (t < SHRT_MIN) {
t = SHRT_MIN;
}
*((short *)dst_ptr) = (short) t;
if (dst_swap) {
miswap2(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 1:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = rint(*(double *) src_ptr);
if (t > CHAR_MAX) {
t = CHAR_MAX;
}
else if (t < CHAR_MIN) {
t = CHAR_MIN;
}
*((char *)src_ptr) = (char) t;
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
default:
/* Can't handle this! */
break;
}
}
else {
switch (dst_nb) {
case 4:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = rint(*(double *)src_ptr);
if (t > UINT_MAX) {
t = UINT_MAX;
}
else if (t < 0) {
t = 0;
}
*((unsigned int *)dst_ptr) = (unsigned int) t;
if (dst_swap) {
miswap4(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 2:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = rint(*(double *)src_ptr);
if (t > USHRT_MAX) {
t = USHRT_MAX;
}
else if (t < 0) {
t = 0;
}
*((unsigned short *)dst_ptr) = (unsigned short) t;
if (dst_swap) {
miswap2(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 1:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = rint(*(double *)src_ptr);
if (t > UCHAR_MAX) {
t = UCHAR_MAX;
}
else if (t < 0) {
t = 0;
}
*((unsigned char *)dst_ptr) = (unsigned char) t;
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
default:
/* Can't handle any other values */
break;
}
}
}
else {
if (dst_sg == H5T_SGN_2) {
switch (dst_nb) {
case 4:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = (int)(*(double *) src_ptr);
if (t > INT_MAX) {
t = INT_MAX;
}
else if (t < INT_MIN) {
t = INT_MIN;
}
*((int *)dst_ptr) = (int) t;
if (dst_swap) {
miswap4(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 2:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = (int)(*(double *) src_ptr);
if (t > SHRT_MAX) {
t = SHRT_MAX;
}
else if (t < SHRT_MIN) {
t = SHRT_MIN;
}
*((short *)dst_ptr) = (short) t;
if (dst_swap) {
miswap4(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 1:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = (int)(*(double *) src_ptr);
if (t > CHAR_MAX) {
t = CHAR_MAX;
}
else if (t < CHAR_MIN) {
t = CHAR_MIN;
}
*((char *)src_ptr) = (char) t;
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
default:
/* Can't handle this! */
break;
}
}
else {
switch (dst_nb) {
case 4:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = (int)(*(double *)src_ptr);
if (t > UINT_MAX) {
t = UINT_MAX;
}
else if (t < 0) {
t = 0;
}
*((unsigned int *)dst_ptr) = (unsigned int) t;
if (dst_swap) {
miswap4(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 2:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = (int)(*(double *)src_ptr);
if (t > USHRT_MAX) {
t = USHRT_MAX;
}
else if (t < 0) {
t = 0;
}
*((unsigned short *)dst_ptr) = (unsigned short) t;
if (dst_swap) {
miswap2(dst_ptr);
}
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
case 1:
while (nelements-- > 0) {
if (src_swap) {
miswap8(src_ptr);
}
t = (int)(*(double *)src_ptr);
if (t > UCHAR_MAX) {
t = UCHAR_MAX;
}
else if (t < 0) {
t = 0;
}
*((unsigned char *)dst_ptr) = (unsigned char) t;
dst_ptr += dst_cnt;
src_ptr += src_cnt;
}
break;
default:
/* Can't handle any other values */
break;
}
}
}
break;
case H5T_CONV_FREE:
break;
default:
/* Unknown command */
return (-1);
}
return (0);
}
/** Initialize some critical pieces of the library. For now all this does
is install the double-to-integer and integer-to-double conversion functions.
*/
void
miinit(void)
{
H5Tregister(H5T_PERS_SOFT, "i2d", H5T_NATIVE_INT, H5T_NATIVE_DOUBLE,
mi2_int_to_dbl);
H5Tregister(H5T_PERS_SOFT, "d2i", H5T_NATIVE_DOUBLE, H5T_NATIVE_INT,
mi2_dbl_to_int);
}
/** HDF5 type conversion function for converting an arbitrary integer type to
* an arbitrary enumerated type. The beauty part of this is that it is
* not necessary to actually perform any real conversion!
*/
herr_t mi2_null_conv(hid_t src_id,
hid_t dst_id,
H5T_cdata_t *cdata,
H5_NELEMENTS_T nelements,
size_t buf_stride,
size_t bkg_stride,
void *buf_ptr,
void *bkg_ptr,
hid_t dset_xfer_plist)
{
switch (cdata->command) {
case H5T_CONV_INIT:
break;
case H5T_CONV_CONV:
break;
case H5T_CONV_FREE:
break;
default:
/* Unknown command */
return (-1);
}
return (0);
}
/**
This function should be called when a labeled volume is created or opened
in order to facilitate conversions from the integer to the enumerated type.
*/
void
miinit_enum(hid_t type_id)
{
H5Tregister(H5T_PERS_SOFT, "i2e", H5T_NATIVE_INT, type_id,
mi2_null_conv);
H5Tregister(H5T_PERS_SOFT, "e2i", type_id, H5T_NATIVE_INT,
mi2_null_conv);
H5Tregister(H5T_PERS_SOFT, "d2e", H5T_NATIVE_DOUBLE, type_id,
mi2_dbl_to_int);
H5Tregister(H5T_PERS_SOFT, "e2d", type_id, H5T_NATIVE_DOUBLE,
mi2_int_to_dbl);
}
int
minc_create_thumbnail(mihandle_t volume, int grp)
{
char path[MI2_MAX_PATH];
hid_t grp_id;
/* Don't handle negative or overly large numbers!
*/
if (grp <= 0 || grp > MI2_MAX_RESOLUTION_GROUP) {
return (MI_ERROR);
}
sprintf(path, "/minc-2.0/image/%d", grp);
grp_id = H5Gcreate(volume->hdf_id, path, 0);
if (grp_id < 0) {
return (MI_ERROR);
}
H5Gclose(grp_id);
return (MI_NOERROR);
}
/** Function to downsample a single slice of an image.
* \param in_ptr the 3D input slice, scale x isize[1] x isize[2]
* \param out_ptr the 2D output slice, osize[1] x osize[2]
*/
static void
midownsample_slice(double *in_ptr, double *out_ptr, hsize_t isize[],
hsize_t osize[], int scale)
{
int j, k;
int x, y, z;
double d;
hsize_t total;
total = scale * scale * scale;
/* These two loops iterate over all of the voxels in the 2D output
* image.
*/
for (j = 0; j < osize[1]; j++) {
for (k = 0; k < osize[2]; k++) {
/* The three inner loops iterate all scale^3
* voxels in the input image which will be averaged
* to form the output image.
*/
d = 0;
for (x = 0; x < scale; x++) {
for (y = 0; y < scale; y++) {
for (z = 0; z < scale; z++) {
int x1,y1,z1;
double t;
x1 = x;
y1 = y + (j * scale);
z1 = z + (k * scale);
t = in_ptr[((x1 * isize[1]) + y1) * isize[2] + z1];
d += t;
}
}
}
d /= total;
out_ptr[(j * osize[2]) + k] = d;
}
}
}
/** Convert the hyperslab from real to voxel values, calculating and
* returning the minimum and maximum real values for the slab. This
* could form the basis for a public function one day, but for now it
* is considered private.
*/
static void
miconvert_hyperslab_to_voxel(mihandle_t volume, hssize_t start[],
hsize_t count[], double *slab_ptr,
double *max_ptr, double *min_ptr)
{
/* This code is not intended to be a general hyperslab-to-voxel
* converter yet. That is why it is not public.
*/
double real_min, real_max; /* Minimum and maximum values */
hsize_t index;
hsize_t total;
double voxel_range, voxel_offset;
double real_range, real_offset;
int i;
double tmp_val;
real_min = DBL_MAX;
real_max = -DBL_MAX;
total = 1;
for (i = 0; i < volume->number_of_dims; i++) {
total *= count[i];
}
/* Find the global minimum and maximum for this hyperslab.
*/
for (index = 0; index < total; index++) {
tmp_val = slab_ptr[index];
if (tmp_val > real_max) {
real_max = tmp_val;
}
if (tmp_val < real_min) {
real_min = tmp_val;
}
}
voxel_range = volume->valid_max - volume->valid_min;
voxel_offset = volume->valid_min;
real_range = real_max - real_min;
real_offset = real_min;
for (index = 0; index < total; index++) {
tmp_val = slab_ptr[index];
tmp_val = (tmp_val - real_offset) / real_range;
tmp_val = (tmp_val * voxel_range) + voxel_offset;
slab_ptr[index] = rint(tmp_val);
}
if (min_ptr != NULL) {
*min_ptr = real_min;
}
if (max_ptr != NULL) {
*max_ptr = real_max;
}
}
/** Convert the hyperslab from voxel to real values. This could form
* the basis for a public function one day, but for now it is
* considered private.
*/
static void
miconvert_hyperslab_to_real(mihandle_t volume, hssize_t start[],
hsize_t count[], double *slab_ptr)
{
/* This code is not intended to be a general hyperslab-to-real
* converter yet. That is why it is not public.
*/
double real_min, real_max; /* Minimum and maximum values */
hsize_t index;
hsize_t total;
double voxel_range, voxel_offset;
double real_range, real_offset;
int i;
double tmp_val;
unsigned long pos[MI2_MAX_VAR_DIMS];
int r;
total = 1;
for (i = 0; i < volume->number_of_dims; i++) {
total *= count[i];
pos[i] = start[i];
}
voxel_offset = volume->valid_min;
voxel_range = volume->valid_max - volume->valid_min;
/* Get the initial real minimum & maximum.
*/
r = miget_slice_range(volume, pos, volume->number_of_dims,
&real_max, &real_min);
if (r == MI_ERROR) {
fprintf(stderr, "Oops - failed to get slice range\n");
}
real_offset = real_min;
real_range = real_max - real_min;
for (index = 0; index < total; index++) {
/* Since this calculation may cross slice boundaries, I need to
* grab the correct real minimum and maximum for the coordinates
* I happen to be in at the time.
*
* This next loop attempts to keep track of the current position,
* and reloads the minimum and maximum whenever we change any other
* than the fastest-varying dimension.
*/
for (i = volume->number_of_dims - 1; i >= 0; i--) {
pos[i]++;
if (pos[i] >= count[i]) {
pos[i] = start[i];
r = miget_slice_range(volume, pos, volume->number_of_dims,
&real_max, &real_min);
if (r == MI_ERROR) {
fprintf(stderr, "Oops - failed to get slice range\n");
}
real_offset = real_min;
real_range = real_max - real_min;
}
else {
break;
}
}
tmp_val = (slab_ptr[index] - voxel_offset) / voxel_range;
slab_ptr[index] = (tmp_val * real_range) + real_offset;
}
}
/** Update an individual thumbnail for the \a volume. Updates group
* number \a ogrp from source group \a igrp. The whole image tree must
* be rooted at \a loc_id.
*/
int
minc_update_thumbnail(mihandle_t volume, hid_t loc_id, int igrp, int ogrp)
{
hsize_t isize[MI2_MAX_VAR_DIMS];
hsize_t osize[MI2_MAX_VAR_DIMS];
hsize_t count[MI2_MAX_VAR_DIMS];
H5_START_T start[MI2_MAX_VAR_DIMS];
hid_t idst_id; /* Input dataset */
hid_t odst_id; /* Output dataset */
hid_t ifspc_id; /* Input "file" dataspace */
hid_t ofspc_id; /* Output "file" dataspace */
hid_t typ_id; /* Type ID */
hid_t imspc_id; /* Input memory dataspace */
hid_t omspc_id; /* Output memory dataspace */
char path[MI2_MAX_PATH];
int ndims; /* Number of dimensions in the image */
int scale;
int i; /* Generic loop counter */
double *in_ptr;
double *out_ptr;
int slice;
int in_bytes;
int out_bytes;
double smax, smin; /* Slice minimum and maximum */
hid_t omax_id; /* Output image-max dataset */
hid_t omin_id; /* Output image-min dataset */
hid_t tfspc_id; /* Dimensionality of image-max/image-min */
hid_t tmspc_id;
hid_t dcpl_id; /* Dataset creation property list */
miinit();
/* Check arguments for basic validity. */
if (ogrp <= igrp) {
return (MI_ERROR);
}
/* Calculate scale factor (always a power of 2) */
for (i = igrp, scale = 1; i < ogrp; i++, scale <<= 1)
;
/* Open the input path.
*/
sprintf(path, "%d/image", igrp);
idst_id = H5Dopen(loc_id, path);
if (idst_id < 0) {
return (MI_ERROR);
}
/* Get the input type.
*/
typ_id = H5Dget_type(idst_id);
/* Get the input dataspace.
*/
ifspc_id = H5Dget_space(idst_id);
/* Get the input dataset creation property list
*/
dcpl_id = H5Dget_create_plist(idst_id);
ndims = H5Sget_simple_extent_ndims(ifspc_id);
H5Sget_simple_extent_dims(ifspc_id, isize, NULL);
/* Calculate the size of the new thumbnail.
*/
for (i = 0; i < ndims; i++) {
osize[i] = isize[i] / scale;
if (osize[i] == 0) { /* Too small? */
return (MI_ERROR);
}
}
/* Create dataspace for new resolution
*/
ofspc_id = H5Screate_simple(ndims, osize, NULL);
sprintf(path, "%d/image", ogrp);
H5E_BEGIN_TRY {
odst_id = H5Dcreate(loc_id, path, typ_id, ofspc_id, H5P_DEFAULT);
} H5E_END_TRY;
if (odst_id < 0) {
odst_id = H5Dopen(loc_id, path);
if (odst_id < 0) {
return (MI_ERROR);
}
}
H5Pclose(dcpl_id); /* No longer needed. */
if (volume->volume_class == MI_CLASS_REAL) {
/* TODO: This is a bit of a hack - I need a better way to get the
* dimensionality of the source image-min and image-max.
*/
tfspc_id = H5Screate_simple(1, &osize[0], NULL);
/* Create a simple scalar dataspace. */
tmspc_id = H5Screate(H5S_SCALAR);
sprintf(path, "%d/image-max", ogrp);
H5E_BEGIN_TRY {
omax_id = H5Dcreate(loc_id, path, H5T_IEEE_F64LE, tfspc_id,
H5P_DEFAULT);
} H5E_END_TRY;
if (omax_id < 0) {
omax_id = H5Dopen(loc_id, path);
}
sprintf(path, "%d/image-min", ogrp);
H5E_BEGIN_TRY {
omin_id = H5Dcreate(loc_id, path, H5T_IEEE_F64LE, tfspc_id,
H5P_DEFAULT);
} H5E_END_TRY;
if (omin_id < 0) {
omin_id = H5Dopen(loc_id, path);
}
}
/* Calculate the input buffer size - scale slices.
*/
in_bytes = scale * isize[1] * isize[2] * sizeof(double);
in_ptr = malloc(in_bytes);
out_bytes = osize[1] * osize[2] * sizeof(double);
out_ptr = malloc(out_bytes);
count[0] = scale;
count[1] = isize[1];
count[2] = isize[2];
imspc_id = H5Screate_simple(ndims, count, NULL);
count[0] = 1;
count[1] = osize[1];
count[2] = osize[2];
omspc_id = H5Screate_simple(ndims, count, NULL);
/*
* read image & TODO: convert to "real" range.
*/
for (slice = 0; slice < osize[0]; slice++) {
start[0] = slice * scale;
start[1] = 0;
start[2] = 0;
count[0] = scale;
count[1] = isize[1];
count[2] = isize[2];
H5Sselect_hyperslab(ifspc_id, H5S_SELECT_SET, start, NULL, count,
NULL);
H5Dread(idst_id, H5T_NATIVE_DOUBLE, imspc_id, ifspc_id, H5P_DEFAULT,
in_ptr);
/* Scale slice from voxel to real values. */
miconvert_hyperslab_to_real(volume, start, count, in_ptr);
midownsample_slice(in_ptr, out_ptr, isize, osize, scale);
start[0] = slice;
start[1] = 0;
start[2] = 0;
count[0] = 1;
count[1] = osize[1];
count[2] = osize[2];
H5Sselect_hyperslab(ofspc_id, H5S_SELECT_SET, start, NULL, count,
NULL);
miconvert_hyperslab_to_voxel(volume, start, count, out_ptr,
&smax, &smin);
H5Dwrite(odst_id, H5T_NATIVE_DOUBLE, omspc_id, ofspc_id, H5P_DEFAULT,
out_ptr);
if (volume->volume_class == MI_CLASS_REAL) {
/* Select the right point in tfspc_id */
H5Sselect_elements(tfspc_id, H5S_SELECT_SET, 1,
(const H5_START_T **) &start[0]);
H5Dwrite(omax_id, H5T_NATIVE_DOUBLE, tmspc_id, tfspc_id,
H5P_DEFAULT, &smax);
H5Dwrite(omin_id, H5T_NATIVE_DOUBLE, tmspc_id, tfspc_id,
H5P_DEFAULT, &smin);
}
}
free(in_ptr);
free(out_ptr);
H5Dclose(omax_id);
H5Dclose(omin_id);
H5Sclose(tfspc_id);
H5Sclose(tmspc_id);
H5Sclose(omspc_id);
H5Sclose(imspc_id);
H5Dclose(odst_id);
H5Tclose(typ_id);
H5Sclose(ofspc_id);
H5Sclose(ifspc_id);
return (MI_NOERROR);
}
/** Cycle through and update each of the lower-resolution images in
* the file.
*/
int
minc_update_thumbnails(mihandle_t volume)
{
int grp_no, prv_grp_no;
hid_t grp_id;
hsize_t n;
hsize_t i;
char name[MI2_MAX_PATH];
grp_id = H5Gopen(volume->hdf_id, "/minc-2.0/image");
if (grp_id < 0) {
return (MI_ERROR); /* Error opening group. */
}
if (H5Gget_num_objs(grp_id, &n) < 0) {
return (MI_ERROR); /* Error getting object count. */
}
grp_no = -1;
for (i = 0; i < n; i++) {
if (H5Gget_objname_by_idx(grp_id, i, name, MI2_MAX_PATH) < 0) {
return (MI_ERROR);
}
prv_grp_no = grp_no;
grp_no = atoi(name);
if (grp_no != 0) {
minc_update_thumbnail(volume, grp_id, prv_grp_no, grp_no);
}
}
H5Gclose(grp_id);
return (MI_NOERROR);
}
double *
alloc1d(int n)
{
return ((double *) malloc(sizeof(double) * n));
}
double **
alloc2d(int n, int m)
{
double **mat;
int i;
mat = (double **) malloc(n * sizeof(double *));
if (mat == NULL) {
return NULL;
}
for (i = 0; i < n; i++) {
mat[i] = (double *) malloc(m * sizeof(double));
if (mat[i] == NULL) {
return NULL;
}
}
return (mat);
}
void
free2d(int n, double **mat)
{
int i;
for (i = 0; i < n; i++) {
free(mat[i]);
}
free(mat);
}
/** Performs scaled maximal pivoting gaussian elimination as a
numerically robust method to solve systems of linear equations.
*/
int
scaled_maximal_pivoting_gaussian_elimination(int n,
int row[],
double **a,
int n_values,
double **solution )
{
int i, j, k, p, v, tmp;
double *s, val, best_val, m, scale_factor;
int success;
s = alloc1d( n );
for ( i = 0; i < n; i++ )
row[i] = i;
for ( i = 0; i < n; i++ ) {
s[i] = fabs( a[i][0] );
for ( j = 1; j < n; j++ ) {
if ( fabs(a[i][j]) > s[i] )
s[i] = fabs(a[i][j]);
}
if ( s[i] == 0.0 ) {
free( s );
return ( FALSE );
}
}
success = TRUE;
for ( i = 0; i < n-1; i++ ) {
p = i;
best_val = a[row[i]][i] / s[row[i]];
best_val = fabs( best_val );
for ( j = i+1; j < n; j++ ) {
val = a[row[j]][i] / s[row[j]];
val = fabs( val );
if ( val > best_val ) {
best_val = val;
p = j;
}
}
if ( a[row[p]][i] == 0.0 ) {
success = FALSE;
break;
}
if ( i != p ) {
tmp = row[i];
row[i] = row[p];
row[p] = tmp;
}
for ( j = i+1; j < n; j++ ) {
if ( a[row[i]][i] == 0.0 ) {
success = FALSE;
break;
}
m = a[row[j]][i] / a[row[i]][i];
for ( k = i+1; k < n; k++ )
a[row[j]][k] -= m * a[row[i]][k];
for ( v = 0; v < n_values; v++ )
solution[row[j]][v] -= m * solution[row[i]][v];
}
if ( !success )
break;
}
if ( success && a[row[n-1]][n-1] == 0.0 )
success = FALSE;
if ( success ) {
for ( i = n-1; i >= 0; --i ) {
for ( j = i+1; j < n; j++ ) {
scale_factor = a[row[i]][j];
for ( v = 0; v < n_values; v++ )
solution[row[i]][v] -= scale_factor * solution[row[j]][v];
}
for ( v = 0; v < n_values; v++ )
solution[row[i]][v] /= a[row[i]][i];
}
}
free( s );
return( success );
}
int
scaled_maximal_pivoting_gaussian_elimination_real(int n,
double **coefs,
int n_values,
double **values )
{
int i, j, v, *row;
double **a, **solution;
int success;
row = (int *) alloc1d( n ); /* Ugly and wasteful but OK for 1D array */
a = alloc2d( n, n );
solution = alloc2d( n, n_values );
for (i = 0; i < n; i++) {
for ( j = 0; j < n; j++ )
a[i][j] = coefs[i][j];
for ( v = 0; v < n_values; v++ )
solution[i][v] = values[v][i];
}
success = scaled_maximal_pivoting_gaussian_elimination( n, row, a,
n_values,
solution );
if ( success ) {
for ( i = 0; i < n; i++ ) {
for ( v = 0; v < n_values; v++ ) {
values[v][i] = solution[row[i]][v];
}
}
}
free(a);
free(solution);
free(row);
return ( success );
}
/** Computes the inverse of a square matrix.
*/
static int
invert_4x4_matrix(double matrix[4][4], /**< Input matrix */
double inverse[4][4]) /**< Output (inverted) matrix */
{
double tmp;
int result;
int i, j;
double **mtmp;
double **itmp;
mtmp = alloc2d(4, 4);
itmp = alloc2d(4, 4);
/* Start off with the identity matrix. */
for ( i = 0; i < 4; i++ ) {
for ( j = 0; j < 4; j++ ) {
itmp[i][j] = 0.0;
mtmp[i][j] = matrix[i][j];
}
itmp[i][i] = 1.0;
}
result = scaled_maximal_pivoting_gaussian_elimination_real( 4, mtmp,
4, itmp );
if ( result ) {
for ( i = 0; i < 4; i++) {
for ( j = 0; j < 4; j++ ) {
inverse[i][j] = itmp[j][i];
}
}
}
free(mtmp);
free(itmp);
return (result ? MI_NOERROR : MI_ERROR);
}
/** Fills in the transform with the identity matrix.
*/
static void
mimake_identity_transform(mi_lin_xfm_t transform)
{
int i, j;
for (i = 0; i < MI2_LIN_XFM_SIZE; i++) {
for (j = 0; j < MI2_LIN_XFM_SIZE; j++) {
transform[i][j] = 0.0;
}
transform[i][i] = 1.0;
}
}
/** Function for inverting a MINC linear transform.
*/
int
miinvert_transform(mi_lin_xfm_t transform, mi_lin_xfm_t inverse )
{
int result;
result = invert_4x4_matrix( transform, inverse );
if (result != MI_NOERROR) {
mimake_identity_transform(inverse);
}
return ( result );
}
/** Function to read a scalar variable of HDF5 type \a type_id from the given
* \a path relative to the HDF5 file or group \a loc_id.
*/
int
miget_scalar(hid_t loc_id, hid_t type_id, const char *path, void *data)
{
hid_t dset_id;
hid_t spc_id;
int result = MI_ERROR;
H5E_BEGIN_TRY {
dset_id = H5Dopen(loc_id, path);
} H5E_END_TRY;
if (dset_id >= 0) {
spc_id = H5Dget_space(dset_id);
if (spc_id >= 0) {
if (H5Sget_simple_extent_ndims(spc_id) == 0) {
if (H5Dread(dset_id, type_id, H5S_ALL, H5S_ALL, H5P_DEFAULT,
data) >= 0) {
result = MI_NOERROR;
}
}
H5Sclose(spc_id);
}
H5Dclose(dset_id);
}
return (result);
}
syntax highlighted by Code2HTML, v. 0.9.1