--- 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[])