SCM Repository

[diderot] Diff of /branches/pure-cfg/src/lib/common/image.c

[diderot] / branches / pure-cfg / src / lib / common / image.c

Diff of /branches/pure-cfg/src/lib/common/image.c

-revision 1212, Fri May 13 08:00:33 2011 UTC
+revision 1213, Fri May 13 14:29:16 2011 UTC
 Line 1
  /*! \file image.c
   *
   * \author John Reppy
+  *
+  * \brief code to support the loading of image data from Nrrd files.
   */
  /*
-  * COPYRIGHT (c) 2010 The Diderot Project (http://diderot-language.cs.uchicago.edu)
+  * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
   * All rights reserved.
   */
  #include "Diderot/diderot.h"
  #include <teem/nrrd.h>
+ /* FIXME: eventally we should change the naming conventions in the header files to be generic */
+ #if defined(DIDEROT_SINGLE_PRECISION)
+ #define vec2(a,b)               vec2f(a,b)
+ #define vec3(a,b,c)             vec3f(a,b,c)
+ #define mulMat2x2Vec2(m,v)      mulMat2x2Vec2f(m,v)
+ #define mulMat3x3Vec3(m,v)      mulMat3x3Vec3f(m,v)
+ #else
+ #define vec2(a,b)               vec2d(a,b)
+ #define vec3(a,b,c)             vec3d(a,b,c)
+ #define mulMat2x2Vec2(m,v)      mulMat2x2Vec2d(m,v)
+ #define mulMat3x3Vec3(m,v)      mulMat3x3Vec3d(m,v)
+ #endif
+ // we use double precision representations of the transforms here to
+ // improve robustness, even when the results are stored as single-precision
+ typedef double Matrix2x2_t[4];
+ typedef double Matrix3x3_t[9];
+ typedef double Matrix4x4_t[16];
+ typedef double Vec2_t[2];
+ typedef double Vec3_t[3];
- static void PrintMat2x2 (const char *title, Diderot_Mat2x2_t m)
+ static void PrintMat2x2 (const char *title, Matrix2x2_t m)
  {
- #define M(i)    m[(i)/2].r[(i)%2]
+ #define M(i)    m[i]
      printf ("%s @ %p:\n", title, (void *)m);
      for (int i = 0;  i < 2;  i++) {
          int j = 2*i;
-         printf ("  [ %6f  %6f ]\n", M(j+0), M(j+1));
+         printf ("  [ %6lf  %6lf ]\n", M(j+0), M(j+1));
      }
  #undef M
  }
- static void PrintMat3x3 (const char *title, Diderot_Mat3x3_t m)
+ static void PrintMat3x3 (const char *title, Matrix3x3_t m)
  {
- #define M(i)    m[(i)/3].r[(i)%3]
+ #define M(i)    m[i]
      printf ("%s @ %p:\n", title, (void *)m);
      for (int i = 0;  i < 3;  i++) {
          int j = 3*i;
-         printf ("  [ %6f  %6f  %6f ]\n", M(j+0), M(j+1), M(j+2));
+         printf ("  [ %6lf  %6lf  %6lf ]\n", M(j+0), M(j+1), M(j+2));
      }
  #undef M
  }
- static void PrintMat4x4 (const char *title, Diderot_Mat4x4_t m)
+ static void PrintMat4x4 (const char *title, Matrix4x4_t m)
  {
- #define M(i)    m[(i)>>2].r[(i)&3]
+ #define M(i)    m[i]
      printf ("%s @ %p:\n", title, (void *)m);
      for (int i = 0;  i < 4;  i++) {
          int j = 4*i;
-         printf ("  [ %6f  %6f  %6f  %6f ]\n", M(j+0), M(j+1), M(j+2), M(j+3));
+         printf ("  [ %6lf  %6lf  %6lf  %6lf ]\n", M(j+0), M(j+1), M(j+2), M(j+3));
      }
  #undef M
-Line 50
+Line 73
  /* Loading transfomation information from Nrrd files */
- static void LoadTransform1D (Nrrd *nin, float *s, float *t)
+ static void LoadTransform1D (Nrrd *nin, double *s, double *t)
  {
- #define M(i,j)  m[i].r[j]
    // compute the offset to the first space axis
      int base = nin->dim - nin->spaceDim;
-Line 62
+Line 83
    // Image location
      *t = nin->spaceOrigin[0];
- #undef M
  }
- static void LoadTransform2D (Nrrd *nin, Diderot_Mat2x2_t m, vec2f_t *t)
+ static void LoadTransform2D (Nrrd *nin, Matrix2x2_t m, Vec2_t t)
  {
- #define M(i,j)  m[i].r[j]
+ #define M(i,j)  m[(i)*2+(j)]
    // compute the offset to the first space axis
      int base = nin->dim - nin->spaceDim;
-Line 80
+Line 99
      M(1,1) = nin->axis[base+1].spaceDirection[1];
    // Image location
-     *t = vec2f(
+     t[0] = nin->spaceOrigin[0];
-         nin->spaceOrigin[0],
+     t[1] = nin->spaceOrigin[1];
-         nin->spaceOrigin[1]);
  #undef M
  }
- static void LoadTransform3D (Nrrd *nin, Diderot_Mat3x3_t m, vec3f_t *t)
+ static void LoadTransform3D (Nrrd *nin, Matrix3x3_t m, Vec3_t t)
  {
- #define M(i,j)  m[i].r[j]
+ #define M(i,j)  m[3*(i)+(j)]
    // compute the offset to the first space axis
      int base = nin->dim - nin->spaceDim;
-Line 106
+Line 124
      M(2,2) = nin->axis[base+2].spaceDirection[2];
    // Image location
-     *t = vec3f(
+     t[0] = nin->spaceOrigin[0];
-         nin->spaceOrigin[0],
+     t[1] = nin->spaceOrigin[1];
-         nin->spaceOrigin[1],
+     t[2] = nin->spaceOrigin[2];
-         nin->spaceOrigin[2]);
  #undef M
  }
-Line 117
+Line 134
  /* Transformation matrix operations */
- STATIC_INLINE Diderot_real_t Det2 (
+ STATIC_INLINE double Det2 (
-     Diderot_real_t a, Diderot_real_t b,
+     double a, double b,
-     Diderot_real_t c, Diderot_real_t d)
+     double c, double d)
  {
      return (a*d - b*c);
  }
- STATIC_INLINE Diderot_real_t Det3 (
+ STATIC_INLINE double Det3 (
-     Diderot_real_t a, Diderot_real_t b, Diderot_real_t c,
+     double a, double b, double c,
-     Diderot_real_t d, Diderot_real_t e, Diderot_real_t f,
+     double d, double e, double f,
-     Diderot_real_t g, Diderot_real_t h, Diderot_real_t i)
+     double g, double h, double i)
  {
      return (a*e*i
            + d*h*c
-Line 137
+Line 154
            - a*h*f);
  }
- static Diderot_real_t DetM3x3 (Diderot_Mat3x3_t m)
+ static double DetM3x3 (Matrix3x3_t m)
  {
- #define M(i)    m[(i)/3].r[(i)%3]
+ #define M(i)    m[i]
      return Det3(M(0), M(1), M(2),
                  M(3), M(4), M(5),
                  M(6), M(7), M(8));
  #undef M
  }
- static Diderot_real_t DetM4x4 (Diderot_Mat4x4_t m)
+ static double DetM4x4 (Matrix4x4_t m)
  {
- #define M(i)    m[(i)>>2].r[(i)&3]
+ #define M(i)    m[i]
          return  (M( 0) * Det3(M( 5), M( 6), M( 7),
                                M( 9), M(10), M(11),
                                M(13), M(14), M(15))
-Line 169
+Line 186
   *  \param m the matrix to invert
   *  \param i the inverted matrix
   */
- static void InvertM2x2 (Diderot_Mat2x2_t i, Diderot_Mat2x2_t m)
+ static void InvertM2x2 (Matrix2x2_t i, Matrix2x2_t m)
  {
- #define M(i)    m[(i)/2].r[(i)%2]
+ #define M(i)    m[i]
- #define I(j)    i[(j)/2].r[(j)%2]
+ #define I(j)    i[j]
-         Diderot_real_t scale = 1.0 / Det2(M(0), M(1), M(2), M(3));
+         double scale = 1.0 / Det2(M(0), M(1), M(2), M(3));
          I(0) =  scale * M(3);
          I(1) = -scale * M(1);
-Line 187
+Line 204
   *  \param m the matrix to invert
   *  \param i the inverted matrix
   */
- static void InvertM3x3 (Diderot_Mat3x3_t i, Diderot_Mat3x3_t m)
+ static void InvertM3x3 (Matrix3x3_t i, Matrix3x3_t m)
  {
-         Diderot_real_t scale = 1.0 / DetM3x3(m);
+         double scale = 1.0 / DetM3x3(m);
- #define M(i)    m[(i)/3].r[(i)%3]
+ #define M(i)    m[i]
- #define I(j)    i[(j)/3].r[(j)%3]
+ #define I(j)    i[j]
          I(0) = scale * Det2 (M(4), M(5), M(7), M(8));
          I(1) = scale * Det2 (M(2), M(1), M(8), M(7));
          I(2) = scale * Det2 (M(1), M(2), M(4), M(5));
-Line 210
+Line 227
   *  \param m the matrix to invert
   *  \param i the inverted matrix
   */
- static void InvertM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m)
+ static void InvertM4x4 (Matrix4x4_t i, Matrix4x4_t m)
  {
-         Diderot_real_t scale = 1.0 / DetM4x4(m);
+         double scale = 1.0 / DetM4x4(m);
- #define M(i)    m[(i)>>2].r[(i)&3]
+ #define M(i)    m[i]
- #define I(j)    i[(j)>>2].r[(j)&3]
+ #define I(j)    i[j]
          I( 0) =  Det3(M( 5),M( 6),M( 7),
                        M( 9),M(10),M(11),
                        M(13),M(14),M(15)) * scale;
-Line 338
+Line 355
  name, img->size[base+0]);
    // from the Nrrd file, we load the scaling and translation
-     float s, t;
+     double s, t;
      LoadTransform1D (nin, &s, &t);
-     img->s = 1.0 / s;
+     img->s = (Diderot_real_t)(1.0 / s);
-     img->t = -t / s;
+     img->t = (Diderot_real_t)(-t / s);
      *imgOut = img;
      return DIDEROT_OK;
-Line 380
+Line 397
  name, img->size[0], img->size[1]);
    // from the Nrrd file, we load the affine image-to-world transform matrix
-     Diderot_Mat2x2_t m;         // rotation and scaling
+     Matrix2x2_t m;              // rotation and scaling
-     vec2f_t t;                  // translation part of the transform
+     Vec2_t t;                   // translation part of the transform
-     LoadTransform2D (nin, m, &t);
+     LoadTransform2D (nin, m, t);
  PrintMat2x2("m", m);
- printf("t = <%f, %f>\n", ((union2f_t)t).r[0], ((union2f_t)t).r[1]);
+ printf("t = <%f, %f>\n", t[0], t[1]);
    // compute inverse of m, which is the transform from the world basis to the image basis
-     InvertM2x2 (img->w2i, m);
+     Matrix2x2_t mInv;
-     img->tVec = -mulMat2x2Vec2f(img->w2i, t);  // inv([M t]) = [inv(M) -inv(M)t]
+     InvertM2x2 (mInv, m);
- PrintMat2x2("w2i", img->w2i);
+ PrintMat2x2("w2i", mInv);
- printf("tVec = <%f, %f>\n", ((union2f_t)img->tVec).r[0], ((union2f_t)img->tVec).r[1]);
-   // compute transpose of w2i
+   // copy results into the image data structure
+     for (int i = 0;  i < 2;  i++) {
+         for (int j = 0;  j < 2;  j++) {
+             img->w2i[i].r[j] = mInv[2*i+j];
+         }
+     }
      transpose2x2f (img->w2iT, img->w2i);
+     Diderot_vec2_t tVec = vec2(t[0], t[1]);
+     img->tVec = -mulMat2x2Vec2(img->w2i, tVec);  // inv([M t]) = [inv(M) -inv(M)t]
+ printf("tVec = <%f, %f>\n", ((union2f_t)img->tVec).r[0], ((union2f_t)img->tVec).r[1]);
      *imgOut = img;
      return DIDEROT_OK;
-Line 435
+Line 459
  name, img->size[0], img->size[1], img->size[2]);
    // from the Nrrd file, we load the affine image-to-world transform matrix
-     Diderot_Mat3x3_t m;         // rotation and scaling
+     Matrix3x3_t m;              // rotation and scaling
-     vec3f_t t;                  // translation part of the transform
+     Vec3_t t;                   // translation part of the transform
-     LoadTransform3D (nin, m, &t);
+     LoadTransform3D (nin, m, t);
  PrintMat3x3("m", m);
- printf("t = <%f, %f, %f>\n", ((union3f_t)t).r[0], ((union3f_t)t).r[1], ((union3f_t)t).r[2]);
+ printf("t = <%f, %f, %f>\n", t[0], t[1], t[2]);
    // compute inverse of m, which is the transform from the world basis to the image basis
-     InvertM3x3 (img->w2i, m);
+     Matrix3x3_t mInv;
-     img->tVec = -mulMat3x3Vec3f(img->w2i, t);  // inv([M t]) = [inv(M) -inv(M)t]
+     InvertM3x3 (mInv, m);
- PrintMat3x3("w2i", img->w2i);
+ PrintMat3x3("w2i", mInv);
- printf("tVec = <%f, %f, %f>\n", ((union3f_t)img->tVec).r[0], ((union3f_t)img->tVec).r[1], ((union3f_t)img->tVec).r[2]);
-   // compute transpose of w2i
+   // copy results into the image data structure
+     for (int i = 0;  i < 3;  i++) {
+         for (int j = 0;  j < 3;  j++) {
+             img->w2i[i].r[j] = mInv[3*i+j];
+         }
+     }
      transpose3x3f (img->w2iT, img->w2i);
+     Diderot_vec3_t tVec = vec3(t[0], t[1], t[2]);
+     img->tVec = -mulMat3x3Vec3(img->w2i, tVec);  // inv([M t]) = [inv(M) -inv(M)t]
+ printf("tVec = <%f, %f, %f>\n", ((Diderot_union3_t)img->tVec).r[0], ((Diderot_union3_t)img->tVec).r[1], ((Diderot_union3_t)img->tVec).r[2]);
      *imgOut = img;
      return DIDEROT_OK;
-Line 461
+Line 492
      return DIDEROT_FAIL;
  }
- Status_t Diderot_Inputf (const char *name, float *v, bool hasDflt)
+ Status_t Diderot_InputReal (const char *name, Diderot_real_t *v, bool hasDflt)
  {
      char        buf[256];
-     float       f;
+     double      f;
      while (true) {
          if (hasDflt)
-             printf("Enter value for %s (default %f): ", name, *v);
+             printf("Enter value for %s (default %lf): ", name, (double)*v);
          else
              printf("Enter value for %s: ", name);
          fflush (stdout);
-Line 476
+Line 507
          char *p = fgets(buf, sizeof(buf), stdin);
          if (p == NULL)
              return DIDEROT_FAIL;        // EOF
-         int n = sscanf(buf, "%f\n", &f);
+         int n = sscanf(buf, "%lf\n", &f);
          if (n == 1) {
-             *v = f;
+             *v = (Diderot_real_t)f;
              return DIDEROT_OK;;
          }
          else if (hasDflt)
-Line 487
+Line 518
  }
- Status_t Diderot_InputVec3f (const char *name, vec3f_t *v, bool hasDflt)
+ Status_t Diderot_InputVec3 (const char *name, Diderot_vec3_t *v, bool hasDflt)
  {
      char        buf[256];
-     float       f1, f2, f3;
+     double      f1, f2, f3;
      while (true) {
          if (hasDflt) {
-             union3f_t u;
+             Diderot_union3_t u;
              u.v = *v;
              printf("Enter value for %s (default %f %f %f): ",
-                 name, u.r[0], u.r[1], u.r[2]);
+                 name, (double)(u.r[0]), (double)(u.r[1]), (double)(u.r[2]));
          }
          else
              printf("Enter value for %s: ", name);
-Line 506
+Line 537
          char *p = fgets(buf, sizeof(buf), stdin);
          if (p == NULL)
              return DIDEROT_FAIL;        // EOF
-         int n = sscanf(buf, "%f %f %f\n", &f1, &f2, &f3);
+         int n = sscanf(buf, "%lf %lf %lf\n", &f1, &f2, &f3);
          if (n == 3) {
-             *v = vec3f(f1, f2, f3);
+             *v = vec3((Diderot_real_t)f1, (Diderot_real_t)f2, (Diderot_real_t)f3);
              return DIDEROT_OK;;
          }
          else if (hasDflt)

 Legend:



Removed from v.1212
 


changed lines


 
Added in v.1213
 Legend:



Removed from v.1212
 


changed lines


 
Added in v.1213
-Removed from v.1212
+Added in v.1213

root@smlnj-gforge.cs.uchicago.edu	ViewVC Help
Powered by ViewVC 1.0.0