Mercurial > hg > minc-tools
changeset 1863:f8c0ff2d4cdd
Implement miset_world_origin, change names of miconvert_3D_world_to_voxel() and miconvert_3D_voxel_to_world() to drop the _3D_ part
author | bert <bert> |
---|---|
date | Wed, 11 Aug 2004 20:40:13 +0000 |
parents | 0ca40a05853a |
children | d4b7f4d2e643 |
files | libsrc2/convert.c |
diffstat | 1 files changed, 271 insertions(+), 78 deletions(-) [+] |
line wrap: on
line diff
--- a/libsrc2/convert.c +++ b/libsrc2/convert.c @@ -6,6 +6,10 @@ * data, and to convert coordinates between "voxel" and "world" systems. */ +/** + * \defgroup mi2Cvt MINC 2.0 Coordinate and Voxel Conversion Functions + */ + #include <stdlib.h> #include <hdf5.h> #include <math.h> @@ -14,10 +18,22 @@ #include "minc2.h" #include "minc2_private.h" -/** convert values between real (scaled) values and voxel (unscaled) +/** Convert values between real (scaled) values and voxel (unscaled) * values. The voxel value is the unscaled value, and corresponds to the * value actually stored in the file, whereas the "real" value is the * value at the given location after scaling has been applied. + * + * The \a coords parameter specifies the location at which the + * conversion is performed. This is needed because MINC supports + * per-slice scaling, therefore a conversion performed at one location + * may differ from that performed at another location. + * + * \param volume A volume handle + * \param coords The position for which to perform the conversion. + * \param ncoords The length of the \a coords array. + * \param real_value The original real value, to be converted to voxel. + * \param voxel_value_ptr A pointer to the converted voxel value. + * \ingroup mi2Cvt */ int miconvert_real_to_voxel(mihandle_t volume, @@ -55,10 +71,22 @@ } -/** convert values between real (scaled) values and voxel (unscaled) +/** Convert values between real (scaled) values and voxel (unscaled) * values. The voxel value is the unscaled value, and corresponds to the * value actually stored in the file, whereas the "real" value is the * value at the given location after scaling has been applied. + * + * The \a coords parameter specifies the location at which the + * conversion is performed. This is needed because MINC supports + * per-slice scaling, therefore a conversion performed at one location + * may differ from that performed at another location. + * + * \param volume A volume handle + * \param coords The position for which to perform the conversion. + * \param ncoords The length of the \a coords array. + * \param voxel_value The original voxel value, to be converted to real. + * \param real_value_ptr A pointer to the converted real value. + * \ingroup mi2Cvt */ int miconvert_voxel_to_real(mihandle_t volume, @@ -93,126 +121,103 @@ return (result); } -/** +/** \internal */ static void -mireorder_voxel_to_xyz(mihandle_t volume, const double voxel[MI2_3D], - double xyz[MI2_3D]) +mireorder_voxel_to_xyz(mihandle_t volume, + const double voxel[], + double xyz[MI2_3D], + midimclass_t dimclass) { int i, axis; for (i = 0; i < volume->number_of_dims; i++) { - axis = volume->world_indices[i]; - if ( axis >= 0 ) - xyz[i] = voxel[axis]; - else - xyz[i] = 0.0; + midimhandle_t hdim = volume->dim_handles[i]; + axis = hdim->world_index; + if ( axis >= 0 && dimclass == hdim->class) { + xyz[axis] = voxel[i]; + } } } -/** +/** \internal */ static void -mireorder_xyz_to_voxel(mihandle_t volume, const double xyz[MI2_3D], - double voxel[MI2_3D]) +mireorder_xyz_to_voxel(mihandle_t volume, + const double xyz[MI2_3D], + double voxel[], + midimclass_t dimclass) { int i, axis; for (i = 0; i < volume->number_of_dims; i++) { - axis = volume->world_indices[i]; - if ( axis >= 0 ) { - voxel[axis] = xyz[i]; + midimhandle_t hdim = volume->dim_handles[i]; + axis = hdim->world_index; + if ( axis >= 0 && dimclass == hdim->class) { + voxel[i] = xyz[axis]; } } } -/** Converts a 3-dimensional spatial position in voxel coordinates into a +/** Converts an N-dimensional spatial position in voxel coordinates into a * 3-dimensional spatial position in world coordinates. * * The returned world coordinate vector is in a standardized order, with * the X position first (at index 0), followed by the Y and Z coordinates. * The voxel coordinate vector is in the native order appropriate to the * file. + * + * \ingroup mi2Cvt */ int -miconvert_3D_voxel_to_world(mihandle_t volume, - const double voxel[MI2_3D], - double world[MI2_3D]) +miconvert_voxel_to_world(mihandle_t volume, + const double voxel[], + double world[MI2_3D]) { double temp[MI2_3D]; - mireorder_voxel_to_xyz(volume, voxel, temp); + mireorder_voxel_to_xyz(volume, voxel, temp, MI_DIMCLASS_SPATIAL); mitransform_coord(world, volume->v2w_transform, temp); return (MI_NOERROR); } /** Converts a 3-dimensional spatial position in world coordinates into a - * 3-dimensional spatial position in voxel coordinates. + * N-dimensional spatial position in voxel coordinates. * * The input world coordinate vector is in a standardized order, with * the X position first (at index 0), followed by the Y and Z coordinates. * The voxel coordinate vector is in the native order appropriate to the * file. + * + * \ingroup mi2Cvt */ -int miconvert_3D_world_to_voxel(mihandle_t volume, - const double world[MI2_3D], - double voxel[MI2_3D]) +int +miconvert_world_to_voxel(mihandle_t volume, + const double world[MI2_3D], + double voxel[]) { double temp[MI2_3D]; + int i; + + for (i = 0; i < volume->number_of_dims; i++) { + voxel[i] = 0.0; + } mitransform_coord(temp, volume->w2v_transform, world); - mireorder_xyz_to_voxel(volume, temp, voxel); - return (MI_NOERROR); -} - -/** Convert a */ -int -miconvert_3D_voxel_to_spatial_frequency(mihandle_t volume, - const double voxel[MI2_3D], - double world[MI2_3D]) - -{ - return (MI_NOERROR); -} - -int -miconvert_3D_spatial_frequency_to_voxel(mihandle_t volume, - const double world[MI2_3D], - double voxel[MI2_3D]) - -{ - return (MI_NOERROR); -} - -/** - * This function converts coordinates in voxel coordinates to world - * coordinates. Only the spatial coordinates are converted. - */ -int -miconvert_voxel_to_world(midimhandle_t dimensions[], - int ndims, - const double voxel[], - double world[]) -{ - return (MI_NOERROR); -} - -/** - * This function converts coordinates in world coordinates to voxel - * coordinates. Only the spatial coordinates are converted. - */ -int -miconvert_world_to_voxel(midimhandle_t dimensions[], - int ndims, - const double voxel[], - double world[]) -{ + mireorder_xyz_to_voxel(volume, temp, voxel, MI_DIMCLASS_SPATIAL); return (MI_NOERROR); } /** This function retrieves the real values of a position in the * MINC volume. The "real" value is the value at the given location * after scaling has been applied. + * + * \param volume A volume handle + * \param coords The voxel position to retrieve + * \param ndims The number of values in the \a coords array + * \param value_ptr Pointer to a double variable to hold the returned value. + * + * \ingroup mi2Cvt */ int miget_real_value(mihandle_t volume, @@ -234,6 +239,13 @@ /** This function sets the real value of a position in the MINC * volume. The "real" value is the value at the given location * after scaling has been applied. + * + * \param volume A volume handle + * \param coords The voxel position to retrieve + * \param ndims The number of values in the \a coords array + * \param value The value to save at this location. + * + * \ingroup mi2Cvt */ int miset_real_value(mihandle_t volume, @@ -243,19 +255,37 @@ { double voxel; + if ((volume->mode & MI2_OPEN_RDWR) == 0) { + return (MI_ERROR); + } + miconvert_real_to_voxel(volume, coords, ndims, value, &voxel); return (miset_voxel_value(volume, coords, ndims, voxel)); } +/** + * This function calculates the start values for the volume dimensions, + * assuming that the spatial origin is relocated to the given world + * coordinate. + * + * \ingroup mi2Cvt + */ int -miconvert_spatial_origin_to_start( mihandle_t volume, - double world[3], - double starts[3]) +miconvert_world_origin_to_start( mihandle_t volume, + double world[3], + double starts[3]) { return (MI_NOERROR); } +/** + * This function calculates the start values for the volume dimensions, + * assuming that the spatial origin is relocated to the given world + * coordinate. + * + * \ingroup mi2Cvt + */ int miconvert_spatial_frequency_origin_to_start( mihandle_t volume, double world[3], @@ -264,23 +294,178 @@ return (MI_NOERROR); } +static double +dot_vectors(int n, double v1[], double v2[]) +{ + int i; + double d; + + d = 0.0; + for ( i = 0; i < n; i++ ) { + d += v1[i] * v2[i]; + } + return ( d ); +} + int -miset_spatial_origin_to_start( mihandle_t volume, - double world[3]) +solve_linear_system( int n, double **coefs, double *values, double *solution) +{ + int i; + + for ( i = 0; i < n; i++ ) { + solution[i] = values[i]; + } + + return (scaled_maximal_pivoting_gaussian_elimination_real(n, coefs, 1, + &solution )); +} + +static void +convert_transform_origin_to_starts(mihandle_t hvol, + double origin[], + double starts[] ) { + int axis; + int which[MI2_3D]; + int n_axes, i, j; + double o_dot_c, c_dot_c; + double x_dot_x, x_dot_y, x_dot_v, y_dot_y, y_dot_v, bottom; + double **matrix, solution[MI2_3D]; + + for ( i = 0; i < hvol->number_of_dims; i++) { + starts[i] = 0.0; + } + + /*--- get the list of valid axes (which) */ + + n_axes = 0; + for ( i = 0; i < hvol->number_of_dims; i++ ) { + axis = hvol->dim_handles[i]->world_index; + if ( axis >= 0 ) { + which[axis] = i; + ++n_axes; + } + } + + /*--- get the starts: computed differently for 1, 2, or 3 axes */ + + switch (n_axes) { + case 1: + o_dot_c = dot_vectors(MI2_3D, + origin, + hvol->dim_handles[which[0]]->direction_cosines); + c_dot_c = dot_vectors(MI2_3D, + hvol->dim_handles[which[0]]->direction_cosines, + hvol->dim_handles[which[0]]->direction_cosines); + if ( c_dot_c != 0.0 ) { + starts[which[0]] = o_dot_c / c_dot_c; + } + break; + + case 2: + x_dot_x = dot_vectors(MI2_3D, + hvol->dim_handles[which[0]]->direction_cosines, + hvol->dim_handles[which[0]]->direction_cosines ); + x_dot_v = dot_vectors(MI2_3D, + hvol->dim_handles[which[0]]->direction_cosines, + origin ); + x_dot_y = dot_vectors(MI2_3D, + hvol->dim_handles[which[0]]->direction_cosines, + hvol->dim_handles[which[1]]->direction_cosines ); + y_dot_y = dot_vectors(MI2_3D, + hvol->dim_handles[which[1]]->direction_cosines, + hvol->dim_handles[which[1]]->direction_cosines ); + y_dot_v = dot_vectors(MI2_3D, + hvol->dim_handles[which[1]]->direction_cosines, + origin ); + + bottom = x_dot_x * y_dot_y - x_dot_y * x_dot_y; + + if ( bottom != 0.0 ) { + starts[which[0]] = (x_dot_v * y_dot_y - x_dot_y * y_dot_v) / bottom; + starts[which[1]] = (y_dot_v * x_dot_x - x_dot_y * x_dot_v) / bottom; + } + break; + + case 3: + /*--- this is the usual case, solve the equations to find what + starts give the desired origin */ + + matrix = alloc2d(MI2_3D, MI2_3D); + + for ( i = 0; i < MI2_3D; i++) { + for ( j = 0; j < hvol->number_of_dims; j++ ) { + matrix[i][j] = hvol->dim_handles[j]->direction_cosines[i]; + } + } + + if ( solve_linear_system( MI2_3D, matrix, origin, solution ) ) { + starts[which[0]] = solution[0]; + starts[which[1]] = solution[1]; + starts[which[2]] = solution[2]; + } + + free(matrix); + + break; + } +} +/** + * This function sets the world coordinates of the point (0,0,0) in voxel + * coordinates. This changes the constant offset of the two coordinate + * systems. + * + * \ingroup mi2Cvt + */ +int +miset_world_origin(mihandle_t volume, /**< A volume handle */ + double world[MI2_3D]) /**< The world coordinates of voxel origin */ +{ + double starts[MI2_3D]; + int i; + + convert_transform_origin_to_starts(volume, world, starts); + for (i = 0; i < volume->number_of_dims; i++) { + midimhandle_t hdim = volume->dim_handles[i]; + if (hdim->class == MI_DIMCLASS_SPATIAL || + hdim->class == MI_DIMCLASS_SFREQUENCY) { + hdim->start = starts[hdim->world_index]; + } + } + + /* Get the voxel to world transform for the volume + */ + miget_voxel_to_world(volume, volume->v2w_transform); + + /* Calculate the inverse transform */ + miinvert_transform(volume->v2w_transform, volume->w2v_transform); + return (MI_NOERROR); } +/** + * This function sets the world coordinates of the point (0,0,0) in voxel + * coordinates. This changes the constant offset of the two coordinate + * systems. + * + * \ingroup mi2Cvt + */ int -miset_spatial_frequency_origin_to_start(mihandle_t volume, - double world[3]) +miset_spatial_frequency_origin(mihandle_t volume, + double world[3]) { + if ((volume->mode & MI2_OPEN_RDWR) == 0) { + return (MI_ERROR); + } + return (MI_NOERROR); } /** This function retrieves the voxel values of a position in the * MINC volume. The voxel value is the unscaled value, and corresponds * to the value actually stored in the file. + * + * \ingroup mi2Cvt */ int miget_voxel_value(mihandle_t volume, @@ -303,6 +488,8 @@ /** This function sets the voxel value of a position in the MINC * volume. The voxel value is the unscaled value, and corresponds to the * value actually stored in the file. + * + * \ingroup mi2Cvt */ int miset_voxel_value(mihandle_t volume, @@ -314,6 +501,10 @@ unsigned long count[MI2_MAX_VAR_DIMS]; int i; + if ((volume->mode & MI2_OPEN_RDWR) == 0) { + return (MI_ERROR); + } + for (i = 0; i < ndims; i++) { count[i] = 1; } @@ -325,6 +516,8 @@ /** Get the absolute minimum and maximum values of a volume. + * + * \ingroup mi2Cvt */ int miget_volume_real_range(mihandle_t volume, double real_range[])