Medical Imaging Interaction Toolkit  2018.4.99-6a3ea89d
Medical Imaging Interaction Toolkit
mitkComputeContourSetNormalsFilter.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 
14 
15 #include "mitkIOUtil.h"
17 
19  : m_SegmentationBinaryImage(nullptr),
20  m_MaxSpacing(5),
21  m_NegativeNormalCounter(0),
22  m_PositiveNormalCounter(0),
23  m_UseProgressBar(false),
24  m_ProgressStepSize(1)
25 {
27  this->SetNthOutput(0, output.GetPointer());
28 }
29 
31 {
32 }
33 
35 {
36  unsigned int numberOfInputs = this->GetNumberOfIndexedInputs();
37 
38  // Iterating over each input
39  for (unsigned int i = 0; i < numberOfInputs; i++)
40  {
41  // Getting the inputs polydata and polygons
42  auto *currentSurface = this->GetInput(i);
43  vtkPolyData *polyData = currentSurface->GetVtkPolyData();
44 
45  vtkSmartPointer<vtkCellArray> existingPolys = polyData->GetPolys();
46 
47  vtkSmartPointer<vtkPoints> existingPoints = polyData->GetPoints();
48 
49  existingPolys->InitTraversal();
50 
51  vtkIdType *cell(nullptr);
52  vtkIdType cellSize(0);
53 
54  // The array that contains all the vertex normals of the current polygon
55  vtkSmartPointer<vtkDoubleArray> normals = vtkSmartPointer<vtkDoubleArray>::New();
56  normals->SetNumberOfComponents(3);
57  normals->SetNumberOfTuples(polyData->GetNumberOfPoints());
58 
59  // If the current contour is an inner contour then the direction is -1
60  // A contour lies inside another one if the pixel values in the direction of the normal is 1
61  m_NegativeNormalCounter = 0;
62  m_PositiveNormalCounter = 0;
63  vtkIdType offSet(0);
64 
65  // Iterating over each polygon
66  for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);)
67  {
68  if (cellSize < 3)
69  continue;
70 
71  // First we calculate the current polygon's normal
72  double polygonNormal[3] = {0.0};
73 
74  double p1[3];
75  double p2[3];
76 
77  double v1[3];
78  double v2[3];
79 
80  existingPoints->GetPoint(cell[0], p1);
81  unsigned int index = cellSize * 0.5;
82  existingPoints->GetPoint(cell[index], p2);
83 
84  v1[0] = p2[0] - p1[0];
85  v1[1] = p2[1] - p1[1];
86  v1[2] = p2[2] - p1[2];
87 
88  for (vtkIdType k = 2; k < cellSize; k++)
89  {
90  index = cellSize * 0.25;
91  existingPoints->GetPoint(cell[index], p1);
92  index = cellSize * 0.75;
93  existingPoints->GetPoint(cell[index], p2);
94 
95  v2[0] = p2[0] - p1[0];
96  v2[1] = p2[1] - p1[1];
97  v2[2] = p2[2] - p1[2];
98 
99  vtkMath::Cross(v1, v2, polygonNormal);
100  if (vtkMath::Norm(polygonNormal) != 0)
101  break;
102  }
103 
104  vtkMath::Normalize(polygonNormal);
105 
106  // Now we start computing the normal for each vertex
107 
108  double vertexNormalTemp[3];
109  existingPoints->GetPoint(cell[0], p1);
110  existingPoints->GetPoint(cell[1], p2);
111 
112  v1[0] = p2[0] - p1[0];
113  v1[1] = p2[1] - p1[1];
114  v1[2] = p2[2] - p1[2];
115 
116  vtkMath::Cross(v1, polygonNormal, vertexNormalTemp);
117 
118  vtkMath::Normalize(vertexNormalTemp);
119 
120  double vertexNormal[3];
121 
122  for (vtkIdType j = 0; j < cellSize - 2; j++)
123  {
124  existingPoints->GetPoint(cell[j + 1], p1);
125  existingPoints->GetPoint(cell[j + 2], p2);
126 
127  v1[0] = p2[0] - p1[0];
128  v1[1] = p2[1] - p1[1];
129  v1[2] = p2[2] - p1[2];
130 
131  vtkMath::Cross(v1, polygonNormal, vertexNormal);
132 
133  vtkMath::Normalize(vertexNormal);
134 
135  double finalNormal[3];
136 
137  finalNormal[0] = (vertexNormal[0] + vertexNormalTemp[0]) * 0.5;
138  finalNormal[1] = (vertexNormal[1] + vertexNormalTemp[1]) * 0.5;
139  finalNormal[2] = (vertexNormal[2] + vertexNormalTemp[2]) * 0.5;
140  vtkMath::Normalize(finalNormal);
141 
142  // Here we determine the direction of the normal
143  if (m_SegmentationBinaryImage)
144  {
145  Point3D worldCoord;
146  worldCoord[0] = p1[0] + finalNormal[0] * m_MaxSpacing;
147  worldCoord[1] = p1[1] + finalNormal[1] * m_MaxSpacing;
148  worldCoord[2] = p1[2] + finalNormal[2] * m_MaxSpacing;
149 
150  double val = 0.0;
151 
152  itk::Index<3> idx;
153  m_SegmentationBinaryImage->GetGeometry()->WorldToIndex(worldCoord, idx);
154  try
155  {
156  if (m_SegmentationBinaryImage->GetImageDescriptor()
157  ->GetChannelDescriptor()
158  .GetPixelType()
159  .GetComponentType() == itk::ImageIOBase::UCHAR)
160  {
161  mitk::ImagePixelReadAccessor<unsigned char> readAccess(m_SegmentationBinaryImage);
162  val = readAccess.GetPixelByIndexSafe(idx);
163  }
164  else if (m_SegmentationBinaryImage->GetImageDescriptor()
165  ->GetChannelDescriptor()
166  .GetPixelType()
167  .GetComponentType() == itk::ImageIOBase::USHORT)
168  {
169  mitk::ImagePixelReadAccessor<unsigned short> readAccess(m_SegmentationBinaryImage);
170  val = readAccess.GetPixelByIndexSafe(idx);
171  }
172  }
173  catch (const mitk::Exception &e)
174  {
175  // If value is outside the image's region ignore it
176  MITK_WARN << e.what();
177  }
178 
179  if (val == 0.0)
180  {
181  // MITK_INFO << "val equals zero.";
182  ++m_PositiveNormalCounter;
183  }
184  else
185  {
186  // MITK_INFO << "val does not equal zero.";
187  ++m_NegativeNormalCounter;
188  }
189  }
190 
191  vertexNormalTemp[0] = vertexNormal[0];
192  vertexNormalTemp[1] = vertexNormal[1];
193  vertexNormalTemp[2] = vertexNormal[2];
194 
195  vtkIdType id = cell[j + 1];
196  normals->SetTuple(id, finalNormal);
197  }
198 
199  existingPoints->GetPoint(cell[0], p1);
200  existingPoints->GetPoint(cell[1], p2);
201 
202  v1[0] = p2[0] - p1[0];
203  v1[1] = p2[1] - p1[1];
204  v1[2] = p2[2] - p1[2];
205 
206  vtkMath::Cross(v1, polygonNormal, vertexNormal);
207 
208  vtkMath::Normalize(vertexNormal);
209 
210  vertexNormal[0] = (vertexNormal[0] + vertexNormalTemp[0]) * 0.5;
211  vertexNormal[1] = (vertexNormal[1] + vertexNormalTemp[1]) * 0.5;
212  vertexNormal[2] = (vertexNormal[2] + vertexNormalTemp[2]) * 0.5;
213  vtkMath::Normalize(vertexNormal);
214 
215  vtkIdType id = cell[0];
216  normals->SetTuple(id, vertexNormal);
217  id = cell[cellSize - 1];
218  normals->SetTuple(id, vertexNormal);
219 
220  if (m_NegativeNormalCounter > m_PositiveNormalCounter)
221  {
222  for (vtkIdType n = 0; n < cellSize; n++)
223  {
224  double normal[3];
225  normals->GetTuple(offSet + n, normal);
226  normal[0] = (-1) * normal[0];
227  normal[1] = (-1) * normal[1];
228  normal[2] = (-1) * normal[2];
229  normals->SetTuple(offSet + n, normal);
230  }
231  }
232 
233  m_NegativeNormalCounter = 0;
234  m_PositiveNormalCounter = 0;
235  offSet += cellSize;
236 
237  } // end for all cells
238 
239  Surface::Pointer surface = this->GetOutput(i);
240  surface->GetVtkPolyData()->GetCellData()->SetNormals(normals);
241  } // end for all inputs
242 
243  // Setting progressbar
244  if (this->m_UseProgressBar)
245  mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize);
246 }
247 
249 {
250  // Just for debugging:
251  vtkSmartPointer<vtkPolyData> newPolyData = vtkSmartPointer<vtkPolyData>::New();
252  vtkSmartPointer<vtkCellArray> newLines = vtkSmartPointer<vtkCellArray>::New();
253  vtkSmartPointer<vtkPoints> newPoints = vtkSmartPointer<vtkPoints>::New();
254  unsigned int idCounter(0);
255  // Debug end
256 
257  for (unsigned int i = 0; i < this->GetNumberOfIndexedOutputs(); i++)
258  {
259  auto *currentSurface = this->GetOutput(i);
260  vtkPolyData *polyData = currentSurface->GetVtkPolyData();
261 
262  vtkSmartPointer<vtkDoubleArray> currentCellNormals =
263  vtkDoubleArray::SafeDownCast(polyData->GetCellData()->GetNormals());
264 
265  vtkSmartPointer<vtkCellArray> existingPolys = polyData->GetPolys();
266 
267  vtkSmartPointer<vtkPoints> existingPoints = polyData->GetPoints();
268 
269  existingPolys->InitTraversal();
270 
271  vtkIdType *cell(nullptr);
272  vtkIdType cellSize(0);
273 
274  for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);)
275  {
276  for (vtkIdType j = 0; j < cellSize; j++)
277  {
278  double currentNormal[3];
279  currentCellNormals->GetTuple(cell[j], currentNormal);
280  vtkSmartPointer<vtkLine> line = vtkSmartPointer<vtkLine>::New();
281  line->GetPointIds()->SetNumberOfIds(2);
282  double newPoint[3];
283  double p0[3];
284  existingPoints->GetPoint(cell[j], p0);
285  newPoint[0] = p0[0] + currentNormal[0];
286  newPoint[1] = p0[1] + currentNormal[1];
287  newPoint[2] = p0[2] + currentNormal[2];
288 
289  line->GetPointIds()->SetId(0, idCounter);
290  newPoints->InsertPoint(idCounter, p0);
291  idCounter++;
292 
293  line->GetPointIds()->SetId(1, idCounter);
294  newPoints->InsertPoint(idCounter, newPoint);
295  idCounter++;
296 
297  newLines->InsertNextCell(line);
298  } // end for all points
299  } // end for all cells
300  } // end for all outputs
301 
302  newPolyData->SetPoints(newPoints);
303  newPolyData->SetLines(newLines);
304  newPolyData->BuildCells();
305 
307  surface->SetVtkPolyData(newPolyData);
308 
309  return surface;
310 }
311 
313 {
314  m_MaxSpacing = maxSpacing;
315 }
316 
318 {
319  Superclass::GenerateOutputInformation();
320 }
321 
323 {
324  for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++)
325  {
326  this->PopBackInput();
327  }
328  this->SetNumberOfIndexedInputs(0);
329  this->SetNumberOfIndexedOutputs(0);
330 
332  this->SetNthOutput(0, output.GetPointer());
333 }
334 
336 {
337  this->m_UseProgressBar = status;
338 }
339 
341 {
342  this->m_ProgressStepSize = stepSize;
343 }
void Progress(unsigned int steps=1)
Sets the current amount of progress to current progress + steps.
static char * line
Definition: svm.cpp:2870
float k(1.0)
Gives locked and index-based read access for a particular image part. The class provides several set-...
OutputType * GetOutput()
virtual const mitk::Surface * GetInput()
void SetProgressStepSize(unsigned int stepSize)
Set the stepsize which the progress bar should proceed.
static ProgressBar * GetInstance()
static method to get the GUI dependent ProgressBar-instance so the methods for steps to do and progre...
const TPixel & GetPixelByIndexSafe(const itk::Index< VDimension > &idx) const
void SetUseProgressBar(bool)
Set whether the mitkProgressBar should be used.
#define MITK_WARN
Definition: mitkLogMacros.h:19
An object of this class represents an exception of MITK. Please don&#39;t instantiate exceptions manually...
Definition: mitkException.h:45
void Normalize(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer im2, mitk::Image::Pointer mask1, std::string output)
Definition: CLBrainMask.cpp:36
static Pointer New()