Home My Page Projects Code Snippets Project Openings diderot
Summary Activity Tracker Tasks SCM

SCM Repository

[diderot] Diff of /trunk/src/lib/common/image.c
ViewVC logotype

Diff of /trunk/src/lib/common/image.c

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1231, Mon May 16 13:49:17 2011 UTC revision 1232, Mon May 16 23:37:52 2011 UTC
# Line 1  Line 1 
1  /*! \file image.c  /*! \file image.c
2   *   *
3   * \author John Reppy   * \author John Reppy
4     *
5     * \brief code to support the loading of image data from Nrrd files.
6   */   */
7    
8  /*  /*
9   * COPYRIGHT (c) 2010 The Diderot Project (http://diderot-language.cs.uchicago.edu)   * COPYRIGHT (c) 2011 The Diderot Project (http://diderot-language.cs.uchicago.edu)
10   * All rights reserved.   * All rights reserved.
11   */   */
12    
13  #include "Diderot/diderot.h"  #include "Diderot/diderot.h"
14  #include <teem/nrrd.h>  #include <teem/nrrd.h>
15    
16    /* FIXME: eventally we should change the naming conventions in the header files to be generic */
17    #if defined(DIDEROT_SINGLE_PRECISION)
18    #define vec2(a,b)               vec2f(a,b)
19    #define vec3(a,b,c)             vec3f(a,b,c)
20    #define mulMat2x2Vec2(m,v)      mulMat2x2Vec2f(m,v)
21    #define mulMat3x3Vec3(m,v)      mulMat3x3Vec3f(m,v)
22    #else
23    #define vec2(a,b)               vec2d(a,b)
24    #define vec3(a,b,c)             vec3d(a,b,c)
25    #define mulMat2x2Vec2(m,v)      mulMat2x2Vec2d(m,v)
26    #define mulMat3x3Vec3(m,v)      mulMat3x3Vec3d(m,v)
27    #endif
28    
29    // we use double precision representations of the transforms here to
30    // improve robustness, even when the results are stored as single-precision
31    typedef double Matrix2x2_t[4];
32    typedef double Matrix3x3_t[9];
33    typedef double Matrix4x4_t[16];
34    
35    typedef double Vec2_t[2];
36    typedef double Vec3_t[3];
37    
38  static void PrintMat2x2 (const char *title, Diderot_Mat2x2_t m)  static void PrintMat2x2 (const char *title, Matrix2x2_t m)
39  {  {
40  #define M(i)    m[(i)/2].r[(i)%2]  #define M(i)    m[i]
41      printf ("%s @ %p:\n", title, (void *)m);      printf ("%s @ %p:\n", title, (void *)m);
42      for (int i = 0;  i < 2;  i++) {      for (int i = 0;  i < 2;  i++) {
43          int j = 2*i;          int j = 2*i;
44          printf ("  [ %6f  %6f ]\n", M(j+0), M(j+1));          printf ("  [ %6lf  %6lf ]\n", M(j+0), M(j+1));
45      }      }
46  #undef M  #undef M
47    
48  }  }
49    
50  static void PrintMat3x3 (const char *title, Diderot_Mat3x3_t m)  static void PrintMat3x3 (const char *title, Matrix3x3_t m)
51  {  {
52  #define M(i)    m[(i)/3].r[(i)%3]  #define M(i)    m[i]
53      printf ("%s @ %p:\n", title, (void *)m);      printf ("%s @ %p:\n", title, (void *)m);
54      for (int i = 0;  i < 3;  i++) {      for (int i = 0;  i < 3;  i++) {
55          int j = 3*i;          int j = 3*i;
56          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));
57      }      }
58  #undef M  #undef M
59    
60  }  }
61    
62  static void PrintMat4x4 (const char *title, Diderot_Mat4x4_t m)  static void PrintMat4x4 (const char *title, Matrix4x4_t m)
63  {  {
64  #define M(i)    m[(i)>>2].r[(i)&3]  #define M(i)    m[i]
65      printf ("%s @ %p:\n", title, (void *)m);      printf ("%s @ %p:\n", title, (void *)m);
66      for (int i = 0;  i < 4;  i++) {      for (int i = 0;  i < 4;  i++) {
67          int j = 4*i;          int j = 4*i;
68          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));
69      }      }
70  #undef M  #undef M
71    
# Line 50  Line 73 
73    
74  /* Loading transfomation information from Nrrd files */  /* Loading transfomation information from Nrrd files */
75    
76  static void LoadTransform1D (Nrrd *nin, float *s, float *t)  static void LoadTransform1D (Nrrd *nin, double *s, double *t)
77  {  {
 #define M(i,j)  m[i].r[j]  
   
78    // compute the offset to the first space axis    // compute the offset to the first space axis
79      int base = nin->dim - nin->spaceDim;      int base = nin->dim - nin->spaceDim;
80    
# Line 62  Line 83 
83    
84    // Image location    // Image location
85      *t = nin->spaceOrigin[0];      *t = nin->spaceOrigin[0];
   
 #undef M  
86  }  }
87    
88  static void LoadTransform2D (Nrrd *nin, Diderot_Mat2x2_t m, vec2f_t *t)  static void LoadTransform2D (Nrrd *nin, Matrix2x2_t m, Vec2_t t)
89  {  {
90  #define M(i,j)  m[i].r[j]  #define M(i,j)  m[(i)*2+(j)]
91    
92    // compute the offset to the first space axis    // compute the offset to the first space axis
93      int base = nin->dim - nin->spaceDim;      int base = nin->dim - nin->spaceDim;
# Line 80  Line 99 
99      M(1,1) = nin->axis[base+1].spaceDirection[1];      M(1,1) = nin->axis[base+1].spaceDirection[1];
100    
101    // Image location    // Image location
102      *t = vec2f(      t[0] = nin->spaceOrigin[0];
103          nin->spaceOrigin[0],      t[1] = nin->spaceOrigin[1];
         nin->spaceOrigin[1]);  
104    
105  #undef M  #undef M
106  }  }
107    
108  static void LoadTransform3D (Nrrd *nin, Diderot_Mat3x3_t m, vec3f_t *t)  static void LoadTransform3D (Nrrd *nin, Matrix3x3_t m, Vec3_t t)
109  {  {
110  #define M(i,j)  m[i].r[j]  #define M(i,j)  m[3*(i)+(j)]
111    
112    // compute the offset to the first space axis    // compute the offset to the first space axis
113      int base = nin->dim - nin->spaceDim;      int base = nin->dim - nin->spaceDim;
# Line 106  Line 124 
124      M(2,2) = nin->axis[base+2].spaceDirection[2];      M(2,2) = nin->axis[base+2].spaceDirection[2];
125    
126    // Image location    // Image location
127      *t = vec3f(      t[0] = nin->spaceOrigin[0];
128          nin->spaceOrigin[0],      t[1] = nin->spaceOrigin[1];
129          nin->spaceOrigin[1],      t[2] = nin->spaceOrigin[2];
         nin->spaceOrigin[2]);  
130    
131  #undef M  #undef M
132  }  }
# Line 117  Line 134 
134    
135  /* Transformation matrix operations */  /* Transformation matrix operations */
136    
137  STATIC_INLINE Diderot_real_t Det2 (  STATIC_INLINE double Det2 (
138      Diderot_real_t a, Diderot_real_t b,      double a, double b,
139      Diderot_real_t c, Diderot_real_t d)      double c, double d)
140  {  {
141      return (a*d - b*c);      return (a*d - b*c);
142  }  }
143    
144  STATIC_INLINE Diderot_real_t Det3 (  STATIC_INLINE double Det3 (
145      Diderot_real_t a, Diderot_real_t b, Diderot_real_t c,      double a, double b, double c,
146      Diderot_real_t d, Diderot_real_t e, Diderot_real_t f,      double d, double e, double f,
147      Diderot_real_t g, Diderot_real_t h, Diderot_real_t i)      double g, double h, double i)
148  {  {
149      return (a*e*i      return (a*e*i
150            + d*h*c            + d*h*c
# Line 137  Line 154 
154            - a*h*f);            - a*h*f);
155  }  }
156    
157  static Diderot_real_t DetM3x3 (Diderot_Mat3x3_t m)  static double DetM3x3 (Matrix3x3_t m)
158  {  {
159  #define M(i)    m[(i)/3].r[(i)%3]  #define M(i)    m[i]
160      return Det3(M(0), M(1), M(2),      return Det3(M(0), M(1), M(2),
161                  M(3), M(4), M(5),                  M(3), M(4), M(5),
162                  M(6), M(7), M(8));                  M(6), M(7), M(8));
163  #undef M  #undef M
164  }  }
165    
166  static Diderot_real_t DetM4x4 (Diderot_Mat4x4_t m)  static double DetM4x4 (Matrix4x4_t m)
167  {  {
168  #define M(i)    m[(i)>>2].r[(i)&3]  #define M(i)    m[i]
169          return  (M( 0) * Det3(M( 5), M( 6), M( 7),          return  (M( 0) * Det3(M( 5), M( 6), M( 7),
170                                M( 9), M(10), M(11),                                M( 9), M(10), M(11),
171                                M(13), M(14), M(15))                                M(13), M(14), M(15))
# Line 169  Line 186 
186   *  \param m the matrix to invert   *  \param m the matrix to invert
187   *  \param i the inverted matrix   *  \param i the inverted matrix
188   */   */
189  static void InvertM2x2 (Diderot_Mat2x2_t i, Diderot_Mat2x2_t m)  static void InvertM2x2 (Matrix2x2_t i, Matrix2x2_t m)
190  {  {
191  #define M(i)    m[(i)/2].r[(i)%2]  #define M(i)    m[i]
192  #define I(j)    i[(j)/2].r[(j)%2]  #define I(j)    i[j]
193          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));
194    
195          I(0) =  scale * M(3);          I(0) =  scale * M(3);
196          I(1) = -scale * M(1);          I(1) = -scale * M(1);
# Line 187  Line 204 
204   *  \param m the matrix to invert   *  \param m the matrix to invert
205   *  \param i the inverted matrix   *  \param i the inverted matrix
206   */   */
207  static void InvertM3x3 (Diderot_Mat3x3_t i, Diderot_Mat3x3_t m)  static void InvertM3x3 (Matrix3x3_t i, Matrix3x3_t m)
208  {  {
209          Diderot_real_t scale = 1.0 / DetM3x3(m);          double scale = 1.0 / DetM3x3(m);
210    
211  #define M(i)    m[(i)/3].r[(i)%3]  #define M(i)    m[i]
212  #define I(j)    i[(j)/3].r[(j)%3]  #define I(j)    i[j]
213          I(0) = scale * Det2 (M(4), M(5), M(7), M(8));          I(0) = scale * Det2 (M(4), M(5), M(7), M(8));
214          I(1) = scale * Det2 (M(2), M(1), M(8), M(7));          I(1) = scale * Det2 (M(2), M(1), M(8), M(7));
215          I(2) = scale * Det2 (M(1), M(2), M(4), M(5));          I(2) = scale * Det2 (M(1), M(2), M(4), M(5));
# Line 210  Line 227 
227   *  \param m the matrix to invert   *  \param m the matrix to invert
228   *  \param i the inverted matrix   *  \param i the inverted matrix
229   */   */
230  static void InvertM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m)  static void InvertM4x4 (Matrix4x4_t i, Matrix4x4_t m)
231  {  {
232          Diderot_real_t scale = 1.0 / DetM4x4(m);          double scale = 1.0 / DetM4x4(m);
233    
234  #define M(i)    m[(i)>>2].r[(i)&3]  #define M(i)    m[i]
235  #define I(j)    i[(j)>>2].r[(j)&3]  #define I(j)    i[j]
236          I( 0) =  Det3(M( 5),M( 6),M( 7),          I( 0) =  Det3(M( 5),M( 6),M( 7),
237                        M( 9),M(10),M(11),                        M( 9),M(10),M(11),
238                        M(13),M(14),M(15)) * scale;                        M(13),M(14),M(15)) * scale;
# Line 338  Line 355 
355  name, img->size[base+0]);  name, img->size[base+0]);
356    
357    // from the Nrrd file, we load the scaling and translation    // from the Nrrd file, we load the scaling and translation
358      float s, t;      double s, t;
359      LoadTransform1D (nin, &s, &t);      LoadTransform1D (nin, &s, &t);
360    
361      img->s = 1.0 / s;      img->s = (Diderot_real_t)(1.0 / s);
362      img->t = -t / s;      img->t = (Diderot_real_t)(-t / s);
363    
364      *imgOut = img;      *imgOut = img;
365      return DIDEROT_OK;      return DIDEROT_OK;
# Line 380  Line 397 
397  name, img->size[0], img->size[1]);  name, img->size[0], img->size[1]);
398    
399    // from the Nrrd file, we load the affine image-to-world transform matrix    // from the Nrrd file, we load the affine image-to-world transform matrix
400      Diderot_Mat2x2_t m;         // rotation and scaling      Matrix2x2_t m;              // rotation and scaling
401      vec2f_t t;                  // translation part of the transform      Vec2_t t;                   // translation part of the transform
402      LoadTransform2D (nin, m, &t);      LoadTransform2D (nin, m, t);
403    
404  PrintMat2x2("m", m);  PrintMat2x2("m", m);
405  printf("t = <%f, %f>\n", ((union2f_t)t).r[0], ((union2f_t)t).r[1]);  printf("t = <%f, %f>\n", t[0], t[1]);
406    
407    // compute inverse of m, which is the transform from the world basis to the image basis    // compute inverse of m, which is the transform from the world basis to the image basis
408      InvertM2x2 (img->w2i, m);      Matrix2x2_t mInv;
409      img->tVec = -mulMat2x2Vec2f(img->w2i, t);  // inv([M t]) = [inv(M) -inv(M)t]      InvertM2x2 (mInv, m);
410  PrintMat2x2("w2i", img->w2i);  PrintMat2x2("w2i", mInv);
 printf("tVec = <%f, %f>\n", ((union2f_t)img->tVec).r[0], ((union2f_t)img->tVec).r[1]);  
411    
412    // compute transpose of w2i    // copy results into the image data structure
413        for (int i = 0;  i < 2;  i++) {
414            for (int j = 0;  j < 2;  j++) {
415                img->w2i[i].r[j] = mInv[2*i+j];
416            }
417        }
418      transpose2x2f (img->w2iT, img->w2i);      transpose2x2f (img->w2iT, img->w2i);
419        Diderot_vec2_t tVec = vec2(t[0], t[1]);
420        img->tVec = -mulMat2x2Vec2(img->w2i, tVec);  // inv([M t]) = [inv(M) -inv(M)t]
421    printf("tVec = <%f, %f>\n", ((union2f_t)img->tVec).r[0], ((union2f_t)img->tVec).r[1]);
422    
423      *imgOut = img;      *imgOut = img;
424      return DIDEROT_OK;      return DIDEROT_OK;
# Line 435  Line 459 
459  name, img->size[0], img->size[1], img->size[2]);  name, img->size[0], img->size[1], img->size[2]);
460    
461    // from the Nrrd file, we load the affine image-to-world transform matrix    // from the Nrrd file, we load the affine image-to-world transform matrix
462      Diderot_Mat3x3_t m;         // rotation and scaling      Matrix3x3_t m;              // rotation and scaling
463      vec3f_t t;                  // translation part of the transform      Vec3_t t;                   // translation part of the transform
464      LoadTransform3D (nin, m, &t);      LoadTransform3D (nin, m, t);
465    
466  PrintMat3x3("m", m);  PrintMat3x3("m", m);
467  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]);
468    
469    // compute inverse of m, which is the transform from the world basis to the image basis    // compute inverse of m, which is the transform from the world basis to the image basis
470      InvertM3x3 (img->w2i, m);      Matrix3x3_t mInv;
471      img->tVec = -mulMat3x3Vec3f(img->w2i, t);  // inv([M t]) = [inv(M) -inv(M)t]      InvertM3x3 (mInv, m);
472  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]);  
473    
474    // compute transpose of w2i    // copy results into the image data structure
475        for (int i = 0;  i < 3;  i++) {
476            for (int j = 0;  j < 3;  j++) {
477                img->w2i[i].r[j] = mInv[3*i+j];
478            }
479        }
480      transpose3x3f (img->w2iT, img->w2i);      transpose3x3f (img->w2iT, img->w2i);
481        Diderot_vec3_t tVec = vec3(t[0], t[1], t[2]);
482        img->tVec = -mulMat3x3Vec3(img->w2i, tVec);  // inv([M t]) = [inv(M) -inv(M)t]
483    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]);
484    
485      *imgOut = img;      *imgOut = img;
486      return DIDEROT_OK;      return DIDEROT_OK;
487  }  }
   
 /* functions to get input-parameter values */  
 Status_t Diderot_InputString (const char *name, const char **v, bool hasDflt)  
 {  
     return DIDEROT_FAIL;  
 }  
   
 Status_t Diderot_Inputf (const char *name, float *v, bool hasDflt)  
 {  
     char        buf[256];  
     float       f;  
   
     while (true) {  
         if (hasDflt)  
             printf("Enter value for %s (default %f): ", name, *v);  
         else  
             printf("Enter value for %s: ", name);  
         fflush (stdout);  
   
         char *p = fgets(buf, sizeof(buf), stdin);  
         if (p == NULL)  
             return DIDEROT_FAIL;        // EOF  
         int n = sscanf(buf, "%f\n", &f);  
         if (n == 1) {  
             *v = f;  
             return DIDEROT_OK;;  
         }  
         else if (hasDflt)  
             return DIDEROT_OK;;  
     }  
   
 }  
   
 Status_t Diderot_InputVec3f (const char *name, vec3f_t *v, bool hasDflt)  
 {  
     char        buf[256];  
     float       f1, f2, f3;  
   
     while (true) {  
         if (hasDflt) {  
             union3f_t u;  
             u.v = *v;  
             printf("Enter value for %s (default %f %f %f): ",  
                 name, u.r[0], u.r[1], u.r[2]);  
         }  
         else  
             printf("Enter value for %s: ", name);  
         fflush (stdout);  
   
         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);  
         if (n == 3) {  
             *v = vec3f(f1, f2, f3);  
             return DIDEROT_OK;;  
         }  
         else if (hasDflt)  
             return DIDEROT_OK;;  
     }  
   
 }  

Legend:
Removed from v.1231  
changed lines
  Added in v.1232

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