Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkCorrectorAlgorithm.cpp
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 #include "mitkCorrectorAlgorithm.h"
18 #include "mitkContourUtils.h"
19 #include "mitkITKImageImport.h"
20 #include "mitkImageAccessByItk.h"
21 #include "mitkImageCast.h"
22 #include "mitkImageDataItem.h"
23 #include "mitkLegacyAdaptors.h"
24 
25 #include <mitkContourModelUtils.h>
26 
27 #include "itkCastImageFilter.h"
28 #include "itkImageDuplicator.h"
29 #include "itkImageRegionIterator.h"
30 
32 {
33 }
34 
36 {
37 }
38 
39 template <typename TPixel, unsigned int VDimensions>
41  itk::Image<TPixel, VDimensions> *reference,
42  mitk::Image::Pointer target,
44 {
45  typedef itk::Image<mitk::CorrectorAlgorithm::DefaultSegmentationDataType, 2> InputImageType;
46  typedef itk::Image<TPixel, 2> OutputImageType;
47  typedef itk::CastImageFilter<InputImageType, OutputImageType> CastImageFilterType;
48 
49  typename CastImageFilterType::Pointer castImageFilter = CastImageFilterType::New();
50  castImageFilter->SetInput(segmentationPixelTypeImage);
51  castImageFilter->Update();
52  typename OutputImageType::Pointer tempItkImage = castImageFilter->GetOutput();
53  tempItkImage->DisconnectPipeline();
54  mitk::CastToMitkImage(tempItkImage, target);
55 }
56 
58 {
59  Image::Pointer inputImage = const_cast<Image *>(ImageToImageFilter::GetInput(0));
60 
61  if (inputImage.IsNull() || inputImage->GetDimension() != 2)
62  {
63  itkExceptionMacro("CorrectorAlgorithm needs a 2D image as input.");
64  }
65 
66  if (m_Contour.IsNull())
67  {
68  itkExceptionMacro("CorrectorAlgorithm needs a Contour object as input.");
69  }
70 
71  // copy the input (since m_WorkingImage will be changed later)
72  m_WorkingImage = inputImage;
73 
74  TimeGeometry::Pointer originalGeometry = NULL;
75 
76  if (inputImage->GetTimeGeometry())
77  {
78  originalGeometry = inputImage->GetTimeGeometry()->Clone();
79  m_WorkingImage->SetTimeGeometry(originalGeometry);
80  }
81  else
82  {
83  itkExceptionMacro("Original image does not have a 'Time sliced geometry'! Cannot copy.");
84  }
85 
86  Image::Pointer temporarySlice;
87 
88  // Convert to DefaultSegmentationDataType (because TobiasHeimannCorrectionAlgorithm relys on that data type)
89  {
91  CastToItkImage(m_WorkingImage, correctPixelTypeImage);
92  assert(correctPixelTypeImage.IsNotNull());
93 
94  // possible bug in CastToItkImage ?
95  // direction maxtrix is wrong/broken/not working after CastToItkImage, leading to a failed assertion in
96  // mitk/Core/DataStructures/mitkSlicedGeometry3D.cpp, 479:
97  // virtual void mitk::SlicedGeometry3D::SetSpacing(const mitk::Vector3D&): Assertion `aSpacing[0]>0 && aSpacing[1]>0
98  // && aSpacing[2]>0' failed
99  // solution here: we overwrite it with an unity matrix
100  itk::Image<DefaultSegmentationDataType, 2>::DirectionType imageDirection;
101 
102  imageDirection.SetIdentity();
103  // correctPixelTypeImage->SetDirection(imageDirection);
104 
105  temporarySlice = this->GetOutput();
106  // temporarySlice = ImportItkImage( correctPixelTypeImage );
107  // m_FillColor = 1;
108  m_EraseColor = 0;
109 
110  ImprovedHeimannCorrectionAlgorithm(correctPixelTypeImage);
111 
112  // this is suboptimal, needs to be kept synchronous to DefaultSegmentationDataType
113  if (inputImage->GetChannelDescriptor().GetPixelType().GetComponentType() == itk::ImageIOBase::USHORT)
114  { // the cast at the beginning did not copy the data
115  CastToMitkImage(correctPixelTypeImage, temporarySlice);
116  }
117  else
118  { // it did copy the data and cast the pixel type
119  AccessByItk_n(m_WorkingImage, ConvertBackToCorrectPixelType, (temporarySlice, correctPixelTypeImage));
120  }
121  }
122  temporarySlice->SetTimeGeometry(originalGeometry);
123 }
124 
125 template <typename ScalarType>
126 itk::Index<2> mitk::CorrectorAlgorithm::ensureIndexInImage(ScalarType i0, ScalarType i1)
127 {
128  itk::Index<2> toReturn;
129 
130  itk::Size<5> size = m_WorkingImage->GetLargestPossibleRegion().GetSize();
131 
132  toReturn[0] = std::min((ScalarType)(size[0] - 1), std::max((ScalarType)0.0, i0));
133  toReturn[1] = std::min((ScalarType)(size[1] - 1), std::max((ScalarType)0.0, i1));
134 
135  return toReturn;
136 }
137 
140 {
167  ContourModel::Pointer projectedContour =
168  mitk::ContourModelUtils::ProjectContourTo2DSlice(m_WorkingImage, m_Contour, true, false);
169 
170  bool firstPointIsFillingColor = false;
171 
172  if (projectedContour.IsNull() || projectedContour->GetNumberOfVertices() < 2)
173  {
174  return false;
175  }
176 
177  // Read the first point of the contour
178  ContourModel::VertexIterator contourIter = projectedContour->Begin();
179  if (contourIter == projectedContour->End())
180  return false;
181  itk::Index<2> previousIndex;
182 
183  previousIndex = ensureIndexInImage((*contourIter)->Coordinates[0], (*contourIter)->Coordinates[1]);
184  ++contourIter;
185 
186  int currentColor = (pic->GetPixel(previousIndex) == m_FillColor);
187  firstPointIsFillingColor = currentColor;
188  TSegData currentSegment;
189  int countOfSegments = 1;
190 
191  bool firstSegment = true;
192  ContourModel::VertexIterator contourEnd = projectedContour->End();
193  for (; contourIter != contourEnd; ++contourIter)
194  {
195  // Get current point
196  itk::Index<2> currentIndex;
197  currentIndex = ensureIndexInImage((*contourIter)->Coordinates[0] + 0.5, (*contourIter)->Coordinates[1] + 0.5);
198 
199  // Calculate length and slope
200  double slopeX = currentIndex[0] - previousIndex[0];
201  double slopeY = currentIndex[1] - previousIndex[1];
202  double length = std::sqrt(slopeX * slopeX + slopeY * slopeY);
203  double deltaX = slopeX / length;
204  double deltaY = slopeY / length;
205 
206  for (double i = 0; i <= length && length > 0; i += 1)
207  {
208  itk::Index<2> temporaryIndex;
209  temporaryIndex = ensureIndexInImage(previousIndex[0] + deltaX * i, previousIndex[1] + deltaY * i);
210 
211  if (!pic->GetLargestPossibleRegion().IsInside(temporaryIndex))
212  continue;
213  if ((pic->GetPixel(temporaryIndex) == m_FillColor) != currentColor)
214  {
215  currentSegment.points.push_back(temporaryIndex);
216  if (!firstSegment)
217  {
218  ModifySegment(currentSegment, pic);
219  }
220  else
221  {
222  firstSegment = false;
223  }
224  currentSegment = TSegData();
225  ++countOfSegments;
226  currentColor = (pic->GetPixel(temporaryIndex) == m_FillColor);
227  }
228  currentSegment.points.push_back(temporaryIndex);
229  }
230  previousIndex = currentIndex;
231  }
232 
233  return true;
234 }
235 
236 void mitk::CorrectorAlgorithm::ColorSegment(
237  const mitk::CorrectorAlgorithm::TSegData &segment,
239 {
240  int colorMode = (pic->GetPixel(segment.points[0]) == m_FillColor);
241  int color = 0;
242  if (colorMode)
243  color = m_EraseColor;
244  else
245  color = m_FillColor;
246 
247  std::vector<itk::Index<2>>::const_iterator indexIterator;
248  std::vector<itk::Index<2>>::const_iterator indexEnd;
249 
250  indexIterator = segment.points.begin();
251  indexEnd = segment.points.end();
252 
253  for (; indexIterator != indexEnd; ++indexIterator)
254  {
255  pic->SetPixel(*indexIterator, color);
256  }
257 }
258 
261 {
262  typedef itk::Image<mitk::CorrectorAlgorithm::DefaultSegmentationDataType, 2> ItkImageType;
263  typedef itk::ImageDuplicator<ItkImageType> DuplicatorType;
265  duplicator->SetInputImage(pic);
266  duplicator->Update();
267 
268  return duplicator->GetOutput();
269 }
270 
271 itk::Index<2> mitk::CorrectorAlgorithm::GetFirstPoint(
272  const mitk::CorrectorAlgorithm::TSegData &segment,
274 {
275  int colorMode = (pic->GetPixel(segment.points[0]) == m_FillColor);
276 
277  std::vector<itk::Index<2>>::const_iterator indexIterator;
278  std::vector<itk::Index<2>>::const_iterator indexEnd;
279 
280  indexIterator = segment.points.begin();
281  indexEnd = segment.points.end();
282 
283  itk::Index<2> index;
284 
285  for (; indexIterator != indexEnd; ++indexIterator)
286  {
287  for (int xOffset = -1; xOffset < 2; ++xOffset)
288  {
289  for (int yOffset = -1; yOffset < 2; ++yOffset)
290  {
291  index = ensureIndexInImage((*indexIterator)[0] - xOffset, (*indexIterator)[1] - yOffset);
292 
293  if ((pic->GetPixel(index) == m_FillColor) != colorMode)
294  {
295  return index;
296  }
297  }
298  }
299  }
300  mitkThrow() << "No Starting point is found next to the curve.";
301 }
302 
303 std::vector<itk::Index<2>> mitk::CorrectorAlgorithm::FindSeedPoints(
304  const mitk::CorrectorAlgorithm::TSegData &segment,
306 {
307  typedef itk::Image<mitk::CorrectorAlgorithm::DefaultSegmentationDataType, 2> ItkImageType;
309 
310  std::vector<itk::Index<2>> seedPoints;
311 
312  try
313  {
314  itk::Index<2> firstPoint = GetFirstPoint(segment, pic);
315  seedPoints.push_back(firstPoint);
316  }
317  catch (mitk::Exception e)
318  {
319  return seedPoints;
320  }
321 
322  if (segment.points.size() < 4)
323  return seedPoints;
324 
325  std::vector<itk::Index<2>>::const_iterator indexIterator;
326  std::vector<itk::Index<2>>::const_iterator indexEnd;
327 
328  indexIterator = segment.points.begin();
329  indexEnd = segment.points.end();
330 
331  ItkImagePointerType listOfPoints = CloneImage(pic);
332  listOfPoints->FillBuffer(0);
333  listOfPoints->SetPixel(seedPoints[0], 1);
334  for (; indexIterator != indexEnd; ++indexIterator)
335  {
336  listOfPoints->SetPixel(*indexIterator, 2);
337  }
338  indexIterator = segment.points.begin();
339  indexIterator++;
340  indexIterator++;
341  indexEnd--;
342  indexEnd--;
343  for (; indexIterator != indexEnd; ++indexIterator)
344  {
345  bool pointFound = true;
346  while (pointFound)
347  {
348  pointFound = false;
349  itk::Index<2> index;
350  itk::Index<2> index2;
351  for (int xOffset = -1; xOffset < 2; ++xOffset)
352  {
353  for (int yOffset = -1; yOffset < 2; ++yOffset)
354  {
355  index = ensureIndexInImage((*indexIterator)[0] - xOffset, (*indexIterator)[1] - yOffset);
356  index2 = index;
357 
358  if (listOfPoints->GetPixel(index2) > 0)
359  continue;
360 
361  index[0] = index[0] - 1;
362  index = ensureIndexInImage(index[0], index[1]);
363  if (listOfPoints->GetPixel(index) == 1)
364  {
365  pointFound = true;
366  seedPoints.push_back(index2);
367  listOfPoints->SetPixel(index2, 1);
368  continue;
369  }
370 
371  index[0] = index[0] + 2;
372  index = ensureIndexInImage(index[0], index[1]);
373  if (listOfPoints->GetPixel(index) == 1)
374  {
375  pointFound = true;
376  seedPoints.push_back(index2);
377  listOfPoints->SetPixel(index2, 1);
378  continue;
379  }
380 
381  index[0] = index[0] - 1;
382  index[1] = index[1] - 1;
383  index = ensureIndexInImage(index[0], index[1]);
384  if (listOfPoints->GetPixel(index) == 1)
385  {
386  pointFound = true;
387  seedPoints.push_back(index2);
388  listOfPoints->SetPixel(index2, 1);
389  continue;
390  }
391 
392  index[1] = index[1] + 2;
393  index = ensureIndexInImage(index[0], index[1]);
394  if (listOfPoints->GetPixel(index) == 1)
395  {
396  pointFound = true;
397  seedPoints.push_back(index2);
398  listOfPoints->SetPixel(index2, 1);
399  continue;
400  }
401  }
402  }
403  }
404  }
405  return seedPoints;
406 }
407 
408 int mitk::CorrectorAlgorithm::FillRegion(
409  const std::vector<itk::Index<2>> &seedPoints,
411 {
412  int numberOfPixel = 0;
413  int mode = (pic->GetPixel(seedPoints[0]) == m_FillColor);
414  int drawColor = m_FillColor;
415  if (mode)
416  {
417  drawColor = m_EraseColor;
418  }
419 
420  std::vector<itk::Index<2>> workPoints;
421  workPoints = seedPoints;
422  // workPoints.push_back(seedPoints[0]);
423  while (workPoints.size() > 0)
424  {
425  itk::Index<2> currentIndex = workPoints.back();
426  workPoints.pop_back();
427  if ((pic->GetPixel(currentIndex) == m_FillColor) == mode)
428  ++numberOfPixel;
429  pic->SetPixel(currentIndex, drawColor);
430 
431  currentIndex = ensureIndexInImage(currentIndex[0] - 1, currentIndex[1]);
432  if (pic->GetLargestPossibleRegion().IsInside(currentIndex) && (pic->GetPixel(currentIndex) == m_FillColor) == mode)
433  workPoints.push_back(currentIndex);
434 
435  currentIndex = ensureIndexInImage(currentIndex[0] + 2, currentIndex[1]);
436  if (pic->GetLargestPossibleRegion().IsInside(currentIndex) && (pic->GetPixel(currentIndex) == m_FillColor) == mode)
437  workPoints.push_back(currentIndex);
438 
439  currentIndex = ensureIndexInImage(currentIndex[0] - 1, currentIndex[1] - 1);
440 
441  if (pic->GetLargestPossibleRegion().IsInside(currentIndex) && (pic->GetPixel(currentIndex) == m_FillColor) == mode)
442  workPoints.push_back(currentIndex);
443 
444  currentIndex = ensureIndexInImage(currentIndex[0], currentIndex[1] + 2);
445  if (pic->GetLargestPossibleRegion().IsInside(currentIndex) && (pic->GetPixel(currentIndex) == m_FillColor) == mode)
446  workPoints.push_back(currentIndex);
447  }
448 
449  return numberOfPixel;
450 }
451 
452 void mitk::CorrectorAlgorithm::OverwriteImage(
455 {
456  typedef itk::Image<mitk::CorrectorAlgorithm::DefaultSegmentationDataType, 2> ItkImageType;
457  typedef itk::ImageRegionIterator<ItkImageType> ImageIteratorType;
458 
459  ImageIteratorType sourceIter(source, source->GetLargestPossibleRegion());
460  ImageIteratorType targetIter(target, target->GetLargestPossibleRegion());
461  while (!sourceIter.IsAtEnd())
462  {
463  targetIter.Set(sourceIter.Get());
464  ++sourceIter;
465  ++targetIter;
466  }
467 }
468 
471 {
472  typedef itk::Image<DefaultSegmentationDataType, 2>::Pointer ItkImagePointerType;
473 
474  ItkImagePointerType firstSideImage = CloneImage(pic);
475  ColorSegment(segment, firstSideImage);
476  ItkImagePointerType secondSideImage = CloneImage(firstSideImage);
477 
478  std::vector<itk::Index<2>> seedPoints = FindSeedPoints(segment, firstSideImage);
479  if (seedPoints.size() < 1)
480  return false;
481  int firstSidePixel = FillRegion(seedPoints, firstSideImage);
482 
483  std::vector<itk::Index<2>> secondSeedPoints = FindSeedPoints(segment, firstSideImage);
484  if (secondSeedPoints.size() < 1)
485  return false;
486  int secondSidePixel = FillRegion(secondSeedPoints, secondSideImage);
487 
488  if (firstSidePixel < secondSidePixel)
489  {
490  OverwriteImage(firstSideImage, pic);
491  }
492  else
493  {
494  OverwriteImage(secondSideImage, pic);
495  }
496  return true;
497 }
itk::SmartPointer< Self > Pointer
double ScalarType
mitk::ContourElement::VertexIterator VertexIterator
void ConvertBackToCorrectPixelType(itk::Image< TPixel, VDimensions > *reference, mitk::Image::Pointer target, itk::Image< mitk::CorrectorAlgorithm::DefaultSegmentationDataType, 2 >::Pointer segmentationPixelTypeImage)
#define AccessByItk_n(mitkImage, itkImageTypeFunction, va_tuple)
Access a MITK image by an ITK image with one or more parameters.
bool ModifySegment(const TSegData &segment, itk::Image< DefaultSegmentationDataType, 2 >::Pointer pic)
bool ImprovedHeimannCorrectionAlgorithm(itk::Image< DefaultSegmentationDataType, 2 >::Pointer pic)
itk::Image< double, 3 > InputImageType
An object of this class represents an exception of MITK. Please don't instantiate exceptions manually...
Definition: mitkException.h:49
Calculated difference image.
#define mitkThrow()
Image class for storing images.
Definition: mitkImage.h:76
static ContourModel::Pointer ProjectContourTo2DSlice(Image *slice, ContourModel *contourIn3D, bool correctionForIpSegmentation, bool constrainToInside)
Projects a contour onto an image point by point. Converts from world to index coordinates.
static T max(T x, T y)
Definition: svm.cpp:70
Superclass of all classes having one or more Images as input and generating Images as output...
static T min(T x, T y)
Definition: svm.cpp:67
InputImageType * GetInput(void)
void CastToMitkImage(const itk::SmartPointer< ItkOutputImageType > &itkimage, itk::SmartPointer< mitk::Image > &mitkoutputimage)
Cast an itk::Image (with a specific type) to an mitk::Image.
Definition: mitkImageCast.h:78
void MITKCORE_EXPORT CastToItkImage(const mitk::Image *mitkImage, itk::SmartPointer< ItkOutputImageType > &itkOutputImage)
Cast an mitk::Image to an itk::Image with a specific type.
std::vector< itk::Index< 2 > > points
virtual void GenerateData() override
A version of GenerateData() specific for image processing filters.
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.