Medical Imaging Interaction Toolkit  2018.4.99-ef453c4b
Medical Imaging Interaction Toolkit
mitkContourModelUtils.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 (DKFZ)
6 All rights reserved.
7 
8 Use of this source code is governed by a 3-clause BSD license that can be
9 found in the LICENSE file.
10 
11 ============================================================================*/
12 
13 #include <mitkContourModelUtils.h>
14 
16 #include <mitkLabelSetImage.h>
17 #include <mitkSurface.h>
18 #include <vtkImageStencil.h>
19 #include <vtkPointData.h>
20 #include <vtkPolyData.h>
21 #include <vtkPolyDataToImageStencil.h>
22 
24 {
25 }
26 
28 {
29 }
30 
32  Image *slice, ContourModel *contourIn3D, bool, bool)
33 {
34  if (nullptr == slice || nullptr == contourIn3D)
35  return nullptr;
36 
37  auto projectedContour = ContourModel::New();
38  projectedContour->Initialize(*contourIn3D);
39 
40  auto sliceGeometry = slice->GetGeometry();
41  auto numberOfTimesteps = static_cast<int>(contourIn3D->GetTimeSteps());
42 
43  for (decltype(numberOfTimesteps) t = 0; t < numberOfTimesteps; ++t)
44  {
45  auto iter = contourIn3D->Begin(t);
46  auto end = contourIn3D->End(t);
47 
48  while (iter != end)
49  {
50  const auto &currentPointIn3D = (*iter)->Coordinates;
51 
52  Point3D projectedPointIn2D;
53  projectedPointIn2D.Fill(0.0);
54 
55  sliceGeometry->WorldToIndex(currentPointIn3D, projectedPointIn2D);
56 
57  projectedContour->AddVertex(projectedPointIn2D, t);
58  ++iter;
59  }
60  }
61 
62  return projectedContour;
63 }
64 
66  const BaseGeometry *sliceGeometry, ContourModel *contourIn2D, bool)
67 {
68  if (nullptr == sliceGeometry || nullptr == contourIn2D)
69  return nullptr;
70 
71  auto worldContour = ContourModel::New();
72  worldContour->Initialize(*contourIn2D);
73 
74  auto numberOfTimesteps = static_cast<int>(contourIn2D->GetTimeSteps());
75 
76  for (decltype(numberOfTimesteps) t = 0; t < numberOfTimesteps; ++t)
77  {
78  auto iter = contourIn2D->Begin(t);
79  auto end = contourIn2D->End(t);
80 
81  while (iter != end)
82  {
83  const auto &currentPointIn2D = (*iter)->Coordinates;
84 
85  Point3D worldPointIn3D;
86  worldPointIn3D.Fill(0.0);
87 
88  sliceGeometry->IndexToWorld(currentPointIn2D, worldPointIn3D);
89 
90  worldContour->AddVertex(worldPointIn3D, t);
91  ++iter;
92  }
93  }
94 
95  return worldContour;
96 }
97 
99  ContourModel *projectedContour, Image *sliceImage, Image::Pointer workingImage, int paintingPixelValue)
100 {
101  FillContourInSlice(projectedContour, 0, sliceImage, workingImage, paintingPixelValue);
102 }
103 
105  ContourModel *projectedContour, unsigned int t, Image *sliceImage, Image::Pointer workingImage, int paintingPixelValue)
106 {
107  auto contourModelFilter = mitk::ContourModelToSurfaceFilter::New();
108  contourModelFilter->SetInput(projectedContour);
109  contourModelFilter->Update();
110 
111  auto surface = mitk::Surface::New();
112  surface = contourModelFilter->GetOutput();
113 
114  if (nullptr == surface->GetVtkPolyData(t))
115  {
116  MITK_WARN << "Could not create surface from contour model.";
117  return;
118  }
119 
120  auto surface2D = vtkSmartPointer<vtkPolyData>::New();
121  surface2D->SetPoints(surface->GetVtkPolyData(t)->GetPoints());
122  surface2D->SetLines(surface->GetVtkPolyData(t)->GetLines());
123 
124  auto image = vtkSmartPointer<vtkImageData>::New();
125  image->DeepCopy(sliceImage->GetVtkImageData());
126 
127  const double FOREGROUND_VALUE = 255.0;
128  const double BACKGROUND_VALUE = 0.0;
129 
130  vtkIdType count = image->GetNumberOfPoints();
131  for (decltype(count) i = 0; i < count; ++i)
132  image->GetPointData()->GetScalars()->SetTuple1(i, FOREGROUND_VALUE);
133 
134  auto polyDataToImageStencil = vtkSmartPointer<vtkPolyDataToImageStencil>::New();
135 
136  // Set a minimal tolerance, so that clipped pixels will be added to contour as well.
137  polyDataToImageStencil->SetTolerance(mitk::eps);
138  polyDataToImageStencil->SetInputData(surface2D);
139  polyDataToImageStencil->Update();
140 
141  auto imageStencil = vtkSmartPointer<vtkImageStencil>::New();
142 
143  imageStencil->SetInputData(image);
144  imageStencil->SetStencilConnection(polyDataToImageStencil->GetOutputPort());
145  imageStencil->ReverseStencilOff();
146  imageStencil->SetBackgroundValue(BACKGROUND_VALUE);
147  imageStencil->Update();
148 
149  vtkSmartPointer<vtkImageData> filledImage = imageStencil->GetOutput();
150  vtkSmartPointer<vtkImageData> resultImage = sliceImage->GetVtkImageData();
151  FillSliceInSlice(filledImage, resultImage, workingImage, paintingPixelValue);
152 
153  sliceImage->SetVolume(resultImage->GetScalarPointer());
154 }
155 
157  vtkSmartPointer<vtkImageData> filledImage, vtkSmartPointer<vtkImageData> resultImage, mitk::Image::Pointer image, int paintingPixelValue)
158 {
159  auto labelImage = dynamic_cast<LabelSetImage *>(image.GetPointer());
160  auto numberOfPoints = filledImage->GetNumberOfPoints();
161 
162  if (nullptr == labelImage)
163  {
164  for (decltype(numberOfPoints) i = 0; i < numberOfPoints; ++i)
165  {
166  if (1 < filledImage->GetPointData()->GetScalars()->GetTuple1(i))
167  resultImage->GetPointData()->GetScalars()->SetTuple1(i, paintingPixelValue);
168  }
169  }
170  else
171  {
172  auto backgroundValue = labelImage->GetExteriorLabel()->GetValue();
173 
174  if (paintingPixelValue != backgroundValue)
175  {
176  for (decltype(numberOfPoints) i = 0; i < numberOfPoints; ++i)
177  {
178  if (1 < filledImage->GetPointData()->GetScalars()->GetTuple1(i))
179  {
180  auto existingValue = resultImage->GetPointData()->GetScalars()->GetTuple1(i);
181 
182  if (!labelImage->GetLabel(existingValue, labelImage->GetActiveLayer())->GetLocked())
183  resultImage->GetPointData()->GetScalars()->SetTuple1(i, paintingPixelValue);
184  }
185  }
186  }
187  else
188  {
189  auto activePixelValue = labelImage->GetActiveLabel(labelImage->GetActiveLayer())->GetValue();
190 
191  for (decltype(numberOfPoints) i = 0; i < numberOfPoints; ++i)
192  {
193  if (1 < filledImage->GetPointData()->GetScalars()->GetTuple1(i))
194  {
195  if (resultImage->GetPointData()->GetScalars()->GetTuple1(i) == activePixelValue)
196  resultImage->GetPointData()->GetScalars()->SetTuple1(i, paintingPixelValue);
197  }
198  }
199  }
200  }
201 }
202 
204 {
205  if (nullptr == contour)
206  return nullptr;
207 
208  auto resultContour = ContourModel::New();
209  resultContour->Expand(t + 1);
210 
211  std::for_each(contour->Begin(), contour->End(), [&resultContour, t](ContourElement::VertexType *vertex) {
212  resultContour->AddVertex(vertex, t);
213  });
214 
215  return resultContour;
216 }
static ContourModel::Pointer BackProjectContourFrom2DSlice(const BaseGeometry *sliceGeometry, ContourModel *contourIn2D, bool correctionForIpSegmentation=false)
Projects a slice index coordinates of a contour back into world coordinates.
ContourModel is a structure of linked vertices defining a contour in 3D space. The vertices are store...
void IndexToWorld(const mitk::Vector3D &vec_units, mitk::Vector3D &vec_mm) const
Convert (continuous or discrete) index coordinates of a vector vec_units to world coordinates (in mm)...
static void FillContourInSlice(ContourModel *projectedContour, Image *sliceImage, mitk::Image::Pointer workingImage, int paintingPixelValue=1)
Fill a contour in a 2D slice with a specified pixel value at time step 0.
virtual vtkImageData * GetVtkImageData(int t=0, int n=0)
Get a volume at a specific time t of channel n as a vtkImageData.
Definition: mitkImage.cpp:217
mitk::Label * GetExteriorLabel()
Gets the mitk::Label which is used as default exterior label.
VertexIterator Begin(int timestep=0) const
Returns a const VertexIterator at the start element of the contour.
virtual bool SetVolume(const void *data, int t=0, int n=0)
Set data as volume at time t in channel n. It is in the responsibility of the caller to ensure that t...
Definition: mitkImage.cpp:669
#define MITK_WARN
Definition: mitkLogMacros.h:19
Image class for storing images.
Definition: mitkImage.h:72
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.
mitk::Image::Pointer image
PixelType GetValue() const
Definition: mitkLabel.cpp:169
VertexIterator End(int timestep=0) const
Returns a const VertexIterator at the end element of the contour.
LabelSetImage class for handling labels and layers in a segmentation session.
MITKCORE_EXPORT const ScalarType eps
unsigned int GetTimeSteps() const
Get the number of time steps from the TimeGeometry As the base data has not a data vector given by it...
Definition: mitkBaseData.h:355
static Pointer New()
static void FillSliceInSlice(vtkSmartPointer< vtkImageData > filledImage, vtkSmartPointer< vtkImageData > resultImage, mitk::Image::Pointer image, int paintingPixelValue)
Fills a image (filledImage) into another image (resultImage) by incorporating the rules of LabelSet-I...
mitk::BaseGeometry * GetGeometry(int t=0) const
Return the geometry, which is a TimeGeometry, of the data as non-const pointer.
Definition: mitkBaseData.h:138
static Pointer New()
static ContourModel::Pointer MoveZerothContourTimeStep(const ContourModel *contour, unsigned int timeStep)
Move the contour in time step 0 to to a new contour model at the given time step. ...
Represents a single vertex of contour.
BaseGeometry Describes the geometry of a data object.