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

SCM Repository

[diderot] Diff of /branches/pure-cfg/src/lib/diderot.c
ViewVC logotype

Diff of /branches/pure-cfg/src/lib/diderot.c

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

revision 990, Tue Apr 26 23:15:12 2011 UTC revision 991, Wed Apr 27 04:53:08 2011 UTC
# Line 48  Line 48 
48    
49  }  }
50    
51    /* Loading transfomation information from Nrrd files */
52    
53    static void LoadTransform1D (Nrrd *nin, float *s, float *t)
54    {
55    #define M(i,j)  m[i].r[j]
56    
57      // compute the offset to the first space axis
58        int base = nin->dim - nin->spaceDim;
59    
60      // Image axis Scaling and Rotation
61        *s = nin->axis[base+0].spaceDirection[0];
62    
63      // Image location
64        *t = nin->spaceOrigin[0];
65    
66    #undef M
67    }
68    
69    static void LoadTransform2D (Nrrd *nin, Diderot_Mat2x2_t m, vec2f_t *t)
70    {
71    #define M(i,j)  m[i].r[j]
72    
73      // compute the offset to the first space axis
74        int base = nin->dim - nin->spaceDim;
75    
76      // Image axis Scaling and Rotation
77        M(0,0) = nin->axis[base+0].spaceDirection[0];
78        M(1,0) = nin->axis[base+0].spaceDirection[1];
79        M(0,1) = nin->axis[base+1].spaceDirection[0];
80        M(1,1) = nin->axis[base+1].spaceDirection[1];
81    
82      // Image location
83        *t = vec2f(
84            nin->spaceOrigin[0],
85            nin->spaceOrigin[1]);
86    
87    #undef M
88    }
89    
90    static void LoadTransform3D (Nrrd *nin, Diderot_Mat3x3_t m, vec3f_t *t)
91    {
92    #define M(i,j)  m[i].r[j]
93    
94      // compute the offset to the first space axis
95        int base = nin->dim - nin->spaceDim;
96    
97      // Image axis Scaling and Rotation
98        M(0,0) = nin->axis[base+0].spaceDirection[0];
99        M(1,0) = nin->axis[base+0].spaceDirection[1];
100        M(2,0) = nin->axis[base+0].spaceDirection[2];
101        M(0,1) = nin->axis[base+1].spaceDirection[0];
102        M(1,1) = nin->axis[base+1].spaceDirection[1];
103        M(2,1) = nin->axis[base+1].spaceDirection[2];
104        M(0,2) = nin->axis[base+2].spaceDirection[0];
105        M(1,2) = nin->axis[base+2].spaceDirection[1];
106        M(2,2) = nin->axis[base+2].spaceDirection[2];
107    
108      // Image location
109        *t = vec3f(
110            nin->spaceOrigin[0],
111            nin->spaceOrigin[1],
112            nin->spaceOrigin[2]);
113    
114    #undef M
115    }
116    
117    
118  /* Transformation matrix operations */  /* Transformation matrix operations */
119    
120  STATIC_INLINE Diderot_real_t Det2 (  STATIC_INLINE Diderot_real_t Det2 (
# Line 102  Line 169 
169   *  \param m the matrix to invert   *  \param m the matrix to invert
170   *  \param i the inverted matrix   *  \param i the inverted matrix
171   */   */
172  static void InvertM2x2 (Diderot_Mat2x2_t m, Diderot_Mat2x2_t i)  static void InvertM2x2 (Diderot_Mat2x2_t i, Diderot_Mat2x2_t m)
173  {  {
174  #define M(i)    m[(i)/2].r[(i)%2]  #define M(i)    m[(i)/2].r[(i)%2]
175  #define I(j)    i[(j)/2].r[(j)%2]  #define I(j)    i[(j)/2].r[(j)%2]
# Line 120  Line 187 
187   *  \param m the matrix to invert   *  \param m the matrix to invert
188   *  \param i the inverted matrix   *  \param i the inverted matrix
189   */   */
190  static void InvertTransposeM2x2 (Diderot_Mat2x2_t m, Diderot_Mat2x2_t i)  static void InvertM3x3 (Diderot_Mat3x3_t i, Diderot_Mat3x3_t m)
 {  
 #define M(i)    m[(i)/2].r[(i)%2]  
 #define I(j)    i[(j)/2].r[(j)%2]  
         Diderot_real_t scale = 1.0 / Det2(M(0), M(1), M(2), M(3));  
   
         I(0) =  scale * M(3);  
         I(2) = -scale * M(1);  
         I(1) = -scale * M(2);  
         I(3) =  scale * M(0);  
 #undef M  
 #undef I  
 }  
   
 static void LoadTransformM2x2 (Nrrd *nin, Diderot_Mat2x2_t m)  
 {  
 #define M(i,j)  m[i].r[j]  
   
   // compute the offset to the first space axis  
     int base = nin->dim - nin->spaceDim;  
   
   // Image axis Scaling and Rotation  
     M(0,0) = nin->axis[base+0].spaceDirection[0];  
   
   // Image location  
     M(0,1) = nin->spaceOrigin[0];  
   
   // bottom row is 0 1  
     M(1,0) = 0;  
     M(1,1) = 1;  
   
 #undef M  
 }  
   
 /*! \brief compute the inverse of \arg m, storing the result in \arg i.  
  *  \param m the matrix to invert  
  *  \param i the inverted matrix  
  */  
 static void InvertM3x3 (Diderot_Mat3x3_t m, Diderot_Mat3x3_t i)  
191  {  {
192          Diderot_real_t scale = 1.0 / DetM3x3(m);          Diderot_real_t scale = 1.0 / DetM3x3(m);
193    
# Line 177  Line 206 
206  #undef I  #undef I
207  }  }
208    
 /*! \brief compute the inverse transpose of \arg m, storing the result in \arg i.  
  *  \param m the matrix to invert  
  *  \param i the inverted matrix  
  */  
 static void InvertTransposeM3x3 (Diderot_Mat3x3_t m, Diderot_Mat3x3_t i)  
 {  
         Diderot_real_t scale = 1.0 / DetM3x3(m);  
   
 #define M(i)    m[(i)/3].r[(i)%3]  
 #define I(j)    i[(j)/3].r[(j)%3]  
         I(0) = scale * Det2 (M(4), M(5), M(7), M(8));  
         I(3) = scale * Det2 (M(2), M(1), M(8), M(7));  
         I(6) = scale * Det2 (M(1), M(2), M(4), M(5));  
         I(1) = scale * Det2 (M(5), M(3), M(8), M(6));  
         I(4) = scale * Det2 (M(0), M(2), M(6), M(8));  
         I(7) = scale * Det2 (M(2), M(0), M(5), M(3));  
         I(2) = scale * Det2 (M(3), M(4), M(6), M(7));  
         I(5) = scale * Det2 (M(1), M(0), M(7), M(6));  
         I(8) = scale * Det2 (M(0), M(1), M(3), M(4));  
 #undef M  
 #undef I  
 }  
   
 static void LoadTransformM3x3 (Nrrd *nin, Diderot_Mat3x3_t m)  
 {  
 #define M(i,j)  m[i].r[j]  
   
   // compute the offset to the first space axis  
     int base = nin->dim - nin->spaceDim;  
   
   // Image axis Scaling and Rotation  
     M(0,0) = nin->axis[base+0].spaceDirection[0];  
     M(1,0) = nin->axis[base+0].spaceDirection[1];  
     M(0,1) = nin->axis[base+1].spaceDirection[0];  
     M(1,1) = nin->axis[base+1].spaceDirection[1];  
   
   // Image location  
     M(0,2) = nin->spaceOrigin[0];  
     M(1,2) = nin->spaceOrigin[1];  
   
   // bottom row is 0 0 1  
     M(2,0) = 0;  
     M(2,1) = 0;  
     M(2,2) = 1;  
   
 #undef M  
 }  
   
209  /*! \brief compute the inverse of \arg m, storing the result in \arg i.  /*! \brief compute the inverse of \arg m, storing the result in \arg i.
210   *  \param m the matrix to invert   *  \param m the matrix to invert
211   *  \param i the inverted matrix   *  \param i the inverted matrix
212   */   */
213  static void InvertM4x4 (Diderot_Mat4x4_t m, Diderot_Mat4x4_t i)  static void InvertM4x4 (Diderot_Mat4x4_t i, Diderot_Mat4x4_t m)
214  {  {
215          Diderot_real_t scale = 1.0 / DetM4x4(m);          Diderot_real_t scale = 1.0 / DetM4x4(m);
216    
# Line 302  Line 283 
283  #undef I  #undef I
284  }  }
285    
 /*! \brief compute the inverse transpose of \arg m, storing the result in \arg i.  
  *  \param m the matrix to invert  
  *  \param i the inverted matrix  
  */  
 static void InvertTransposeM4x4 (Diderot_Mat4x4_t m, Diderot_Mat4x4_t i)  
 {  
         Diderot_real_t scale = 1.0 / DetM4x4(m);  
   
 #define M(i)    m[(i)>>2].r[(i)&3]  
 #define I(j)    i[(j)>>2].r[(j)&3]  
         I( 0) =  Det3(M( 5),M( 6),M( 7),  
                       M( 9),M(10),M(11),  
                       M(13),M(14),M(15)) * scale;  
   
         I( 4) = -Det3(M( 1),M( 2),M( 3),  
                       M( 9),M(10),M(11),  
                       M(13),M(14),M(15)) * scale;  
   
         I( 8) =  Det3(M( 1),M( 2),M( 3),  
                       M( 5),M( 6),M( 7),  
                       M(13),M(14),M(15)) * scale;  
   
         I(12) = -Det3(M( 1),M( 2),M( 3),  
                       M( 5),M( 6),M( 7),  
                       M( 9),M(10),M(11)) * scale;  
   
         I( 1) = -Det3(M( 4),M( 6),M( 7),  
                       M( 8),M(10),M(11),  
                       M(12),M(14),M(15)) * scale;  
   
         I( 5) =  Det3(M( 0),M( 2),M( 3),  
                       M( 8),M(10),M(11),  
                       M(12),M(14),M(15)) * scale;  
   
         I( 9) = -Det3(M( 0),M( 2),M( 3),  
                       M( 4),M( 6),M( 7),  
                       M(12),M(14),M(15)) * scale;  
   
         I( 7) =  Det3(M( 0),M( 2),M( 3),  
                       M( 4),M( 6),M( 7),  
                       M( 8),M(10),M(11)) * scale;  
   
         I( 2) =  Det3(M( 4),M( 5),M( 7),  
                       M( 8),M( 9),M(11),  
                       M(12),M(13),M(15)) * scale;  
   
         I( 6) = -Det3(M( 0),M( 1),M( 3),  
                       M( 8),M( 9),M(11),  
                       M(12),M(13),M(15)) * scale;  
   
         I(10) =  Det3(M( 0),M( 1),M( 3),  
                       M( 4),M( 5),M( 7),  
                       M(12),M(13),M(15)) * scale;  
   
         I(14) = -Det3(M( 0),M( 1),M( 3),  
                       M( 4),M( 5),M( 7),  
                       M( 8),M( 9),M(11)) * scale;  
   
         I( 3) = -Det3(M( 4),M( 5),M( 6),  
                       M( 8),M( 9),M(10),  
                       M(12),M(13),M(14)) * scale;  
   
         I( 7) =  Det3(M( 0),M( 1),M( 2),  
                       M( 8),M( 9),M(10),  
                       M(12),M(13),M(14)) * scale;  
   
         I(11) = -Det3(M( 0),M( 1),M( 2),  
                       M( 4),M( 5),M( 6),  
                       M(12),M(13),M(14)) * scale;  
   
         I(15) =  Det3(M( 0),M( 1),M( 2),  
                       M( 4),M( 5),M( 6),  
                       M( 8),M( 9),M(10)) * scale;  
   
 #undef M  
 #undef I  
 }  
   
 static void LoadTransformM4x4 (Nrrd *nin, Diderot_Mat4x4_t m)  
 {  
 #define M(i,j)  m[i].r[j]  
   
   // compute the offset to the first space axis  
     int base = nin->dim - nin->spaceDim;  
   
   // Image axis Scaling and Rotation  
     M(0,0) = nin->axis[base+0].spaceDirection[0];  
     M(1,0) = nin->axis[base+0].spaceDirection[1];  
     M(2,0) = nin->axis[base+0].spaceDirection[2];  
     M(0,1) = nin->axis[base+1].spaceDirection[0];  
     M(1,1) = nin->axis[base+1].spaceDirection[1];  
     M(2,1) = nin->axis[base+1].spaceDirection[2];  
     M(0,2) = nin->axis[base+2].spaceDirection[0];  
     M(1,2) = nin->axis[base+2].spaceDirection[1];  
     M(2,2) = nin->axis[base+2].spaceDirection[2];  
   
   // Image location  
     M(0,3) = nin->spaceOrigin[0];  
     M(1,3) = nin->spaceOrigin[1];  
     M(2,3) = nin->spaceOrigin[2];  
   
   // bottom row is 0 0 0 1  
     M(3,0) = 0;  
     M(3,1) = 0;  
     M(3,2) = 0;  
     M(3,3) = 1;  
   
 #undef M  
 }  
   
286  static Nrrd *LoadNrrdFile (const char *filename)  static Nrrd *LoadNrrdFile (const char *filename)
287  {  {
288    /* create a nrrd; at this point this is just an empty container */    /* create a nrrd; at this point this is just an empty container */
# Line 465  Line 336 
336    
337  printf("LoadImage \"%s\": space dim = 1, axes = <%d>\n",  printf("LoadImage \"%s\": space dim = 1, axes = <%d>\n",
338  name, img->size[base+0]);  name, img->size[base+0]);
339      LoadTransformM2x2 (nin, img->m);  
340      InvertM2x2 (img->m, img->mInv);    // from the Nrrd file, we load the scaling and translation
341  PrintMat2x2("m", img->m);      float s, t;
342  PrintMat2x2("mInv", img->mInv);      LoadTransform1D (nin, &s, &t);
343    
344        img->s = 1.0 / s;
345        img->t = -t / s;
346    
347      *imgOut = img;      *imgOut = img;
348      return DIDEROT_OK;      return DIDEROT_OK;
# Line 504  Line 378 
378    
379  printf("LoadImage \"%s\": space dim = 2, axes = <%d, %d>\n",  printf("LoadImage \"%s\": space dim = 2, axes = <%d, %d>\n",
380  name, img->size[0], img->size[1]);  name, img->size[0], img->size[1]);
381      LoadTransformM3x3 (nin, img->m);  
382      InvertM3x3 (img->m, img->mInv);    // from the Nrrd file, we load the affine image-to-world transform matrix
383      InvertTransposeM3x3 (img->m, img->mInvT);      Diderot_Mat2x2_t m;         // rotation and scaling
384  PrintMat3x3("m", img->m);      vec2f_t t;                  // translation part of the transform
385  PrintMat3x3("mInv", img->mInv);      LoadTransform2D (nin, m, &t);
386  PrintMat3x3("mInvT", img->mInvT);  
387    PrintMat2x2("m", m);
388    printf("t = <%f, %f>\n", ((union2f_t)t).r[0], ((union2f_t)t).r[1]);
389    
390      // compute inverse of m, which is the transform from the world basis to the image basis
391        InvertM2x2 (img->w2i, m);
392        img->tVec = -mulMat2x2Vec2f(img->w2i, t);  // inv([M t]) = [inv(M) -inv(M)t]
393    PrintMat2x2("w2i", img->w2i);
394    printf("tVec = <%f, %f>\n", ((union2f_t)img->tVec).r[0], ((union2f_t)img->tVec).r[1]);
395    
396      // compute transpose of w2i
397        transpose2x2f (img->w2iT, img->w2i);
398    
399      *imgOut = img;      *imgOut = img;
400      return DIDEROT_OK;      return DIDEROT_OK;
# Line 548  Line 433 
433    
434  printf("LoadImage \"%s\": space dim = 3, axes = <%d, %d, %d>\n",  printf("LoadImage \"%s\": space dim = 3, axes = <%d, %d, %d>\n",
435  name, img->size[0], img->size[1], img->size[2]);  name, img->size[0], img->size[1], img->size[2]);
436      LoadTransformM4x4 (nin, img->m);  
437      InvertM4x4 (img->m, img->mInv);    // from the Nrrd file, we load the affine image-to-world transform matrix
438      InvertTransposeM4x4 (img->m, img->mInvT);      Diderot_Mat3x3_t m;         // rotation and scaling
439  PrintMat4x4("m", img->m);      vec3f_t t;                  // translation part of the transform
440  PrintMat4x4("mInv", img->mInv);      LoadTransform3D (nin, m, &t);
441  PrintMat4x4("mInvT", img->mInvT);  
442    PrintMat3x3("m", m);
443    printf("t = <%f, %f, %f>\n", ((union3f_t)t).r[0], ((union3f_t)t).r[1], ((union3f_t)t).r[2]);
444    
445      // compute inverse of m, which is the transform from the world basis to the image basis
446        InvertM3x3 (img->w2i, m);
447        img->tVec = -mulMat3x3Vec3f(img->w2i, t);  // inv([M t]) = [inv(M) -inv(M)t]
448    PrintMat3x3("w2i", img->w2i);
449    printf("tVec = <%f, %f, %f>\n", ((union3f_t)img->tVec).r[0], ((union3f_t)img->tVec).r[1], ((union3f_t)img->tVec).r[2]);
450    
451      // compute transpose of w2i
452        transpose3x3f (img->w2iT, img->w2i);
453    
454      *imgOut = img;      *imgOut = img;
455      return DIDEROT_OK;      return DIDEROT_OK;

Legend:
Removed from v.990  
changed lines
  Added in v.991

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