Fraxinus  18.10
An IGT application
metaimage.hpp
Go to the documentation of this file.
1 #ifndef METAIMAGE_HPP
2 #define METAIMAGE_HPP
3 
4 #include <vtkSmartPointer.h>
5 #include <vtkMetaImageReader.h>
6 #include <vtkImageData.h>
7 #include "ErrorHandler.hpp"
8 
14 template<typename T>
15 class MetaImage {
16 public:
21  {
22  m_xsize = 0;
23  m_ysize = 0;
24  m_img = vtkSmartPointer<vtkImageData>::New();
25  m_xspacing = 0.0;
26  m_yspacing = 0.0;
27  m_transform = Matrix4::Zero();
28  m_idx = -1;
29  }
30 
32  {
33 
34  }
35 
39  int
40  getXSize() const
41  {
42  return m_xsize;
43  }
47  int
48  getYSize() const
49  {
50  return m_ysize;
51  }
52 
56  double
57  getXSpacing() const
58  {
59  return m_xspacing;
60  }
64  double
65  getYSpacing() const
66  {
67  return m_yspacing;
68  }
73  int
74  makehash(const int x, const int y) const
75  {
76  return x + getXSize()*y;
77  }
78 
79 
83  T*
85  {
86  return (T*)m_img->GetScalarPointer();
87  }
88 
92  const T*
94  {
95  return (T*)m_img->GetScalarPointer();
96  }
97 
104  void
105  toImgCoords(double &x,
106  double &y,
107  const double p[3]) const
108  {
109 
110 
111  // Thing is, the matrices we are given are the opposite transform,
112  // so we have to do the inverse transformation: p = R^Tp - R^Tp_0
113  double offset_x, offset_y;
114  x =
115  m_transform(0,0)*p[0] +
116  m_transform(1,0)*p[1] +
117  m_transform(2,0)*p[2];
118  // Now subtract the transformed offset
119  offset_x =
120  m_transform(0,0)*m_transform(0,3) +
121  m_transform(1,0)*m_transform(1,3) +
122  m_transform(2,0)*m_transform(2,3);
123  x = (x - offset_x)/getXSpacing();
124  y =
125  m_transform(0,1)*p[0] +
126  m_transform(1,1)*p[1] +
127  m_transform(2,1)*p[2];
128  // Now subtract the transformed offset
129  offset_y =
130  m_transform(0,1)*m_transform(0,3) +
131  m_transform(1,1)*m_transform(1,3) +
132  m_transform(2,1)*m_transform(2,3);
133 
134  y = (y - offset_y)/getYSpacing();
135 
136 
137  }
138 
139 
140 
146  Matrix4
148  {
149  return m_transform;
150  }
151 
157  const Matrix4&
158  getTransform() const
159  {
160  return m_transform;
161  }
162 
170  template<typename Dt>
171  void regionGrow(vector<Dt>& ret, int imgx, int imgy) const
172  {
173 
174  vector<pair<int,int> > ptstack;
175 
176  // Get image data
177 
178  const T *imagedata = getPixelPointer();
179 
180  pair<int,int> cur;
181  T curpt;
182 
183  // Initialize stack with seed point
184  ptstack.push_back(make_pair(imgx, imgy));
185  curpt = imagedata[imgx + imgy*getXSize()];
186 
187  ret.push_back(curpt);
188 
189  // Use a bool array instead of unordered_map, it's just faster even though we might need a little more memory
190  vector<bool> visited;
191  visited.resize(makehash( getXSize(), getYSize()));
192  std::fill(visited.begin(), visited.end(), false);
193 
194  // Go!
195  while(!ptstack.empty())
196  {
197  cur = ptstack.back();
198 
199  ptstack.pop_back();
200  const int img_x = cur.first;
201  const int img_y = cur.second;
202 
203  if(!inImage(img_x, img_y))
204  {
205  continue;
206  }
207  curpt = imagedata[img_x + img_y*getXSize()];
208 
209  // Is this point in?
210  if(curpt != 0)
211  {
212  if(visited.at(makehash(img_x, img_y)) )
213  continue;
214 
215  ret.push_back(curpt);
216  // Check the neighbours (4-connected)
217  ptstack.insert(ptstack.end(),
218  { make_pair(img_x-1,img_y),
219  make_pair(img_x+1,img_y),
220  make_pair(img_x,img_y-1),
221  make_pair(img_x,img_y+1)});
222  }
223  visited.at(makehash(img_x, img_y)) = true;
224  }
225 
226  }
227 
232  inline
233  void setIdx(int i)
234  {
235  m_idx = i;
236  }
237 
242  inline
243  int getIdx() const
244  {
245  return m_idx;
246  }
247 
248 
254  static vector<MetaImage>* readImages(const string & prefix)
255  {
256  // Images are on the format prefix$NUMBER.mhd
257  vtkSmartPointer<vtkMetaImageReader> reader= vtkSmartPointer<vtkMetaImageReader>::New();
258  vtkSmartPointer<ErrorObserver> errorObserver = vtkSmartPointer<ErrorObserver>::New();
259  reader->AddObserver(vtkCommand::ErrorEvent,errorObserver);
260 
261  int i = 0;
262  ostringstream ss;
263  ss << prefix << i << ".mhd";
264  string filename = ss.str();
265  reader->SetFileName(filename.c_str());
266  reader->Update();
267 
268  if(!reader->CanReadFile(filename.c_str())){
269  cerr << filename.c_str() << endl;
270  reportError("ERROR: Could not read velocity data \n");
271  }
272 
273  vector<MetaImage> *ret = new vector<MetaImage>();
274  while(reader->CanReadFile(filename.c_str()))
275  {
276  ret->push_back(MetaImage());
277  ret->at(i).setIdx(i);
278  ret->at(i).m_img->DeepCopy(reader->GetOutput());
279 
280 #if VTK_MAJOR_VERSION <= 5
281  ret->at(i).m_img->Update();
282 #else
283 #endif
284 
285  if (errorObserver->GetError())
286  {
287  reportError("ERROR: Could not read velocity data \n"+ errorObserver->GetErrorMessage());
288  }
289  if (errorObserver->GetWarning()){
290  cerr << "Caught warning while reading velocity data! \n " << errorObserver->GetWarningMessage();
291  }
292 
293  if ( reader->GetFileDimensionality() != 2){
294  reportError("ERROR: Can only read 2-D data");
295  }
296 
297 
298  ret->at(i).m_xsize = reader->GetWidth();
299  ret->at(i).m_ysize = reader->GetHeight();
300  double *spacing;
301  spacing = reader->GetPixelSpacing();
302  ret->at(i).m_xspacing = spacing[0];
303  ret->at(i).m_yspacing = spacing[1];
304 
305  std::ifstream infile(filename);
306  std::string line;
307  int found =0;
308  int numToFind =2;
309  while (std::getline(infile, line))
310  {
311  if(line.find("Offset")!=std::string::npos)
312  {
313  string buf;
314  stringstream ss(line);
315  int k=0;
316  ss >> buf;
317  ss >> buf;
318  while (ss >> buf)
319  {
320  ret->at(i).m_transform(k++,3)=std::stod(buf);
321  }
322  ret->at(i).m_transform(3,3)=1;
323  found++;
324  if(found >=numToFind) break;
325 
326  }else if(line.find("TransformMatrix")!=std::string::npos)
327  {
328  string buf;
329  stringstream ss(line);
330  int k=0;
331  int j=0;
332  ss >> buf;
333  ss >> buf;
334  while (ss >> buf)
335  {
336  if(k >2){
337  k=0;
338  j++;
339  }
340  ret->at(i).m_transform(k++,j)=std::stod(buf);
341  }
342  found++;
343  if(found >=numToFind) break;
344  }
345  }
346  ss.clear();
347  ss.str("");
348  ss << prefix << ++i << ".mhd";
349  filename = ss.str();
350  reader->SetFileName(filename.c_str());
351  reader->Update();
352  }
353  return ret;
354  }
355 
362  inline bool
363  inImage(double img_x,double img_y) const
364  {
365  return !(img_x < 0
366  || img_x >= getXSize()
367  || img_y < 0
368  || img_y >= getYSize());
369  }
370 
371 private:
372  vtkSmartPointer<vtkImageData> m_img;
373  int m_idx;
374  int m_xsize;
375  int m_ysize;
376  double m_xspacing;
377  double m_yspacing;
378  Matrix4 m_transform;
379 };
380 
381 #endif //METAIMAGE_HPP
void setIdx(int i)
Definition: metaimage.hpp:233
double getXSpacing() const
Definition: metaimage.hpp:57
const Matrix4 & getTransform() const
Definition: metaimage.hpp:158
int makehash(const int x, const int y) const
Definition: metaimage.hpp:74
int getIdx() const
Definition: metaimage.hpp:243
void regionGrow(vector< Dt > &ret, int imgx, int imgy) const
Definition: metaimage.hpp:171
Eigen::Matrix< double, 4, 4 > Matrix4
Definition: matrix.hpp:5
static vector< MetaImage > * readImages(const string &prefix)
Definition: metaimage.hpp:254
void toImgCoords(double &x, double &y, const double p[3]) const
Definition: metaimage.hpp:105
void fill(Eigen::Affine3d *self, vtkMatrix4x4Ptr m)
void reportError(std::string errMsg)
double getYSpacing() const
Definition: metaimage.hpp:65
int getYSize() const
Definition: metaimage.hpp:48
T * getPixelPointer()
Definition: metaimage.hpp:84
bool inImage(double img_x, double img_y) const
Definition: metaimage.hpp:363
Matrix4 getTransform()
Definition: metaimage.hpp:147
int getXSize() const
Definition: metaimage.hpp:40
const T * getPixelPointer() const
Definition: metaimage.hpp:93