1 /*===================================================================
3 The Medical Imaging Interaction Toolkit (MITK)
5 Copyright (c) German Cancer Research Center,
6 Division of Medical and Biological Informatics.
9 This software is distributed WITHOUT ANY WARRANTY; without
10 even the implied warranty of MERCHANTABILITY or FITNESS FOR
13 See LICENSE.txt or http://www.mitk.org for details.
15 ===================================================================*/
17 /*===================================================================
19 This file is based heavily on a corresponding ITK filter.
21 ===================================================================*/
23 #ifndef __itkContourExtractor2DImageFilter_txx
24 #define __itkContourExtractor2DImageFilter_txx
26 #include "itkConstShapedNeighborhoodIterator.h"
27 #include "itkConstShapedNeighborhoodIterator.h"
28 #include "itkContourExtractor2DImageFilter.h"
29 #include "itkProgressReporter.h"
30 #include "vcl_cmath.h"
35 template <class TInputImage>
36 ContourExtractor2DImageFilter<TInputImage>::ContourExtractor2DImageFilter()
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;
46 template <class TInputImage>
47 ContourExtractor2DImageFilter<TInputImage>::~ContourExtractor2DImageFilter()
51 template <class TInputImage>
52 void ContourExtractor2DImageFilter<TInputImage>::GenerateData()
54 // Make sure the structures for containing, looking up, and numbering the
55 // growing contours are empty and ready.
57 m_ContourStarts.clear();
58 m_ContourEnds.clear();
59 m_NumberOfContoursCreated = 0;
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
69 InputRegionType region = this->GetInput()->GetRequestedRegion();
70 typename InputRegionType::SizeType shrunkSize = region.GetSize();
73 InputRegionType shrunkRegion(region.GetIndex(), shrunkSize);
75 // Set up a progress reporter
76 ProgressReporter progress(this, 0, shrunkRegion.GetNumberOfPixels());
78 // A 1-pixel radius sets up a neighborhood with the following indices:
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);
98 for (it.GoToBegin(); !it.IsAtEnd(); ++it)
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:
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 // -- -- -- -- +- +- +- +-
111 // 8-- 9+- 10-+ 11++ 12-- 13+- 14-+ 15++
112 // -+ -+ -+ -+ ++ ++ ++ ++
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.
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
130 InputPixelType v0, v1, v2, v3;
135 InputIndexType index = it.GetIndex();
136 unsigned char squareCase = 0;
137 if (v0 > m_ContourValue)
139 if (v1 > m_ContourValue)
141 if (v2 > m_ContourValue)
143 if (v3 > m_ContourValue)
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)
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.
165 case 1: // top to left
166 this->AddSegment(TOP_, LEFT_);
168 case 2: // right to top
169 this->AddSegment(RIGHT_, TOP_);
171 case 3: // right to left
172 this->AddSegment(RIGHT_, LEFT_);
174 case 4: // left to bottom
175 this->AddSegment(LEFT_, BOTTOM_);
177 case 5: // top to bottom
178 this->AddSegment(TOP_, BOTTOM_);
181 if (m_VertexConnectHighPixels)
184 this->AddSegment(LEFT_, TOP_);
186 this->AddSegment(RIGHT_, BOTTOM_);
191 this->AddSegment(RIGHT_, TOP_);
193 this->AddSegment(LEFT_, BOTTOM_);
196 case 7: // right to bottom
197 this->AddSegment(RIGHT_, BOTTOM_);
199 case 8: // bottom to right
200 this->AddSegment(BOTTOM_, RIGHT_);
203 if (m_VertexConnectHighPixels)
206 this->AddSegment(TOP_, RIGHT_);
208 this->AddSegment(BOTTOM_, LEFT_);
213 this->AddSegment(TOP_, LEFT_);
215 this->AddSegment(BOTTOM_, RIGHT_);
218 case 10: // bottom to top
219 this->AddSegment(BOTTOM_, TOP_);
221 case 11: // bottom to left
222 this->AddSegment(BOTTOM_, LEFT_);
224 case 12: // left to right
225 this->AddSegment(LEFT_, RIGHT_);
227 case 13: // top to right
228 this->AddSegment(TOP_, RIGHT_);
230 case 14: // left to top
231 this->AddSegment(LEFT_, TOP_);
235 } // switch squareCase
236 progress.CompletedPixel();
239 // Now create the outputs paths from the deques we've been using.
242 m_ContourStarts.clear();
243 m_ContourEnds.clear();
244 m_NumberOfContoursCreated = 0;
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)
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");
264 itkAssertOrThrowMacro(((toOffset[0] == 0 && toOffset[1] == 1) || (toOffset[0] == 1 && toOffset[1] == 0)),
265 "toOffset has unexpected values");
268 (m_ContourValue - static_cast<InputRealType>(fromValue)) / (toValue - static_cast<InputRealType>(fromValue));
270 output[0] = fromIndex[0] + x * toOffset[0];
271 output[1] = fromIndex[1] + x * toOffset[1];
276 template <class TInputImage>
277 void ContourExtractor2DImageFilter<TInputImage>::AddSegment(VertexType from, VertexType to)
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
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);
293 if (newTail != m_ContourStarts.end() && newHead != m_ContourEnds.end())
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");
303 // We've closed a contour. Add the end point, and remove from the maps
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.
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)
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());
323 // Now remove 'tail' from the list and the maps because it has been
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
330 itkWarningMacro(<< "There should be exactly one entry in the hash for that endpoint, but there are "
333 m_Contours.erase(tail); // remove from the master list
335 // Now remove the old end of 'head' from the ends map and add
337 m_ContourEnds.erase(newHead);
338 m_ContourEnds.insert(VertexContourRefPair(head->back(), head));
342 // Copy head to the beginning of tail and remove
343 // head from everything.
344 tail->insert(tail->begin(), head->begin(), head->end());
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());
352 itkWarningMacro(<< "There should be exactly one entry in the hash for that endpoint, but there are "
355 m_Contours.erase(head); // remove from the master list
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));
364 else if (newTail == m_ContourStarts.end() && newHead == m_ContourEnds.end())
366 // No contours found: add a new one.
367 // Make it on the heap. It will be copied into m_Contours.
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);
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));
384 else if (newTail != m_ContourStarts.end() && newHead == m_ContourEnds.end())
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));
395 else if (newTail == m_ContourStarts.end() && newHead != m_ContourEnds.end())
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");
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));
408 template <class TInputImage>
409 void ContourExtractor2DImageFilter<TInputImage>::FillOutputs()
411 this->SetNumberOfIndexedOutputs(m_Contours.size());
413 for (auto it = m_Contours.begin(); it != m_Contours.end(); it++, i++)
415 OutputPathPointer output = this->GetOutput(i);
418 // Static cast is OK because we know PathSource will make its templated
420 output = static_cast<OutputPathType *>(this->MakeOutput(i).GetPointer());
421 this->SetNthOutput(i, output.GetPointer());
423 typename VertexListType::Pointer path = const_cast<VertexListType *>(output->GetVertexList());
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.
430 // Now put all the points from the contour deque into the path and
431 // mark output as modified
433 typedef typename ContourType::const_iterator ConstIteratorType;
434 if (m_ReverseContourOrientation)
436 ConstIteratorType itC = (*it).end();
440 path->push_back(*itC);
441 } while (itC != (*it).begin());
445 ConstIteratorType itC = (*it).begin();
446 while (itC != (*it).end())
448 path->push_back(*itC);
456 template <class TInputImage>
457 void ContourExtractor2DImageFilter<TInputImage>::SetRequestedRegion(const InputRegionType region)
459 itkDebugMacro("setting RequestedRegion to " << region);
460 m_UseCustomRegion = true;
461 if (this->m_RequestedRegion != region)
463 this->m_RequestedRegion = region;
468 template <class TInputImage>
469 void ContourExtractor2DImageFilter<TInputImage>::ClearRequestedRegion()
471 itkDebugMacro("Clearing RequestedRegion.");
472 if (this->m_UseCustomRegion == true)
474 this->m_UseCustomRegion = false;
479 template <class TInputImage>
480 void ContourExtractor2DImageFilter<TInputImage>::GenerateInputRequestedRegion() throw(InvalidRequestedRegionError)
482 InputImageType *input = const_cast<InputImageType *>(this->GetInput());
486 if (m_UseCustomRegion)
488 InputRegionType requestedRegion = m_RequestedRegion;
489 if (requestedRegion.Crop(input->GetLargestPossibleRegion()))
491 input->SetRequestedRegion(requestedRegion);
496 // Couldn't crop the region (requested region is outside the largest
497 // possible region). Throw an exception.
499 // store what we tried to request (prior to trying to crop)
500 input->SetRequestedRegion(requestedRegion);
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);
512 input->SetRequestedRegion(input->GetLargestPossibleRegion());
517 * Standard "PrintSelf" method
519 template <class TInputImage>
520 void ContourExtractor2DImageFilter<TInputImage>::PrintSelf(std::ostream &os, Indent indent) const
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)
530 os << indent << "Custom region: " << m_RequestedRegion << std::endl;
533 typedef typename NumericTraits<InputRealType>::PrintType InputRealPrintType;
535 os << indent << "Contour value: " << static_cast<InputRealPrintType>(m_ContourValue) << std::endl;
538 } // end namespace itk