Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
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