changeset 1619:00e9e818f8c9

Streamline, use hypercube functions where appropriate
author bert <bert>
date Thu, 08 Jan 2004 20:57:10 +0000
parents 8c06c4fbc183
children 453ef88a2f6b
files libsrc2/convert.c
diffstat 1 files changed, 55 insertions(+), 384 deletions(-) [+]
line wrap: on
line diff
--- a/libsrc2/convert.c
+++ b/libsrc2/convert.c
@@ -92,305 +92,6 @@
     return (result);
 }
 
-/* Stuff from volume_io 
- */
-
-/**  Computes the cross product of 2 3-dimensional vectors.
- */
-static void
-micross_3D_vector(const double v1[MI2_3D], /**< Input vector #1 */
-                  const double v2[MI2_3D], /**< Input vector #2 */
-                  double cross[MI2_3D]) /**< Output vector */
-{
-    cross[0] = v1[1] * v2[2] - v1[2] * v2[1];
-    cross[1] = v1[2] * v2[0] - v1[0] * v2[2];
-    cross[2] = v1[0] * v2[1] - v1[1] * v2[0];
-}
-
-/** 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;
-    }
-}
-
-/** Computes a linear transform from the indices of the spatial
- dimensions, the step sizes, the start offsets, and the direction
- cosines.
-*/
-static void 
-micompute_v2w_transform(double steps[MI2_3D],
-                        double cosines[MI2_3D][MI2_3D],
-                        double starts[MI2_3D],
-                        mi_lin_xfm_t transform)
-{
-    /* These buffers are here so that we can reorder these if need be.
-     * For now, they really serve no purpose.
-     */
-    double steps_3D[MI2_3D];
-    double cosines_3D[MI2_3D][MI2_3D];
-    double starts_3D[MI2_3D];
-    int i, j;
-
-    /* Find how many direction cosines are specified, and set the
-     * 3D steps and starts.
-     * TODO: THis is where we would reorder values into XYZ order if
-     * necessary!
-     */
-    for ( i = 0; i < MI2_3D; i++ ) {
-        steps_3D[i] = steps[i];
-        starts_3D[i] = starts[i];
-        cosines_3D[i][MI2_X] = cosines[i][MI2_X];
-        cosines_3D[i][MI2_Y] = cosines[i][MI2_Y];
-        cosines_3D[i][MI2_Z] = cosines[i][MI2_Z];
-    }
-
-    /* make sure the 3 axes are not a singular system */
-
-    for ( i = 0; i < MI2_3D; i++ ) {
-        double normal[MI2_3D];
-
-        micross_3D_vector( cosines_3D[i], cosines_3D[(i + 1) % MI2_3D],
-                           normal );
-        if ( normal[0] == 0.0 && normal[1] == 0.0 && normal[2] == 0.0 ) {
-            break;              /* Singular */
-        }
-    }
-
-    if ( i < MI2_3D ) {
-        /* Singular system, convert to identity... 
-         */
-        cosines_3D[0][0] = 1.0;
-        cosines_3D[0][1] = 0.0;
-        cosines_3D[0][2] = 0.0;
-        cosines_3D[1][0] = 0.0;
-        cosines_3D[1][1] = 1.0;
-        cosines_3D[1][2] = 0.0;
-        cosines_3D[2][0] = 0.0;
-        cosines_3D[2][1] = 0.0;
-        cosines_3D[2][2] = 1.0;
-    }
-
-    /* Make the linear transformation 
-     */
-    mimake_identity_transform( transform );
-
-    /* Calculate the voxel-to-world transform 
-     */
-    for ( i = 0; i < MI2_3D; i++ ) {
-        for ( j = 0; j < MI2_3D; j++ ) {
-            transform[j][i] = cosines_3D[i][j] * steps_3D[i];
-            transform[j][3] += cosines_3D[i][j] * starts_3D[i];
-        }
-    }
-}
-
-/** This is just a constant which implies that the transform is a homogenous
- * linear transform.
- */
-#define W 1.0
-
-/** Transforms the point \a in[3] by the homogenous transform
-    matrix, resulting in \a out[3]
-*/
-static void 
-mixfm_point(mi_lin_xfm_t transform,
-            const double in[MI2_3D],
-            double out[MI2_3D])
-{
-    double tmp[MI2_LIN_XFM_SIZE];
-    int i;
-
-    for (i = 0; i < MI2_LIN_XFM_SIZE; i++) {
-      tmp[i] = (transform[i][0] * in[MI2_X] +
-                transform[i][1] * in[MI2_Y] +
-                transform[i][2] * in[MI2_Z] +
-                transform[i][3] * W);
-    }
-
-    if ( tmp[3] != 0.0 && tmp[3] != 1.0 ) {
-        for (i = 0; i < 3; i++) {
-            tmp[i] /= tmp[3];
-        }
-    }
-
-    for (i = 0; i < MI2_3D; i++) {
-        out[i] = tmp[i];
-    }
-}
-
-/** Performs scaled maximal pivoting gaussian elimination as a
-    numerically robust method to solve systems of linear equations.
-*/
-
-static int
-scaled_maximal_pivoting_gaussian_elimination(int   row[4],
-                                             double  a[4][4],
-                                             double solution[4][4] )
-{
-    int i, j, k, p, v, tmp;
-    double s[4], val, best_val, m, scale_factor;
-    int result;
-
-    for ( i = 0; i < 4; i++) {
-        row[i] = i;
-    }
-
-    for ( i = 0; i < 4; i++) {
-        s[i] = fabs( a[i][0] );
-        for ( j = 1; j < 4; j++ ) {
-            if ( fabs(a[i][j]) > s[i] ) {
-                s[i] = fabs(a[i][j]);
-            }
-        }
-
-        if ( s[i] == 0.0 ) {
-            return ( MI_ERROR );
-        }
-    }
-
-    result = MI_NOERROR;
-
-    for ( i = 0; i < 4 - 1; i++ ) {
-        p = i;
-        best_val = a[row[i]][i] / s[row[i]];
-        best_val = fabs( best_val );
-        for ( j = i + 1; j < 4; 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 ) {
-            result = MI_ERROR;
-            break;
-        }
-
-        if ( i != p ) {
-            tmp = row[i];
-            row[i] = row[p];
-            row[p] = tmp;
-        }
-
-        for ( j = i + 1; j < 4; j++ ) {
-            if ( a[row[i]][i] == 0.0 ) {
-                result = MI_ERROR;
-                break;
-            }
-
-            m = a[row[j]][i] / a[row[i]][i];
-            for ( k = i + 1; k < 4; k++ )
-                a[row[j]][k] -= m * a[row[i]][k];
-            for( v = 0; v < 4; v++ )
-                solution[row[j]][v] -= m * solution[row[i]][v];
-        }
-
-        if (result != MI_NOERROR)
-            break;
-    }
-
-    if ( result == MI_NOERROR && a[row[4-1]][4-1] == 0.0 )
-        result = MI_ERROR;
-
-    if ( result == MI_NOERROR ) {
-        for ( i = 4-1;  i >= 0;  --i ) {
-            for ( j = i+1; j < 4; j++) {
-                scale_factor = a[row[i]][j];
-                for ( v = 0; v < 4; v++ )
-                    solution[row[i]][v] -= scale_factor * solution[row[j]][v];
-            }
-
-            for( v = 0; v < 4; v++ )
-                solution[row[i]][v] /= a[row[i]][i];
-        }
-    }
-
-    return (result);
-}
-
-static int
-scaled_maximal_pivoting_gaussian_elimination_real(double coefs[4][4],
-                                                  double values[4][4])
-{
-    int i, j, v, row[4];
-    double a[4][4], solution[4][4];
-    int result;
-
-    for ( i = 0; i < 4; i++) {
-        for (j = 0; j < 4; j++)
-            a[i][j] = coefs[i][j];
-        for ( v = 0; v < 4; v++ )
-            solution[i][v] = values[v][i];
-    }
-
-    result = scaled_maximal_pivoting_gaussian_elimination(row, a, solution);
-
-    if ( result == MI_NOERROR ) {
-        for ( i = 0; i < 4; i++ ) {
-            for ( v = 0; v < 4; v++ )
-                values[v][i] = solution[row[i]][v];
-        }
-    }
-
-    return ( result );
-}
-
-/** 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;
-
-    /* Start off with the identity matrix. */
-    for ( i = 0; i < 4; i++ ) {
-        for ( j = 0; j < 4; j++ ) {
-            inverse[i][j] = 0.0;
-        }
-        inverse[i][i] = 1.0;
-    }
-
-    result = scaled_maximal_pivoting_gaussian_elimination_real( matrix,
-                                                                inverse );
-
-    if ( result == MI_NOERROR )  {
-        for ( i = 0; i < 4 - 1; i++) {
-            for ( j = i + 1; j < 4; j++ ) {
-                tmp = inverse[i][j];
-                inverse[i][j] = inverse[j][i];
-                inverse[j][i] = tmp;
-            }
-        }
-    }
-    return (result);
-}
-
-static 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 );
-}
-
 /** Converts a 3-dimensional position in voxel coordinates into a 
  * 3-dimensional position in world coordinates.
  */
@@ -399,27 +100,8 @@
                             const double voxel[MI2_3D],
                             double world[MI2_3D])
 {
-#if 0                           /* DMcD's way */
-    int result = MI_NOERROR;    /* Return code */
-    hid_t file_id;              /* HDF5 file ID */
-    mi_lin_xfm_t v2w_transform; /* Voxel-to-world transform */
-
-    file_id = volume->hdf_id;
-
-    miget_voxel_to_world(file_id, v2w_transform);
-
-    mixfm_point(v2w_transform, voxel, world);
-
-#else  /* Peter's way */
-    int result = MI_NOERROR;
-    hid_t file_id;
-    mi_lin_xfm_t v2w_transform;
-
-    file_id = volume->hdf_id;
-    miget_voxel_to_world(file_id, v2w_transform);
-    mitransform_coord(world, v2w_transform, voxel);
-#endif
-    return (result);
+    mitransform_coord(world, volume->v2w_transform, voxel);
+    return (MI_NOERROR);
 }
 
 /** Converts a 3-dimensional position in world coordinates into a 
@@ -429,39 +111,8 @@
                                 const double world[MI2_3D],
                                 double voxel[MI2_3D])
 {
-    int result = MI_NOERROR;    /* Return code */
-    hid_t file_id;              /* HDF5 file ID */
-    mi_lin_xfm_t v2w_transform; /* Voxel-to-world transform */
-    mi_lin_xfm_t w2v_transform; /* World-to-voxel transform */
-
-    file_id = volume->hdf_id;
-
-    miget_voxel_to_world(file_id, v2w_transform );
-
-    miinvert_transform(v2w_transform, w2v_transform);
-
-    { 
-        int i, j;
-        printf("Voxel to World:\n");
-        for (i = 0; i < 4; i++) {
-            for (j = 0; j < 4; j++) {
-                printf("%f ", v2w_transform[i][j]);
-            }
-            printf("\n");
-        }
-
-        printf("World to Voxel:\n");
-        for (i = 0; i < 4; i++) {
-            for (j = 0; j < 4; j++) {
-                printf("%f ", w2v_transform[i][j]);
-            }
-            printf("\n");
-        }
-    }
-
-    mixfm_point(w2v_transform, world, voxel);
-
-    return (result);
+    mitransform_coord(voxel, volume->w2v_transform, world);
+    return (MI_NOERROR);
 }
 
 int 
@@ -488,15 +139,7 @@
                          const double voxel[],
                          double world[])
 {
-    int result = MI_NOERROR;    /* Return code */
-    hid_t file_id;              /* HDF5 file ID */
-    int length;                 /* Attribute length */
-    mi_lin_xfm_t v2w_transform; /* Voxel-to-world transform */
-    double cosines[MI2_3D][MI2_3D]; /* Direction cosines */
-    double starts[MI2_3D];      /* Start positions */
-    double steps[MI2_3D];       /* Step sizes */
-
-    return (result);
+    return (MI_NOERROR);
 }
 
 int
@@ -593,16 +236,16 @@
                   int ndims,
                   double *voxel_ptr)
 {
-    hid_t file_id;
     int result;
+    unsigned long count[MAX_VAR_DIMS];
+    int i;
 
-    file_id = volume->hdf_id;
-    result = mivarget1(file_id, /* The file handle */
-                       //MI2varid(file_id, MIimage), /* The variable ID */
-                       (long *) coords, /* The voxel coordinates */
-                       NC_DOUBLE, /* The datatype */
-                       MI_SIGNED, /* Ignored for floating point data */
-                       voxel_ptr); /* Location to store result */
+    for (i = 0; i < ndims; i++) {
+        count[i] = 1;
+    }
+    
+    result = miget_voxel_value_hyperslab(volume, MI_TYPE_DOUBLE, 
+                                         coords, count, voxel_ptr);
     return (result);
 }
 
@@ -616,16 +259,16 @@
                   int ndims,
                   double voxel)
 {
-    hid_t file_id;
     int result;
+    unsigned long count[MAX_VAR_DIMS];
+    int i;
 
-    file_id = volume->hdf_id;
-    result = mivarput1(file_id, /* The file handle */
-                       //MI2varid(file_id, MIimage), /* The variable ID */
-                       (long *) coords, /* The voxel coordinates */
-                       NC_DOUBLE, /* The datatype */
-                       MI_SIGNED, /* Ignored for floating point data */
-                       &voxel); /* Location from which to store result */
+    for (i = 0; i < ndims; i++) {
+        count[i] = 1;
+    }
+    
+    result = miset_voxel_value_hyperslab(volume, MI_TYPE_DOUBLE, 
+                                         coords, count, &voxel);
     return (result);
 }
 
@@ -651,7 +294,15 @@
     double r1, r2;
     double f;
 
+    if (argc != 5) {
+        TESTRPT("must specify a file name and position!", 0);
+        exit(-1);
+    }
+
     result = miopen_volume(argv[1], MI2_OPEN_READ, &hvol);
+    if (result < 0) {
+        TESTRPT("miopen_volume error", result);
+    }
 
     coords[0] = atof(argv[2]);
     coords[1] = atof(argv[3]);
@@ -677,32 +328,52 @@
 
     /* Get a voxel value.
      */
-    miget_voxel_value_hyperslab(hvol, MI_TYPE_DOUBLE,
-                                coords, count, &v1);
+    result = miget_voxel_value_hyperslab(hvol, MI_TYPE_DOUBLE,
+                                         coords, count, &v1);
+    if (result < 0) {
+        TESTRPT("miget_voxel_value_hyperslab error", result);
+    }
 
     /* Convert from voxel to real.
      */
-    miconvert_voxel_to_real(hvol, coords, 3, v1, &r1);
+    result = miconvert_voxel_to_real(hvol, coords, 3, v1, &r1);
+    if (result < 0) {
+        TESTRPT("miconvert_voxel_to_real error", result);
+    }
     printf("voxel %f => real %f\n", v1, r1);
 
     /* Convert it back to voxel.
      */
-    miconvert_real_to_voxel(hvol, coords, 3, r1, &v2);
+    result = miconvert_real_to_voxel(hvol, coords, 3, r1, &v2);
+    if (result < 0) {
+        TESTRPT("miconvert_real_to_voxel error", result);
+    }
     printf("real %f => voxel %f\n", r1, v2);
 
     /* Compare to the value read via the ICV
      */
-    miget_real_value(hvol, coords, 3, &r2);
-
+    result = miget_real_value(hvol, coords, 3, &r2);
+    if (result < 0) {
+        TESTRPT("miget_real_value error", result);
+    }
     printf("real from ICV: %f\n", r2);
     printf("\n");
 
+
     if (v1 != v2) {
         TESTRPT("Voxel value mismatch", 0);
     }
     if (r1 != r2) {
         TESTRPT("Real value mismatch", 0);
     }
+
+    result = miget_voxel_value(hvol, coords, 3, &v1);
+    if (result < 0) {
+        TESTRPT("miget_voxel_value error", result);
+    }
+
+    printf("voxel from mivarget1: %f\n", v1);
+
     return (error_cnt);
 }