Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
itkContourExtractor2DImageFilter.txx
Go to the documentation of this file.
1 /*===================================================================
2 
3 The Medical Imaging Interaction Toolkit (MITK)
4 
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
7 All rights reserved.
8 
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
11 A PARTICULAR PURPOSE.
12 
13 See LICENSE.txt or http://www.mitk.org for details.
14 
15 ===================================================================*/
16 
17 /*===================================================================
18 
19 This file is based heavily on a corresponding ITK filter.
20 
21 ===================================================================*/
22 
23 #ifndef __itkContourExtractor2DImageFilter_txx
24 #define __itkContourExtractor2DImageFilter_txx
25 
26 #include "itkConstShapedNeighborhoodIterator.h"
27 #include "itkConstShapedNeighborhoodIterator.h"
28 #include "itkContourExtractor2DImageFilter.h"
29 #include "itkProgressReporter.h"
30 #include "vcl_cmath.h"
31 
32 namespace itk
33 {
34  // Constructor
35  template <class TInputImage>
36  ContourExtractor2DImageFilter<TInputImage>::ContourExtractor2DImageFilter()
37  {
38  this->m_ContourValue = NumericTraits<InputRealType>::Zero;
39  this->m_ReverseContourOrientation = false;
40  this->m_VertexConnectHighPixels = false;
41  this->m_UseCustomRegion = false;
42  this->m_NumberOfContoursCreated = 0;
43  }
44 
45  // Destructor
46  template <class TInputImage>
47  ContourExtractor2DImageFilter<TInputImage>::~ContourExtractor2DImageFilter()
48  {
49  }
50 
51  template <class TInputImage>
52  void ContourExtractor2DImageFilter<TInputImage>::GenerateData()
53  {
54  // Make sure the structures for containing, looking up, and numbering the
55  // growing contours are empty and ready.
56  m_Contours.clear();
57  m_ContourStarts.clear();
58  m_ContourEnds.clear();
59  m_NumberOfContoursCreated = 0;
60 
61  // Set up an iterator to "march the squares" across the image.
62  // We associate each 2px-by-2px square with the pixel in the upper left of
63  // that square. We then iterate across the image, examining these 2x2 squares
64  // and building the contour. By iterating the upper-left pixel of our
65  // "current square" across every pixel in the image except those on the
66  // bottom row and rightmost column, we have visited every valid square in the
67  // image.
68 
69  InputRegionType region = this->GetInput()->GetRequestedRegion();
70  typename InputRegionType::SizeType shrunkSize = region.GetSize();
71  shrunkSize[0] -= 1;
72  shrunkSize[1] -= 1;
73  InputRegionType shrunkRegion(region.GetIndex(), shrunkSize);
74 
75  // Set up a progress reporter
76  ProgressReporter progress(this, 0, shrunkRegion.GetNumberOfPixels());
77 
78  // A 1-pixel radius sets up a neighborhood with the following indices:
79  // 0 1 2
80  // 3 4 5
81  // 6 7 8
82  // We are interested only in the square of 4,5,7,8 which is the 2x2 square
83  // with the center pixel at the top-left. So we only activate the
84  // coresponding offsets, and only query pixels 4, 5, 7, and 8 with the
85  // iterator's GetPixel method.
86  typedef ConstShapedNeighborhoodIterator<InputImageType> SquareIterator;
87  typename SquareIterator::RadiusType radius = {{1, 1}};
88  SquareIterator it(radius, this->GetInput(), shrunkRegion);
89  InputOffsetType none = {{0, 0}};
90  InputOffsetType right = {{1, 0}};
91  InputOffsetType down = {{0, 1}};
92  InputOffsetType diag = {{1, 1}};
93  it.ActivateOffset(none);
94  it.ActivateOffset(right);
95  it.ActivateOffset(down);
96  it.ActivateOffset(diag);
97 
98  for (it.GoToBegin(); !it.IsAtEnd(); ++it)
99  {
100  // There are sixteen different possible square types, diagramed below.
101  // A + indicates that the vertex is above the contour value, and a -
102  // indicates that the vertex is below or equal to the contour value.
103  // The vertices of each square are here numbered:
104  // 01
105  // 23
106  // and treated as a binary value with the bits in that order. Thus each
107  // square can be so numbered:
108  // 0-- 1+- 2-+ 3++ 4-- 5+- 6-+ 7++
109  // -- -- -- -- +- +- +- +-
110  //
111  // 8-- 9+- 10-+ 11++ 12-- 13+- 14-+ 15++
112  // -+ -+ -+ -+ ++ ++ ++ ++
113  //
114  // The position of the line segment that cuts through (or doesn't, in case
115  // 0 and 15) each square is clear, except in cases 6 and 9. In this case,
116  // where the segments are placed is determined by
117  // m_VertexConnectHighPixels. If m_VertexConnectHighPixels is false, then
118  // lines like are drawn through square 6, and lines like are drawn through
119  // square 9. Otherwise, the situation is reversed.
120  // Finally, recall that we draw the lines so that (moving from tail to
121  // head) the lower-valued pixels are on the left of the line. So, for
122  // example, case 1 entails a line slanting from the middle of the top of
123  // the square to the middle of the left side of the square.
124 
125  // (1) Determine what number square we are currently inspecting. Remember
126  // that as far as the neighborhood iterator is concerned, our square
127  // 01 is numbered as 45
128  // 23 78
129 
130  InputPixelType v0, v1, v2, v3;
131  v0 = it.GetPixel(4);
132  v1 = it.GetPixel(5);
133  v2 = it.GetPixel(7);
134  v3 = it.GetPixel(8);
135  InputIndexType index = it.GetIndex();
136  unsigned char squareCase = 0;
137  if (v0 > m_ContourValue)
138  squareCase += 1;
139  if (v1 > m_ContourValue)
140  squareCase += 2;
141  if (v2 > m_ContourValue)
142  squareCase += 4;
143  if (v3 > m_ContourValue)
144  squareCase += 8;
145 
146 // Set up macros to find the ContinuousIndex where the contour intersects
147 // one of the sides of the square. Normally macros should, of course, be
148 // eschewed, but since this is an inner loop not calling the function four
149 // times when two would do is probably worth while. Plus, copy-pasting
150 // these into the switch below is even worse. InterpolateContourPosition
151 // takes the values at two vertices, the index of the first vertex, and the
152 // offset between the two vertices.
153 #define TOP_ this->InterpolateContourPosition(v0, v1, index, right)
154 #define BOTTOM_ this->InterpolateContourPosition(v2, v3, index + down, right)
155 #define LEFT_ this->InterpolateContourPosition(v0, v2, index, down)
156 #define RIGHT_ this->InterpolateContourPosition(v1, v3, index + right, down)
157 
158  // (2) Add line segments to the growing contours as defined by the cases.
159  // AddSegment takes a "from" vertex and a "to" vertex, and adds it to the
160  // a growing contour, creates a new contour, or merges two together.
161  switch (squareCase)
162  {
163  case 0: // no line
164  break;
165  case 1: // top to left
166  this->AddSegment(TOP_, LEFT_);
167  break;
168  case 2: // right to top
169  this->AddSegment(RIGHT_, TOP_);
170  break;
171  case 3: // right to left
172  this->AddSegment(RIGHT_, LEFT_);
173  break;
174  case 4: // left to bottom
175  this->AddSegment(LEFT_, BOTTOM_);
176  break;
177  case 5: // top to bottom
178  this->AddSegment(TOP_, BOTTOM_);
179  break;
180  case 6:
181  if (m_VertexConnectHighPixels)
182  {
183  // left to top
184  this->AddSegment(LEFT_, TOP_);
185  // right to bottom
186  this->AddSegment(RIGHT_, BOTTOM_);
187  }
188  else
189  {
190  // right to top
191  this->AddSegment(RIGHT_, TOP_);
192  // left to bottom
193  this->AddSegment(LEFT_, BOTTOM_);
194  }
195  break;
196  case 7: // right to bottom
197  this->AddSegment(RIGHT_, BOTTOM_);
198  break;
199  case 8: // bottom to right
200  this->AddSegment(BOTTOM_, RIGHT_);
201  break;
202  case 9:
203  if (m_VertexConnectHighPixels)
204  {
205  // top to right
206  this->AddSegment(TOP_, RIGHT_);
207  // bottom to left
208  this->AddSegment(BOTTOM_, LEFT_);
209  }
210  else
211  {
212  // top to left
213  this->AddSegment(TOP_, LEFT_);
214  // bottom to right
215  this->AddSegment(BOTTOM_, RIGHT_);
216  }
217  break;
218  case 10: // bottom to top
219  this->AddSegment(BOTTOM_, TOP_);
220  break;
221  case 11: // bottom to left
222  this->AddSegment(BOTTOM_, LEFT_);
223  break;
224  case 12: // left to right
225  this->AddSegment(LEFT_, RIGHT_);
226  break;
227  case 13: // top to right
228  this->AddSegment(TOP_, RIGHT_);
229  break;
230  case 14: // left to top
231  this->AddSegment(LEFT_, TOP_);
232  break;
233  case 15: // no line
234  break;
235  } // switch squareCase
236  progress.CompletedPixel();
237  } // iteration
238 
239  // Now create the outputs paths from the deques we've been using.
240  this->FillOutputs();
241  m_Contours.clear();
242  m_ContourStarts.clear();
243  m_ContourEnds.clear();
244  m_NumberOfContoursCreated = 0;
245  }
246 
247  template <class TInputImage>
248  inline typename ContourExtractor2DImageFilter<TInputImage>::VertexType
249  ContourExtractor2DImageFilter<TInputImage>::InterpolateContourPosition(InputPixelType fromValue,
250  InputPixelType toValue,
251  InputIndexType fromIndex,
252  InputOffsetType toOffset)
253  {
254  VertexType output;
255  // Now calculate the fraction of the way from 'from' to 'to' that the contour
256  // crosses. Interpolate linearly: y = v0 + (v1 - v0) * x, and solve for the
257  // x that gives y = m_ContourValue: x = (m_ContourValue - v0) / (v1 - v0).
258  // This assumes that v0 and v1 are separated by exactly ONE unit. So the to
259  // Offset. value must have exactly one component 1 and the other component 0.
260  // Also this assumes that fromValue and toValue are different. Otherwise we
261  // can't interpolate anything!
262  itkAssertOrThrowMacro((fromValue != toValue), "source and destination are the same");
263 
264  itkAssertOrThrowMacro(((toOffset[0] == 0 && toOffset[1] == 1) || (toOffset[0] == 1 && toOffset[1] == 0)),
265  "toOffset has unexpected values");
266 
267  double x =
268  (m_ContourValue - static_cast<InputRealType>(fromValue)) / (toValue - static_cast<InputRealType>(fromValue));
269 
270  output[0] = fromIndex[0] + x * toOffset[0];
271  output[1] = fromIndex[1] + x * toOffset[1];
272 
273  return output;
274  }
275 
276  template <class TInputImage>
277  void ContourExtractor2DImageFilter<TInputImage>::AddSegment(VertexType from, VertexType to)
278  {
279  if (from == to)
280  {
281  // Arc is degenerate: ignore, and the from/two point will be connected
282  // later by other squares. Degeneracy happens when (and only when) a square
283  // has exactly one vertex that is the contour value, and the rest are above
284  // that value.
285  return;
286  }
287 
288  // Try to find an existing contour that starts where the new segment ends.
289  VertexMapIterator newTail = m_ContourStarts.find(to);
290  // Try to find an existing contour that ends where the new segment starts.
291  VertexMapIterator newHead = m_ContourEnds.find(from);
292 
293  if (newTail != m_ContourStarts.end() && newHead != m_ContourEnds.end())
294  {
295  // We need to connect these two contours. The act of connecting them will
296  // add the needed arc.
297  auto tail = newTail->second;
298  itkAssertOrThrowMacro((tail->front() == to), "End doesn't match Beginning");
299  auto head = newHead->second;
300  itkAssertOrThrowMacro((head->back() == from), "Beginning doesn't match End");
301  if (head == tail)
302  {
303  // We've closed a contour. Add the end point, and remove from the maps
304  head->push_back(to);
305  m_ContourStarts.erase(newTail);
306  // erase the front of tail. Because head and tail are the same contour,
307  // don't worry about erasing the front of head!
308  m_ContourEnds.erase(newHead); // erase the end of head/tail.
309  }
310  else
311  {
312  // We have found two distinct contours that need to be joined. Careful
313  // here: we want to keep the first segment in the list when merging so
314  // that contours are always returned in top-to-bottom, right-to-left
315  // order (with regard to the image pixel found to be inside the contour).
316  if (tail->m_ContourNumber > head->m_ContourNumber)
317  {
318  // if tail was created later than head...
319  // Copy tail to the end of head and remove
320  // tail from everything.
321  head->insert(head->end(), tail->begin(), tail->end());
322 
323  // Now remove 'tail' from the list and the maps because it has been
324  // subsumed.
325  m_ContourStarts.erase(newTail);
326  int erased = m_ContourEnds.erase(tail->back());
327  // There should be exactly one entry in the hash for that endpoint
328  if (erased != 1)
329  {
330  itkWarningMacro(<< "There should be exactly one entry in the hash for that endpoint, but there are "
331  << erased);
332  }
333  m_Contours.erase(tail); // remove from the master list
334 
335  // Now remove the old end of 'head' from the ends map and add
336  // the new end.
337  m_ContourEnds.erase(newHead);
338  m_ContourEnds.insert(VertexContourRefPair(head->back(), head));
339  }
340  else
341  {
342  // Copy head to the beginning of tail and remove
343  // head from everything.
344  tail->insert(tail->begin(), head->begin(), head->end());
345 
346  // Now remove 'head' from the list and the maps because
347  // it has been subsumed.
348  m_ContourEnds.erase(newHead);
349  int erased = m_ContourStarts.erase(head->front());
350  if (erased != 1)
351  {
352  itkWarningMacro(<< "There should be exactly one entry in the hash for that endpoint, but there are "
353  << erased);
354  }
355  m_Contours.erase(head); // remove from the master list
356 
357  // Now remove the old start of 'tail' from the starts map and
358  // add the new start.
359  m_ContourStarts.erase(newTail);
360  m_ContourStarts.insert(VertexContourRefPair(tail->front(), tail));
361  }
362  }
363  }
364  else if (newTail == m_ContourStarts.end() && newHead == m_ContourEnds.end())
365  {
366  // No contours found: add a new one.
367  // Make it on the heap. It will be copied into m_Contours.
368  ContourType contour;
369 
370  // Add the endpoints
371  contour.push_front(from);
372  contour.push_back(to);
373  contour.m_ContourNumber = m_NumberOfContoursCreated++;
374  // Add the contour to the end of the list and get a reference to it.
375  m_Contours.push_back(contour);
376 
377  // recall that end() is an iterator to one past the back!
378  auto newContour = --m_Contours.end();
379  // add the endpoints and an iterator pointing to the contour
380  // in the list to the maps.
381  m_ContourStarts.insert(VertexContourRefPair(from, newContour));
382  m_ContourEnds.insert(VertexContourRefPair(to, newContour));
383  }
384  else if (newTail != m_ContourStarts.end() && newHead == m_ContourEnds.end())
385  {
386  // Found a single contour to which the new arc should be prepended.
387  auto tail = newTail->second;
388  itkAssertOrThrowMacro((tail->front() == to), "End doesn't match Beginning");
389  tail->push_front(from);
390  // erase the old start of this contour
391  m_ContourStarts.erase(newTail);
392  // Now add the new start of this contour.
393  m_ContourStarts.insert(VertexContourRefPair(from, tail));
394  }
395  else if (newTail == m_ContourStarts.end() && newHead != m_ContourEnds.end())
396  {
397  // Found a single contour to which the new arc should be appended.
398  auto head = newHead->second;
399  itkAssertOrThrowMacro((head->back() == from), "Beginning doesn't match End");
400  head->push_back(to);
401  // erase the old end of this contour
402  m_ContourEnds.erase(newHead);
403  // Now add the new start of this contour.
404  m_ContourEnds.insert(VertexContourRefPair(to, head));
405  }
406  }
407 
408  template <class TInputImage>
409  void ContourExtractor2DImageFilter<TInputImage>::FillOutputs()
410  {
411  this->SetNumberOfIndexedOutputs(m_Contours.size());
412  int i = 0;
413  for (auto it = m_Contours.begin(); it != m_Contours.end(); it++, i++)
414  {
415  OutputPathPointer output = this->GetOutput(i);
416  if (output.IsNull())
417  {
418  // Static cast is OK because we know PathSource will make its templated
419  // class type
420  output = static_cast<OutputPathType *>(this->MakeOutput(i).GetPointer());
421  this->SetNthOutput(i, output.GetPointer());
422  }
423  typename VertexListType::Pointer path = const_cast<VertexListType *>(output->GetVertexList());
424  path->Initialize();
425  path->reserve(it->size()); // use std::vector version of 'reserve()'
426  // instead of VectorContainer::Reserve() to work around
427  // the fact that the latter is essentially std::vector::resize(),
428  // which is not what we want.
429 
430  // Now put all the points from the contour deque into the path and
431  // mark output as modified
432 
433  typedef typename ContourType::const_iterator ConstIteratorType;
434  if (m_ReverseContourOrientation)
435  {
436  ConstIteratorType itC = (*it).end();
437  do
438  {
439  itC--;
440  path->push_back(*itC);
441  } while (itC != (*it).begin());
442  }
443  else
444  {
445  ConstIteratorType itC = (*it).begin();
446  while (itC != (*it).end())
447  {
448  path->push_back(*itC);
449  itC++;
450  }
451  }
452  output->Modified();
453  }
454  }
455 
456  template <class TInputImage>
457  void ContourExtractor2DImageFilter<TInputImage>::SetRequestedRegion(const InputRegionType region)
458  {
459  itkDebugMacro("setting RequestedRegion to " << region);
460  m_UseCustomRegion = true;
461  if (this->m_RequestedRegion != region)
462  {
463  this->m_RequestedRegion = region;
464  this->Modified();
465  }
466  }
467 
468  template <class TInputImage>
469  void ContourExtractor2DImageFilter<TInputImage>::ClearRequestedRegion()
470  {
471  itkDebugMacro("Clearing RequestedRegion.");
472  if (this->m_UseCustomRegion == true)
473  {
474  this->m_UseCustomRegion = false;
475  this->Modified();
476  }
477  }
478 
479  template <class TInputImage>
480  void ContourExtractor2DImageFilter<TInputImage>::GenerateInputRequestedRegion() throw(InvalidRequestedRegionError)
481  {
482  InputImageType *input = const_cast<InputImageType *>(this->GetInput());
483  if (!input)
484  return;
485 
486  if (m_UseCustomRegion)
487  {
488  InputRegionType requestedRegion = m_RequestedRegion;
489  if (requestedRegion.Crop(input->GetLargestPossibleRegion()))
490  {
491  input->SetRequestedRegion(requestedRegion);
492  return;
493  }
494  else
495  {
496  // Couldn't crop the region (requested region is outside the largest
497  // possible region). Throw an exception.
498 
499  // store what we tried to request (prior to trying to crop)
500  input->SetRequestedRegion(requestedRegion);
501 
502  // build an exception
503  InvalidRequestedRegionError e(__FILE__, __LINE__);
504  e.SetLocation(ITK_LOCATION);
505  e.SetDescription("Requested region is outside the largest possible region.");
506  e.SetDataObject(input);
507  throw e;
508  }
509  }
510  else
511  {
512  input->SetRequestedRegion(input->GetLargestPossibleRegion());
513  }
514  }
515 
516  /**
517  * Standard "PrintSelf" method
518  */
519  template <class TInputImage>
520  void ContourExtractor2DImageFilter<TInputImage>::PrintSelf(std::ostream &os, Indent indent) const
521  {
522  Superclass::PrintSelf(os, indent);
523  os << indent << "ReverseContourOrientation: " << m_ReverseContourOrientation << std::endl;
524  os << indent << "VertexConnectHighPixels: " << m_VertexConnectHighPixels << std::endl;
525  os << indent << "UseCustomRegion: " << m_UseCustomRegion << std::endl;
526  os << indent << "NumericTraits: " << m_UseCustomRegion << std::endl;
527  os << indent << "NumberOfContoursCreated: " << m_NumberOfContoursCreated << std::endl;
528  if (m_UseCustomRegion)
529  {
530  os << indent << "Custom region: " << m_RequestedRegion << std::endl;
531  }
532 
533  typedef typename NumericTraits<InputRealType>::PrintType InputRealPrintType;
534 
535  os << indent << "Contour value: " << static_cast<InputRealPrintType>(m_ContourValue) << std::endl;
536  }
537 
538 } // end namespace itk
539 
540 #endif