/** \file m2util.c * \brief MINC 2.0 Private utility functions * \author Leila Baghdadi, Bert Vincent * ************************************************************************/ #include #include #include #include #include #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); }