--- a/volume_io/Documentation/volume_io.tex
+++ b/volume_io/Documentation/volume_io.tex
@@ -2280,6 +2280,13 @@
 \name{Transform}.
 
 {\bf\begin{verbatim}
+#define   Transform_elem( transform, i, j )
+\end{verbatim}}
+
+\desc{Is used to set or retrieve the value of the \name{i}'th row and
+\name{j}'th column of the transform, where $0 <= i,j <= 3$.}
+
+{\bf\begin{verbatim}
 public  void  make_identity_transform(
     Transform   *transform )
 \end{verbatim}}
@@ -2372,7 +2379,9 @@
 \end{verbatim}}
 
 \desc{Creates a general transform consisting of a single linear
-transform, specified by \name{linear\_transform}.}
+transform, specified by \name{linear\_transform}.  If a \name{NULL} is
+passed as the argument \name{linear\-\_transform}, then an identity
+transform is created.}
 
 {\bf\begin{verbatim}
 public  void  create_thin_plate_transform(
@@ -2409,6 +2418,37 @@
 the general transform.}
 
 {\bf\begin{verbatim}
+typedef  void   (*User_transform_function)( void  *user_data,
+                                            Real  x,
+                                            Real  y,
+                                            Real  z,
+                                            Real  *x_trans,
+                                            Real  *y_trans,
+                                            Real  *z_trans );
+
+public  void  create_user_transform(
+    General_transform         *transform,
+    void                      *user_data,
+    size_t                    size_user_data,
+    User_transform_function   transform_function,
+    User_transform_function   inverse_transform_function )
+\end{verbatim}}
+
+\desc{Creates a user defined transformation, by copying the user data
+(\name{size\-\_user\-\_data} bytes of data starting at
+\name{user\_data}).   Two function pointers are also required, to
+specify the method of transforming points and inversely transforming
+points.  These functions are of the type
+\name{User\-\_transform\-\_function}, which is specified above.}
+
+{\bf\begin{verbatim}
+public  void  delete_general_transform(
+    General_transform   *transform )
+\end{verbatim}}
+
+\desc{Frees up the memory stored in the general transform structure.}
+
+{\bf\begin{verbatim}
 public  Transform_types  get_transform_type(
     General_transform   *transform )
 \end{verbatim}}
@@ -2427,6 +2467,87 @@
 with the routines specific to linear transforms, described earlier.
 Otherwise prints an error message.}
 
+{\bf\begin{verbatim}
+public  void  concat_general_transforms(
+    General_transform   *first,
+    General_transform   *second,
+    General_transform   *result )
+\end{verbatim}}
+
+\desc{Concatenates two general transforms.  If both transforms are of
+type \name{LINEAR}, then the result is also of this type, being the
+matrix product of the two.  Otherwise, the resulting transform is
+simply the concatenation of the two transforms.}
+
+{\bf\begin{verbatim}
+public  int  get_n_concated_transforms(
+    General_transform   *transform )
+\end{verbatim}}
+
+\desc{Returns the number of concatenated transforms in the given
+transform.  If the type of the transform is
+\name{CONCATENATED\_TRANSFORM}, then the number returned is the number
+of transforms, otherwise it is one.}
+
+{\bf\begin{verbatim}
+public  General_transform  *get_nth_general_transform(
+    General_transform   *transform,
+    int                 n )
+\end{verbatim}}
+
+\desc{If the transform is of type \name{CONCATENATED\_TRANSFORM}, then
+a pointer to the \name{n}'th transform is returned, where \name{n}
+is greater than or equal to zero and less than the number of
+transforms in the concatenated transform.}
+
+\section{Using General Transforms}
+
+{\bf\begin{verbatim}
+public  void  general_transform_point(
+    General_transform   *transform,
+    Real                x,
+    Real                y,
+    Real                z,
+    Real                *x_transformed,
+    Real                *y_transformed,
+    Real                *z_transformed )
+\end{verbatim}}
+
+\desc{Transforms a three dimensional point by a general transform,
+passing back the result in the last three arguments.}
+
+{\bf\begin{verbatim}
+public  void  general_inverse_transform_point(
+    General_transform   *transform,
+    Real                x,
+    Real                y,
+    Real                z,
+    Real                *x_transformed,
+    Real                *y_transformed,
+    Real                *z_transformed )
+\end{verbatim}}
+
+\desc{Transforms a three dimensional point by the inverse of the
+general transform, passing back the result in the last three arguments.}
+
+{\bf\begin{verbatim}
+public  void  copy_general_transform(
+    General_transform   *transform,
+    General_transform   *copy )
+\end{verbatim}}
+
+\desc{Creates a copy of the general transform, allocating memory
+within the structure as required.}
+
+{\bf\begin{verbatim}
+public  void  create_inverse_general_transform(
+    General_transform   *transform,
+    General_transform   *inverse )
+\end{verbatim}}
+
+\desc{Creates a new general transform that is the inverse of the given
+one.}
+
 \section{Reading and Writing General Transforms}
 
 General transforms are stored in files in an ascii format devised at

--- a/volume_io/Geometry/gaussian.c
+++ b/volume_io/Geometry/gaussian.c
@@ -10,20 +10,22 @@
     int   n_values,
     Real  **values )
 {
-    int       i, j, k, p, v;
-    Real      **a, *s, val, best_val, m, *tmp, scale_factor;
+    int       i, j, k, p, v, *row, tmp;
+    Real      **a, **solution, *s, val, best_val, m, scale_factor;
     BOOLEAN   success;
 
-    ALLOC( a, n );
-    for_less( i, 0, n )
-        ALLOC( a[i], n );
-
+    ALLOC( row, n );
+    ALLOC2D( a, n, n );
+    ALLOC2D( solution, n, n_values );
     ALLOC( s, n );
 
     for_less( i, 0, n )
     {
+        row[i] = i;
         for_less( j, 0, n )
             a[i][j] = coefs[i][j];
+        for_less( j, 0, n_values )
+            solution[i][j] = values[j][i];
     }
 
     for_less( i, 0, n )
@@ -41,11 +43,11 @@
     for_less( i, 0, n-1 )
     {
         p = i;
-        best_val = a[i][i] / s[i];
+        best_val = a[row[i]][i] / s[row[i]];
         best_val = ABS( best_val );
-        for_less( j, i, n )
+        for_less( j, i+1, n )
         {
-            val = a[j][i] / s[j];
+            val = a[row[j]][i] / s[row[j]];
             val = ABS( val );
             if( val > best_val )
             {
@@ -54,7 +56,7 @@
             }
         }
 
-        if( a[p][i] == 0.0 )
+        if( a[row[p]][i] == 0.0 )
         {
             success = FALSE;
             break;
@@ -62,42 +64,52 @@
 
         if( i != p )
         {
-            tmp = a[i];
-            a[i] = a[p];
-            a[p] = tmp;
+            tmp = row[i];
+            row[i] = row[p];
+            row[p] = tmp;
         }
 
         for_less( j, i+1, n )
         {
-            m = a[j][i] / a[i][i];
+            m = a[row[j]][i] / a[row[i]][i];
             for_less( k, i+1, n )
-                a[j][k] -= m * a[i][k];
+                a[row[j]][k] -= m * a[row[i]][k];
             for_less( v, 0, n_values )
-                values[v][j] -= m * values[v][i];
+                solution[row[j]][v] -= m * solution[row[i]][v];
         }
     }
 
-    if( success && a[n-1][n-1] == 0.0 )
+    if( success && a[row[n-1]][n-1] == 0.0 )
         success = FALSE;
 
-    for( i = n-1;  i >= 0;  --i )
+    if( success )
     {
-        for_less( j, i+1, n )
+        for( i = n-1;  i >= 0;  --i )
         {
-            scale_factor = a[i][j];
+            for_less( j, i+1, n )
+            {
+                scale_factor = a[row[i]][j];
+                for_less( v, 0, n_values )
+                    solution[row[i]][v] -= scale_factor * solution[row[j]][v];
+            }
+
             for_less( v, 0, n_values )
-                values[v][i] -= scale_factor * values[v][j];
+                solution[row[i]][v] /= a[row[i]][i];
         }
 
-        for_less( v, 0, n_values )
-            values[v][i] /= a[i][i];
+        for_less( i, 0, n )
+        {
+            for_less( v, 0, n_values )
+            {
+                values[v][i] = solution[row[i]][v];
+            }
+        }
     }
 
-    for_less( i, 0, n )
-        FREE( a[i] );
-    FREE( a );
-
+    FREE2D( a );
+    FREE2D( solution );
     FREE( s );
+    FREE( row );
 
     return( success );
 }
@@ -122,6 +134,8 @@
     Real  **matrix,
     Real  **inverse )
 {
+    Real      tmp;
+    BOOLEAN   success;
     int       i, j;
 
     for_less( i, 0, n )
@@ -131,6 +145,21 @@
         inverse[i][i] = 1.0;
     }
 
-    return( scaled_maximal_pivoting_gaussian_elimination( n, matrix,
-                                                          n, inverse ) );
+    success = scaled_maximal_pivoting_gaussian_elimination( n, matrix,
+                                                            n, inverse );
+
+    if( success )
+    {
+        for_less( i, 0, n-1 )
+        {
+            for_less( j, i+1, n )
+            {
+                tmp = inverse[i][j];
+                inverse[i][j] = inverse[j][i];
+                inverse[j][i] = tmp;
+            }
+        }
+    }
+
+    return( success );
 }

--- a/volume_io/Geometry/inverse.c
+++ b/volume_io/Geometry/inverse.c
@@ -1,6 +1,6 @@
 
 #ifndef lint
-static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Geometry/inverse.c,v 1.4 1994-11-25 14:19:29 david Exp $";
+static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Geometry/inverse.c,v 1.5 1995-02-20 13:12:07 david Exp $";
 #endif
 
 #include  <internal_volume_io.h>
@@ -22,37 +22,28 @@
     Transform  *transform,
     Transform  *inverse )
 {
-    int        i, j, ind[5];
-    float      **t, **inv, col[5], d;
+    int        i, j;
+    Real       **t, **inv;
     Transform  ident;
+    BOOLEAN    success;
 
     /* --- copy the transform to a numerical recipes type matrix */
 
-    ALLOC2D( t, 5, 5 );
-    ALLOC2D( inv, 5, 5 );
+    ALLOC2D( t, 4, 4 );
+    ALLOC2D( inv, 4, 4 );
 
     for_less( i, 0, 4 )
     {
         for_less( j, 0, 4 )
         {
-            t[i+1][j+1] = Transform_elem(*transform,i,j);
+            t[i][j] = Transform_elem(*transform,i,j);
         }
     }
 
-    ludcmp( t, 4, ind, &d );
+    success = invert_square_matrix( 4, t, inv );
 
-    if( d != 0.0 )
+    if( success )
     {
-        for_inclusive( j, 1, 4 )
-        {
-            for_inclusive( i, 1, 4 )
-                col[i] = 0.0;
-            col[j] = 1.0;
-            lubksb( t, 4, ind, col );
-            for_inclusive( i, 1, 4 )
-                inv[i][j] = col[i];
-        }
-
         /* --- copy the resulting numerical recipes matrix to the
                output argument */
 
@@ -60,7 +51,7 @@
         {
             for_less( j, 0, 4 )
             {
-                Transform_elem(*inverse,i,j) = inv[i+1][j+1];
+                Transform_elem(*inverse,i,j) = inv[i][j];
             }
         }
 
@@ -73,9 +64,11 @@
             print( "Error in compute_transform_inverse\n" );
         }
     }
+    else
+        make_identity_transform( inverse );
 
     FREE2D( t );
     FREE2D( inv );
 
-    return( d != 0.0 );
+    return( success );
 }

--- a/volume_io/MNI_formats/thin_plate_spline.c
+++ b/volume_io/MNI_formats/thin_plate_spline.c
@@ -1,6 +1,6 @@
 
 #ifndef lint
-static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/MNI_formats/thin_plate_spline.c,v 1.4 1994-11-25 14:20:26 david Exp $";
+static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/MNI_formats/thin_plate_spline.c,v 1.5 1995-02-20 13:12:27 david Exp $";
 #endif
 
 /* ----------------------------- MNI Header -----------------------------------
@@ -26,7 +26,7 @@
 private  void mnewt(int ntrial,float x[], int dim, float tolx, float tolf,
                     float **bdefor, float **INVMLY, int num_marks);
 private  void 
-  usrfun(float x[], float **alpha, float bet[], float **bdefor, float **INVMLY, 
+  usrfun(float x[], Real **alpha, Real bet[], float **bdefor, float **INVMLY, 
 	 int num_marks, int dim, float *xout);
 private  float return_r(float *cor1, float *cor2, int dim);
 private  float FU(float r, int dim);
@@ -48,7 +48,7 @@
               INVMLY to 'in' to get 'out'
 @METHOD     : 
 @GLOBALS    : none
-@CALLS      : return_r, vector, free_vector
+@CALLS      : return_r
 @CREATED    : Mon Apr  5 09:00:54 EST 1993
 @MODIFIED   : 
 ---------------------------------------------------------------------------- */
@@ -188,31 +188,28 @@
 	      either summed variable increments tolx or summed function values tolf
               
 @GLOBALS    : none
-@CALLS      : ludcmp, lubksb, usrfun
+@CALLS      : usrfun
 @CREATED    : Mon Apr  5 09:00:54 EST 1993
 @MODIFIED   : 
 ---------------------------------------------------------------------------- */
 
-#define FREERETURN {FREE2D(alpha);FREE(bet);FREE(indx);return;}
+#define FREERETURN {FREE2D(alpha);FREE(bet);return;}
 
 private  void mnewt(int ntrial,float x[], int dim, float tolx, float tolf,
                     float **bdefor, float **INVMLY, int num_marks)
 
 {
-  int k,i,*indx;
-  float 
-    errx,errf,d,
-    *bet,**alpha,
-    xout[N_DIMENSIONS];
+  int   k,i;
+  Real  *bet, **alpha;
+  float errx,errf, xout[N_DIMENSIONS];
   
    xout[0] = x[0];          /* assume that the result is equal to the input, */
    xout[1] = x[1];	    /* as if the transformation was identity         */
    if (dim>2)
      xout[2] = x[2];
 
-   ALLOC( indx, dim + 1 );
-   ALLOC( bet, dim + 1 );
-   ALLOC2D( alpha, dim + 1, dim + 1 );
+   ALLOC( bet, dim );
+   ALLOC2D( alpha, dim, dim );
 
    for (k=1;k<=ntrial;k++) {   
 			/* usrfun will build the matrix coefficients for
@@ -220,14 +217,16 @@
       usrfun(x, alpha, bet, bdefor, INVMLY, num_marks, dim, xout);
 
       errf=0.0;			/* Check for function convergence */
-      for (i=1;i<=dim;i++) errf += fabs(bet[i]);
+      for (i=1;i<=dim;i++) errf += fabs(bet[i-1]);
       if (errf <= tolf) FREERETURN
-      ludcmp(alpha,dim,indx,&d); /* Solve the linear eqs using LU decomposition */
-      lubksb(alpha,dim,indx,bet);
+
+      if( !solve_linear_system( dim, alpha, bet, bet ) )
+          FREERETURN;
+
       errx=0.0;
       for (i=1;i<=dim;i++) {	/* check for root convergence */
-	 errx += fabs(bet[i]);
-	 x[i-1] += bet[i];	/* Update the solution */
+	 errx += fabs(bet[i-1]);
+	 x[i-1] += bet[i-1];	/* Update the solution */
       }
       if (errx <= tolx) FREERETURN
    }
@@ -261,7 +260,7 @@
 ---------------------------------------------------------------------------- */
 
 private  void 
-  usrfun(float x[], float **alpha, float bet[], float **bdefor, float **INVMLY, 
+  usrfun(float x[], Real **alpha, Real bet[], float **bdefor, float **INVMLY, 
 	 int num_marks, int dim, float *xout)
 
 /* NOTE: Now this function will only work for 3-D case */
@@ -271,8 +270,8 @@
 
    ALLOC( d, num_marks );
 
-   for (i=1; i<=dim; i++){	/* build up the matrix of linear parameters */
-      for (j=1; j<=dim; j++){
+   for (i=0; i<dim; i++){	/* build up the matrix of linear parameters */
+      for (j=0; j<dim; j++){
          alpha[i][j] = 0;
       }
    }
@@ -281,25 +280,25 @@
    for (i=0; i<num_marks; i++)  /* build distances array */
        d[i] = return_r( x, bdefor[i], dim );
 
-   for (i = 1; i<= dim; i++){	/* Estimate the linear parameters for the error function */
-      for (j = 1; j <= dim; j++){
+   for (i = 0; i< dim; i++){	/* Estimate the linear parameters for the error function */
+      for (j = 0; j < dim; j++){
          alpha[i][j] = 0;
-         for (k = 1; k <= num_marks; k++){
-            alpha[i][j] += (x[j-1] - bdefor[k-1][j-1])/d[k-1] *INVMLY[k-1][i-1];
+         for (k = 0; k < num_marks; k++){
+            alpha[i][j] += (x[j] - bdefor[k][j])/d[k] *INVMLY[k][i];
          }
-         alpha[i][j] += INVMLY[k+j-1][i-1];
+         alpha[i][j] += INVMLY[k+j][i];
       }
    }
 
-   for (i=1; i<=dim; i++){
+   for (i=0; i<dim; i++){
       bet[i] = 0;
-      for (k=1; k<=num_marks; k++){
-         bet[i] += d[k-1]*INVMLY[k-1][i-1];
+      for (k=0; k<num_marks; k++){
+         bet[i] += d[k]*INVMLY[k][i];
       }
-      bet[i] += INVMLY[k-1][i-1] +
-                INVMLY[k+1-1][i-1]*x[0] +
-                INVMLY[k+2-1][i-1]*x[1] +
-                INVMLY[k+3-1][i-1]*x[2] - xout[i-1];
+      bet[i] += INVMLY[k][i] +
+                INVMLY[k+1][i]*x[0] +
+                INVMLY[k+2][i]*x[1] +
+                INVMLY[k+3][i]*x[2] - xout[i];
       bet[i] = -bet[i];
    }
 

--- a/volume_io/Testing/reorder_volume.c
+++ b/volume_io/Testing/reorder_volume.c
@@ -5,17 +5,16 @@
     char  *argv[] )
 {
     int                 i, a, n_dims, v, n_volumes, sizes[MAX_DIMENSIONS];
-    int                 sizes_2d[MAX_DIMENSIONS];
     char                *input_filename, *output_filename, *string;
     char                *in_dim_names[MAX_DIMENSIONS];
     char                *out_dim_names[MAX_DIMENSIONS];
-    Real                amount_done, value;
+    Real                amount_done;
     Real                real_min, real_max;
     STRING              dim_names[MAX_DIMENSIONS];
     Volume              volume;
     General_transform   transform;
     Minc_file           in_file, out_file;
-    minc_output_options options, *options_ptr;
+    minc_output_options options;
     volume_input_struct volume_input;
     nc_type             data_type;
     BOOLEAN             signed_flag;
@@ -66,8 +65,6 @@
     {
         /*--- determine the output order */
 
-        options_ptr = &options;
-
         for_less( i, 0, N_DIMENSIONS )
         {
             string = argv[a];
@@ -87,6 +84,8 @@
         }
     }
 
+    /*--- open and close volume as 3D to get correct transform */
+
     (void) start_volume_input( input_filename, N_DIMENSIONS,
                                out_dim_names, NC_UNSPECIFIED, FALSE, 0.0, 0.0,
                                TRUE, &volume, (minc_input_options *) NULL,

--- a/volume_io/Volumes/output_mnc.c
+++ b/volume_io/Volumes/output_mnc.c
@@ -2,7 +2,7 @@
 #include  <internal_volume_io.h>
 
 #ifndef lint
-static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/output_mnc.c,v 1.21 1995-02-13 13:14:27 david Exp $";
+static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/output_mnc.c,v 1.22 1995-02-20 13:12:24 david Exp $";
 #endif
 
 #define  INVALID_AXIS   -1
@@ -114,7 +114,7 @@
     Vector              axes[MAX_DIMENSIONS];
     static  STRING      default_dim_names[] = { MIzspace, MIyspace, MIxspace };
     char                *file_dim_names[MAX_VAR_DIMS];
-    Transform           transform;
+    Transform           transform, t, inverse;
     char                **vol_dimension_names;
     minc_output_options default_options;
 
@@ -264,9 +264,16 @@
         }
     }
 
-    start[X] = DOT_POINT_VECTOR( origin, axes[X] );
-    start[Y] = DOT_POINT_VECTOR( origin, axes[Y] );
-    start[Z] = DOT_POINT_VECTOR( origin, axes[Z] );
+    make_identity_transform( &t );
+    set_transform_x_axis( &t, &axes[X] );
+    set_transform_y_axis( &t, &axes[Y] );
+    set_transform_z_axis( &t, &axes[Z] );
+
+    (void) compute_transform_inverse( &t, &inverse );
+
+    transform_point( &inverse,
+                     Point_x(origin), Point_y(origin), Point_z(origin),
+                     &start[X], &start[Y], &start[Z] );
 
     for_less( d, 0, n_dimensions )
     {

--- a/volume_io/Volumes/volumes.c
+++ b/volume_io/Volumes/volumes.c
@@ -3,7 +3,7 @@
 #include  <float.h>
 
 #ifndef lint
-static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/volumes.c,v 1.38 1995-02-17 09:46:24 david Exp $";
+static char rcsid[] = "$Header: /private-cvsroot/minc/volume_io/Volumes/volumes.c,v 1.39 1995-02-20 13:12:21 david Exp $";
 #endif
 
 char   *XYZ_dimension_names[] = { MIxspace, MIyspace, MIzspace };
@@ -997,6 +997,42 @@
     recompute_world_transform( volume );
 }
 
+private  void  reorder_voxel_to_xyz(
+    Volume   volume,
+    Real     voxel[],
+    Real     xyz[] )
+{
+    int   c, axis;
+
+    for_less( c, 0, N_DIMENSIONS )
+    {
+        axis = volume->spatial_axes[c];
+        if( axis >= 0 )
+            xyz[c] = voxel[axis];
+        else
+            xyz[c] = 0.0;
+    }
+}
+
+private  void  reorder_xyz_to_voxel(
+    Volume   volume,
+    Real     xyz[],
+    Real     voxel[] )
+{
+    int   c, axis, n_dims;
+
+    n_dims = get_volume_n_dimensions( volume );
+    for_less( c, 0, n_dims )
+        voxel[c] = 0.0;
+
+    for_less( c, 0, N_DIMENSIONS )
+    {
+        axis = volume->spatial_axes[c];
+        if( axis >= 0 )
+            voxel[axis] = xyz[c];
+    }
+}
+
 /* ----------------------------- MNI Header -----------------------------------
 @NAME       : convert_voxel_to_world
 @INPUT      : volume
@@ -1021,27 +1057,14 @@
     Real     *y_world,
     Real     *z_world )
 {
-    Real   x_voxel, y_voxel, z_voxel;
-
-    if( volume->spatial_axes[0] >= 0 )
-        x_voxel = voxel[volume->spatial_axes[0]];
-    else
-        x_voxel = 0.0;
+    Real   xyz[N_DIMENSIONS];
 
-    if( volume->spatial_axes[1] >= 0 )
-        y_voxel = voxel[volume->spatial_axes[1]];
-    else
-        y_voxel = 0.0;
-
-    if( volume->spatial_axes[2] >= 0 )
-        z_voxel = voxel[volume->spatial_axes[2]];
-    else
-        z_voxel = 0.0;
+    reorder_voxel_to_xyz( volume, voxel, xyz );
 
     /* apply linear transform */
 
     general_transform_point( &volume->voxel_to_world_transform,
-                             x_voxel, y_voxel, z_voxel,
+                             xyz[X], xyz[Y], xyz[Z],
                              x_world, y_world, z_world );
 }
 
@@ -1090,9 +1113,9 @@
 
 /* ----------------------------- MNI Header -----------------------------------@NAME       : convert_voxel_normal_vector_to_world
 @INPUT      : volume
-              x_voxel
-              y_voxel
-              z_voxel
+              voxel_vector0
+              voxel_vector1
+              voxel_vector2
 @OUTPUT     : x_world
               y_world
               z_world
@@ -1105,15 +1128,22 @@
 ---------------------------------------------------------------------------- */
 public  void  convert_voxel_normal_vector_to_world(
     Volume          volume,
-    Real            x_voxel,
-    Real            y_voxel,
-    Real            z_voxel,
+    Real            voxel_vector0,
+    Real            voxel_vector1,
+    Real            voxel_vector2,
     Real            *x_world,
     Real            *y_world,
     Real            *z_world )
 {
+    Real        voxel[N_DIMENSIONS], xyz[N_DIMENSIONS];
     Transform   *inverse;
 
+    voxel[0] = voxel_vector0;
+    voxel[1] = voxel_vector1;
+    voxel[2] = voxel_vector2;
+
+    reorder_voxel_to_xyz( volume, voxel, xyz );
+
     if( get_transform_type( &volume->voxel_to_world_transform ) != LINEAR )
         handle_internal_error( "Cannot get normal vector of nonlinear xforms.");
 
@@ -1122,15 +1152,15 @@
 
     /* transform vector by transpose of inverse transformation */
 
-    *x_world = Transform_elem(*inverse,0,0) * x_voxel +
-               Transform_elem(*inverse,1,0) * y_voxel +
-               Transform_elem(*inverse,2,0) * z_voxel;
-    *y_world = Transform_elem(*inverse,0,1) * x_voxel +
-               Transform_elem(*inverse,1,1) * y_voxel +
-               Transform_elem(*inverse,2,1) * z_voxel;
-    *z_world = Transform_elem(*inverse,0,2) * x_voxel +
-               Transform_elem(*inverse,1,2) * y_voxel +
-               Transform_elem(*inverse,2,2) * z_voxel;
+    *x_world = Transform_elem(*inverse,0,0) * xyz[X] +
+               Transform_elem(*inverse,1,0) * xyz[Y] +
+               Transform_elem(*inverse,2,0) * xyz[Z];
+    *y_world = Transform_elem(*inverse,0,1) * xyz[X] +
+               Transform_elem(*inverse,1,1) * xyz[Y] +
+               Transform_elem(*inverse,2,1) * xyz[Z];
+    *z_world = Transform_elem(*inverse,0,2) * xyz[X] +
+               Transform_elem(*inverse,1,2) * xyz[Y] +
+               Transform_elem(*inverse,2,2) * xyz[Z];
 }
 
 /* ----------------------------- MNI Header -----------------------------------
@@ -1220,26 +1250,15 @@
     Real     z_world,
     Real     voxel[] )
 {
-    int    c;
-    Real   x_voxel, y_voxel, z_voxel;
+    Real   xyz[N_DIMENSIONS];
 
     /* apply linear transform */
 
     general_inverse_transform_point( &volume->voxel_to_world_transform,
-                     x_world, y_world, z_world,
-                     &x_voxel, &y_voxel, &z_voxel );
-
-    for_less( c, 0, get_volume_n_dimensions(volume) )
-        voxel[c] = 0.0;
+                                     x_world, y_world, z_world,
+                                     &xyz[X], &xyz[Y], &xyz[Z] );
 
-    if( volume->spatial_axes[0] >= 0 )
-        voxel[volume->spatial_axes[0]] = x_voxel;
-
-    if( volume->spatial_axes[1] >= 0 )
-        voxel[volume->spatial_axes[1]] = y_voxel;
-
-    if( volume->spatial_axes[2] >= 0 )
-        voxel[volume->spatial_axes[2]] = z_voxel;
+    reorder_xyz_to_voxel( volume, xyz, voxel );
 }
 
 /* ----------------------------- MNI Header -----------------------------------