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

SCM Repository

[diderot] Diff of /branches/vis15/src/lib/common/image.cxx
ViewVC logotype

Diff of /branches/vis15/src/lib/common/image.cxx

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

revision 3928, Sun Jun 5 11:24:36 2016 UTC revision 3929, Sun Jun 5 11:55:07 2016 UTC
# Line 1  Line 1 
1  /*! \file image.cxx  /*! \file image.cxx
2   *   *
3     * Utility code to support images in Diderot.  Most of the actual implementation can be
4     * found in diderot/image-inst.hxx; this file just contains the shared non-template code.
5     *
6   * \author John Reppy   * \author John Reppy
7   */   */
8    
# Line 15  Line 18 
18    
19  namespace diderot {  namespace diderot {
20    
21  /******************** 1D images ********************/      namespace __details {
22    
23  static void load_transform_1d (const Nrrd *nin, double &s, double &t)          void load_transform_1d (const Nrrd *nin, double &s, double &t)
24  {  {
25    // compute the offset to the first space axis    // compute the offset to the first space axis
26      int base = nin->dim - nin->spaceDim;      int base = nin->dim - nin->spaceDim;
# Line 29  Line 32 
32      t = nin->spaceOrigin[0];      t = nin->spaceOrigin[0];
33  }  }
34    
 template <typename REAL>  
 __details::image1d<REAL>::~image1d ()  
 {  
     free (this->_data);  
 }  
   
 template <typename REAL>  
 bool image1d<REAL>::load (struct world_base *wrld, std::string const &name)  
 {  
     Nrrd *nin = load_nrrd_file (wrld, name);  
     if (nin == nullptr)  
         return true;  
     else {  
         bool sts = this->load (wrld, nin);  
         nrrdNix (nin);  
         return sts;  
     }  
 }  
35    
36  template <typename REAL>          void load_transform_2d (const Nrrd *nin, mat2x2<double> &m, real2<double> &t)
 bool image1d<REAL>::load (struct world_base *wrld, const Nrrd *nin)  
 {  
     this->_img = std::make_shared<__details::image1d<REAL>>();  
   
   // compute the offset to the first space axis  
     int base = nin->dim - nin->spaceDim;  
   
 // FIXME: we should also be checking the shape of the voxels and the type of the samples.  
   // check that the image format is what we expect  
     if (nin->spaceDim != 1) {  
         error (wrld, "nrrd has unexpected dimension %d, expected 1\n", nin->spaceDim);  
         return true;  
     }  
     else if (nin->axis[base].kind != nrrdKindSpace) {  
         error (wrld, "nrrd has unexpected axis structure: %s\n",  
             airEnumStr(nrrdKind, nin->axis[base].kind));  
         return true;  
     }  
   
     this->_img->_dim     = 1;  
     this->_img->_size[0] = nin->axis[base+0].size;  
     this->_img->_data    = reinterpret_cast<REAL *>(nin->data);  
     this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);  
   
 //printf("SetImage \"%s\": space dim = 1, axes = <%d>\n",  
 //name, img->size[base+0]);  
   
   // from the Nrrd file, we load the scaling and translation  
     double s, t;  
   
     load_transform_1d (nin, s, t);  
   
     this->_img->_s = static_cast<REAL>(1.0 / s);  
     this->_img->_t = static_cast<REAL>(-t / s);  
   
     return false;  
 }  
   
 /******************** 2D images ********************/  
   
 static void load_transform_2d (const Nrrd *nin, mat2x2<double> &m, real2<double> &t)  
37  {  {
38  #define M(i,j)  m[(i)*2+(j)]  #define M(i,j)  m[(i)*2+(j)]
39    
# Line 114  Line 58 
58      return (a*d - b*c);      return (a*d - b*c);
59  }  }
60    
61  /*! \brief compute the inverse of \arg m, storing the result in \arg i.          void invert2x2 (mat2x2<double> &i, mat2x2<double> &m)
  *  \param m the matrix to invert  
  *  \param i the inverted matrix  
  */  
 inline void invert (mat2x2<double> &i, mat2x2<double> &m)  
62  {  {
63          double 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]);
64    
# Line 128  Line 68 
68          i[3] =  scale * m[0];          i[3] =  scale * m[0];
69  }  }
70    
71  template <typename REAL>          void load_transform_3d (const Nrrd *nin, mat3x3<double> &m, real3<double> &t)
 __details::image2d<REAL>::~image2d ()  
 {  
     free (this->_data);  
 }  
   
 template <typename REAL>  
 bool image2d<REAL>::load (struct world_base *wrld, std::string const &name)  
 {  
     Nrrd *nin = load_nrrd_file (wrld, name);  
     if (nin == nullptr)  
         return true;  
     else {  
         bool sts = this->load (wrld, nin);  
         nrrdNix (nin);  
         return sts;  
     }  
 }  
   
 template <typename REAL>  
 bool image2d<REAL>::load (struct world_base *wrld, const Nrrd *nin)  
 {  
     this->_img = std::make_shared<__details::image2d<REAL>>();  
   
   // compute the offset to the first space axis  
     int base = nin->dim - nin->spaceDim;  
   
   // check that the image format is what we expect  
     if (nin->spaceDim != 2) {  
         error (wrld, "nrrd has unexpected dimension %d, expected 2\n", nin->spaceDim);  
         return true;  
     }  
     else if ((nin->axis[base].kind != nrrdKindSpace)  
     ||       (nin->axis[base+1].kind != nrrdKindSpace)) {  
         error (wrld, "nrrd has unexpected axis structure: %s %s\n",  
             airEnumStr(nrrdKind, nin->axis[base].kind),  
             airEnumStr(nrrdKind, nin->axis[base+1].kind));  
         return true;  
     }  
   
     this->_img->_dim     = 2;  
     this->_img->_size[0] = nin->axis[base+0].size;  
     this->_img->_size[1] = nin->axis[base+1].size;  
     this->_img->_data    = reinterpret_cast<REAL *>(nin->data);  
     this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);  
   
   // from the Nrrd file, we load the affine image-to-world transform matrix  
     mat2x2<double> m;                   // rotation and scaling  
     real2<double> t;                    // translation part of the transform  
     load_transform_2d (nin, m, t);  
   
   // compute inverse of m, which is the transform from the world basis to the image basis  
     mat2x2<double> mInv;  
     invert (mInv, m);  
   
   // copy inverted matrix into the image data structure  
     for (int i = 0;  i < 2;  i++) {  
         for (int j = 0;  j < 2;  j++) {  
             REAL r = static_cast<REAL>(mInv[2*i+j]);  
             this->_img->_w2i[2*i+j] = r;  
             this->_img->_w2iT[2*j+i] = r;       // transpose  
         }  
     }  
   
   // transform the translation vector: inv([M t]) = [inv(M) -inv(M)t]  
     for (int i = 0;  i < 2;  i++) {  
         double sum = 0;  
         for (int j = 0;  j < 2;  j++) {  
             sum += mInv[2*i+j] * t[j];  
         }  
         this->_img->_tVec[i] = static_cast<REAL>(-sum);  
     }  
   
     return false;  
 }  
   
 /******************** 3D images ********************/  
   
 static void load_transform_3d (const Nrrd *nin, mat3x3<double> &m, real3<double> &t)  
72  {  {
73  #define M(i,j)  m[3*(i)+(j)]  #define M(i,j)  m[3*(i)+(j)]
74    
# Line 236  Line 98 
98   *  \param m the matrix to invert   *  \param m the matrix to invert
99   *  \param i the inverted matrix   *  \param i the inverted matrix
100   */   */
101  inline void invert (mat3x3<double> &i, mat3x3<double> &m)          void invert3x3 (mat3x3<double> &i, mat3x3<double> &m)
102  {  {
103          double det =          double det =
104                  ( m[0]*m[4]*m[8]                  ( m[0]*m[4]*m[8]
# Line 258  Line 120 
120          i[8] = scale * det2 (m[0], m[1], m[3], m[4]);          i[8] = scale * det2 (m[0], m[1], m[3], m[4]);
121  }  }
122    
123  template <typename REAL>      } // namespace __details
 __details::image3d<REAL>::~image3d ()  
 {  
     free (this->_data);  
 }  
   
 template <typename REAL>  
 bool image3d<REAL>::load (struct world_base *wrld, std::string const &name)  
 {  
     Nrrd *nin = load_nrrd_file (wrld, name);  
     if (nin == nullptr)  
         return true;  
     else {  
         bool sts = this->load (wrld, nin);  
         nrrdNix (nin);  
         return sts;  
     }  
 }  
   
 template <typename REAL>  
 bool image3d<REAL>::load (struct world_base *wrld, const Nrrd *nin)  
 {  
     this->_img = std::make_shared<__details::image3d<REAL>>();  
   
   // compute the offset to the first space axis  
     int base = nin->dim - nin->spaceDim;  
   
   // check that the image format is what we expect  
     if (nin->spaceDim != 3) {  
         error (wrld, "nrrd has unexpected dimension %d, expected 3\n", nin->spaceDim);  
         return true;  
     }  
     else if ((nin->axis[base].kind != nrrdKindSpace)  
     ||       (nin->axis[base+1].kind != nrrdKindSpace)  
     ||       (nin->axis[base+2].kind != nrrdKindSpace)) {  
         error (wrld, "nrrd has unexpected axis structure: %s %s %s\n",  
             airEnumStr(nrrdKind, nin->axis[base].kind),  
             airEnumStr(nrrdKind, nin->axis[base+1].kind),  
             airEnumStr(nrrdKind, nin->axis[base+2].kind));  
         return true;  
     }  
   
     this->_img->_dim     = 3;  
     this->_img->_size[0] = nin->axis[base+0].size;  
     this->_img->_size[1] = nin->axis[base+1].size;  
     this->_img->_size[2] = nin->axis[base+2].size;  
     this->_img->_data    = reinterpret_cast<REAL *>(nin->data);  
     this->_img->_dataSzb = nrrdElementSize(nin) * nrrdElementNumber(nin);  
   
   // from the Nrrd file, we load the affine image-to-world transform matrix  
     mat3x3<double> m;                   // rotation and scaling  
     real3<double> t;                    // translation part of the transform  
     load_transform_3d (nin, m, t);  
   
   // compute inverse of m, which is the transform from the world basis to the image basis  
     mat3x3<double> mInv;  
     invert (mInv, m);  
   
   // copy results into the image data structure  
     for (int i = 0;  i < 3;  i++) {  
         for (int j = 0;  j < 3;  j++) {  
             REAL r = static_cast<REAL>(mInv[2*i+j]);  
             this->_img->_w2i[3*i+j]  = r;  
             this->_img->_w2iT[3*j+i] = r;       // transpose  
         }  
     }  
   
   // transform the translation vector: inv([M t]) = [inv(M) -inv(M)t]  
     for (int i = 0;  i < 3;  i++) {  
         double sum = 0;  
         for (int j = 0;  j < 3;  j++) {  
             sum += mInv[3*i+j] * t[j];  
         }  
         this->_img->_tVec[i] = static_cast<REAL>(-sum);  
     }  
   
     return false;  
 }  
124    
125  } // namespace Diderot  } // namespace Diderot
   
 void foo (std::string name)  
 {  
     diderot::image1d<float> img1;  
     diderot::image2d<float> img2;  
     diderot::image3d<float> img3;  
   
     img1.load(nullptr, name);  
     img2.load(nullptr, name);  
     img3.load(nullptr, name);  
   
     float x = img2.translate()[0];  
   
 }  

Legend:
Removed from v.3928  
changed lines
  Added in v.3929

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