Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkReduceContourSetFilter.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 
18 
20 {
21  m_MaxSegmentLenght = 0;
22  m_StepSize = 10;
23  m_Tolerance = -1;
24  m_ReductionType = DOUGLAS_PEUCKER;
25  m_MaxSpacing = -1;
26  m_MinSpacing = -1;
27  this->m_UseProgressBar = false;
28  this->m_ProgressStepSize = 1;
29  m_NumberOfPointsAfterReduction = 0;
30 
32  this->SetNthOutput(0, output.GetPointer());
33 }
34 
36 {
37 }
38 
39 void mitk::ReduceContourSetFilter::SetInput(unsigned int idx, const mitk::Surface *surface)
40 {
41  this->SetNthInput(idx, const_cast<mitk::Surface *>(surface));
42  this->Modified();
43 }
44 
46 {
47  this->SetInput(0, const_cast<mitk::Surface *>(surface));
48 }
49 
51 {
52  unsigned int numberOfInputs = this->GetNumberOfIndexedInputs();
53  unsigned int numberOfOutputs(0);
54 
55  vtkSmartPointer<vtkPolyData> newPolyData;
56  vtkSmartPointer<vtkCellArray> newPolygons;
57  vtkSmartPointer<vtkPoints> newPoints;
58 
59  // For the purpose of evaluation
60  // unsigned int numberOfPointsBefore (0);
61  m_NumberOfPointsAfterReduction = 0;
62 
63  for (unsigned int i = 0; i < numberOfInputs; i++)
64  {
65  mitk::Surface *currentSurface = const_cast<mitk::Surface *>(this->GetInput(i));
66  vtkSmartPointer<vtkPolyData> polyData = currentSurface->GetVtkPolyData();
67 
68  newPolyData = vtkSmartPointer<vtkPolyData>::New();
69  newPolygons = vtkSmartPointer<vtkCellArray>::New();
70  newPoints = vtkSmartPointer<vtkPoints>::New();
71 
72  vtkSmartPointer<vtkCellArray> existingPolys = polyData->GetPolys();
73 
74  vtkSmartPointer<vtkPoints> existingPoints = polyData->GetPoints();
75 
76  existingPolys->InitTraversal();
77 
78  vtkIdType *cell(nullptr);
79  vtkIdType cellSize(0);
80 
81  for (existingPolys->InitTraversal(); existingPolys->GetNextCell(cellSize, cell);)
82  {
83  bool incorporatePolygon =
84  this->CheckForIntersection(cell, cellSize, existingPoints, /*numberOfIntersections, intersectionPoints, */ i);
85  if (!incorporatePolygon)
86  continue;
87 
88  vtkSmartPointer<vtkPolygon> newPolygon = vtkSmartPointer<vtkPolygon>::New();
89 
90  if (m_ReductionType == NTH_POINT)
91  {
92  this->ReduceNumberOfPointsByNthPoint(cellSize, cell, existingPoints, newPolygon, newPoints);
93  if (newPolygon->GetPointIds()->GetNumberOfIds() != 0)
94  {
95  newPolygons->InsertNextCell(newPolygon);
96  }
97  }
98  else if (m_ReductionType == DOUGLAS_PEUCKER)
99  {
100  this->ReduceNumberOfPointsByDouglasPeucker(cellSize, cell, existingPoints, newPolygon, newPoints);
101  if (newPolygon->GetPointIds()->GetNumberOfIds() > 3)
102  {
103  newPolygons->InsertNextCell(newPolygon);
104  }
105  }
106 
107  // Again for evaluation
108  // numberOfPointsBefore += cellSize;
109  m_NumberOfPointsAfterReduction += newPolygon->GetPointIds()->GetNumberOfIds();
110  }
111 
112  if (newPolygons->GetNumberOfCells() != 0)
113  {
114  newPolyData->SetPolys(newPolygons);
115  newPolyData->SetPoints(newPoints);
116  newPolyData->BuildLinks();
117 
118  this->SetNumberOfIndexedOutputs(numberOfOutputs + 1);
120  this->SetNthOutput(numberOfOutputs, surface.GetPointer());
121  surface->SetVtkPolyData(newPolyData);
122  numberOfOutputs++;
123  }
124  }
125 
126  // MITK_INFO<<"Points before: "<<numberOfPointsBefore<<" ##### Points after: "<<numberOfPointsAfter;
127  this->SetNumberOfIndexedOutputs(numberOfOutputs);
128 
129  if (numberOfOutputs == 0)
130  {
132  tmp_output->SetVtkPolyData(vtkPolyData::New());
133  this->SetNthOutput(0, tmp_output.GetPointer());
134  }
135  // Setting progressbar
136  if (this->m_UseProgressBar)
137  mitk::ProgressBar::GetInstance()->Progress(this->m_ProgressStepSize);
138 }
139 
140 void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByNthPoint(
141  vtkIdType cellSize, vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints)
142 {
143  unsigned int newNumberOfPoints(0);
144  unsigned int mod = cellSize % m_StepSize;
145 
146  if (mod == 0)
147  {
148  newNumberOfPoints = cellSize / m_StepSize;
149  }
150  else
151  {
152  newNumberOfPoints = ((cellSize - mod) / m_StepSize) + 1;
153  }
154 
155  if (newNumberOfPoints <= 3)
156  {
157  return;
158  }
159  reducedPolygon->GetPointIds()->SetNumberOfIds(newNumberOfPoints);
160  reducedPolygon->GetPoints()->SetNumberOfPoints(newNumberOfPoints);
161 
162  for (vtkIdType i = 0; i < cellSize; i++)
163  {
164  if (i % m_StepSize == 0)
165  {
166  double point[3];
167  points->GetPoint(cell[i], point);
168  vtkIdType id = reducedPoints->InsertNextPoint(point);
169  reducedPolygon->GetPointIds()->SetId(i / m_StepSize, id);
170  }
171  }
172  vtkIdType id = cell[0];
173  double point[3];
174  points->GetPoint(id, point);
175  id = reducedPoints->InsertNextPoint(point);
176  reducedPolygon->GetPointIds()->SetId(newNumberOfPoints - 1, id);
177 }
178 
179 void mitk::ReduceContourSetFilter::ReduceNumberOfPointsByDouglasPeucker(
180  vtkIdType cellSize, vtkIdType *cell, vtkPoints *points, vtkPolygon *reducedPolygon, vtkPoints *reducedPoints)
181 {
182  // If the cell is too small to obtain a reduced polygon with the given stepsize return
183  if (cellSize <= static_cast<vtkIdType>(m_StepSize * 3))
184  return;
185 
186  /*
187  What we do now is (see the Douglas Peucker Algorithm):
188 
189  1. Divide the current contour in two line segments (start - middle; middle - end), put them into the stack
190 
191  2. Fetch first line segment and create the following vectors:
192 
193  - v1 = (start;end)
194  - v2 = (start;currentPoint) -> for each point of the current line segment!
195 
196  3. Calculate the distance from the currentPoint to v1:
197 
198  a. Determine the length of the orthogonal projection of v2 to v1 by:
199 
200  l = v2 * (normalized v1)
201 
202  b. There a three possibilities for the distance then:
203 
204  d = sqrt(lenght(v2)^2 - l^2) if l > 0 and l < length(v1)
205  d = lenght(v2-v1) if l > 0 and l > lenght(v1)
206  d = length(v2) if l < 0 because v2 is then pointing in a different direction than v1
207 
208  4. Memorize the point with the biggest distance and create two new line segments with it at the end of the iteration
209  and put it into the stack
210 
211  5. If the distance value D <= m_Tolerance, then add the start and end index and the corresponding points to the
212  reduced ones
213  */
214 
215  // First of all set tolerance if none is specified
216  if (m_Tolerance < 0)
217  {
218  if (m_MaxSpacing > 0)
219  {
220  m_Tolerance = m_MinSpacing;
221  }
222  else
223  {
224  m_Tolerance = 1.5;
225  }
226  }
227 
228  std::stack<LineSegment> lineSegments;
229 
230  // 1. Divide in line segments
231 
232  LineSegment ls2;
233  ls2.StartIndex = cell[cellSize / 2];
234  ls2.EndIndex = cell[cellSize - 1];
235  lineSegments.push(ls2);
236 
237  LineSegment ls1;
238  ls1.StartIndex = cell[0];
239  ls1.EndIndex = cell[cellSize / 2];
240  lineSegments.push(ls1);
241 
242  LineSegment currentSegment;
243  double v1[3];
244  double v2[3];
245  double tempV[3];
246  double lenghtV1;
247 
248  double currentMaxDistance(0);
249  vtkIdType currentMaxDistanceIndex(0);
250 
251  double l;
252  double d;
253 
254  vtkIdType pointId(0);
255  // Add the start index to the reduced points. From now on just the end indices will be added
256  pointId = reducedPoints->InsertNextPoint(points->GetPoint(cell[0]));
257  reducedPolygon->GetPointIds()->InsertNextId(pointId);
258 
259  while (!lineSegments.empty())
260  {
261  currentSegment = lineSegments.top();
262  lineSegments.pop();
263 
264  // 2. Create vectors
265 
266  points->GetPoint(currentSegment.EndIndex, tempV);
267  points->GetPoint(currentSegment.StartIndex, v1);
268 
269  v1[0] = tempV[0] - v1[0];
270  v1[1] = tempV[1] - v1[1];
271  v1[2] = tempV[2] - v1[2];
272 
273  lenghtV1 = vtkMath::Norm(v1);
274 
275  vtkMath::Normalize(v1);
276  int range = currentSegment.EndIndex - currentSegment.StartIndex;
277  for (int i = 1; i < abs(range); ++i)
278  {
279  points->GetPoint(currentSegment.StartIndex + i, tempV);
280  points->GetPoint(currentSegment.StartIndex, v2);
281 
282  v2[0] = tempV[0] - v2[0];
283  v2[1] = tempV[1] - v2[1];
284  v2[2] = tempV[2] - v2[2];
285 
286  // 3. Calculate the distance
287 
288  l = vtkMath::Dot(v2, v1);
289 
290  d = vtkMath::Norm(v2);
291 
292  if (l > 0 && l < lenghtV1)
293  {
294  d = sqrt((d * d - l * l));
295  }
296  else if (l > 0 && l > lenghtV1)
297  {
298  tempV[0] = lenghtV1 * v1[0] - v2[0];
299  tempV[1] = lenghtV1 * v1[1] - v2[1];
300  tempV[2] = lenghtV1 * v1[2] - v2[2];
301  d = vtkMath::Norm(tempV);
302  }
303 
304  // 4. Memorize maximum distance
305  if (d > currentMaxDistance)
306  {
307  currentMaxDistance = d;
308  currentMaxDistanceIndex = currentSegment.StartIndex + i;
309  }
310  }
311 
312  // 4. & 5.
313  if (currentMaxDistance <= m_Tolerance)
314  {
315  // double temp[3];
316  int segmentLenght = currentSegment.EndIndex - currentSegment.StartIndex;
317  if (segmentLenght > (int)m_MaxSegmentLenght)
318  {
319  m_MaxSegmentLenght = (unsigned int)segmentLenght;
320  }
321 
322  // MITK_INFO<<"Lenght: "<<abs(segmentLenght);
323  if (abs(segmentLenght) > 25)
324  {
325  unsigned int newLenght(segmentLenght);
326  while (newLenght > 25)
327  {
328  newLenght = newLenght * 0.5;
329  }
330  unsigned int divisions = abs(segmentLenght) / newLenght;
331  // MITK_INFO<<"Divisions: "<<divisions;
332 
333  for (unsigned int i = 1; i <= divisions; ++i)
334  {
335  // MITK_INFO<<"Inserting MIDDLE: "<<(currentSegment.StartIndex + newLenght*i);
336  pointId = reducedPoints->InsertNextPoint(points->GetPoint(currentSegment.StartIndex + newLenght * i));
337  reducedPolygon->GetPointIds()->InsertNextId(pointId);
338  }
339  }
340  // MITK_INFO<<"Inserting END: "<<currentSegment.EndIndex;
341  pointId = reducedPoints->InsertNextPoint(points->GetPoint(currentSegment.EndIndex));
342  reducedPolygon->GetPointIds()->InsertNextId(pointId);
343  }
344  else
345  {
346  ls2.StartIndex = currentMaxDistanceIndex;
347  ls2.EndIndex = currentSegment.EndIndex;
348  lineSegments.push(ls2);
349 
350  ls1.StartIndex = currentSegment.StartIndex;
351  ls1.EndIndex = currentMaxDistanceIndex;
352  lineSegments.push(ls1);
353  }
354  currentMaxDistance = 0;
355  }
356 }
357 
358 bool mitk::ReduceContourSetFilter::CheckForIntersection(
359  vtkIdType *currentCell,
360  vtkIdType currentCellSize,
361  vtkPoints *currentPoints,
362  /* vtkIdType numberOfIntersections, vtkIdType* intersectionPoints,*/ unsigned int currentInputIndex)
363 {
364  /*
365  If we check the current cell for intersections then we have to consider three possibilies:
366  1. There is another cell among all the other input surfaces which intersects the current polygon:
367  - That means we have to save the intersection points because these points should not be eliminated
368  2. There current polygon exists just because of an intersection of another polygon with the current plane defined by
369  the current polygon
370  - That means the current polygon should not be incorporated and all of its points should be eliminated
371  3. There is no intersection
372  - That mean we can just reduce the current polygons points without considering any intersections
373  */
374 
375  for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++)
376  {
377  // Don't check for intersection with the polygon itself
378  if (i == currentInputIndex)
379  continue;
380 
381  // Get the next polydata to check for intersection
382  vtkSmartPointer<vtkPolyData> poly = const_cast<Surface *>(this->GetInput(i))->GetVtkPolyData();
383  vtkSmartPointer<vtkCellArray> polygonArray = poly->GetPolys();
384  polygonArray->InitTraversal();
385  vtkIdType anotherInputPolygonSize(0);
386  vtkIdType *anotherInputPolygonIDs(nullptr);
387 
388  /*
389  The procedure is:
390  - Create the equation of the plane, defined by the points of next input
391  - Calculate the distance of each point of the current polygon to the plane
392  - If the maximum distance is not bigger than 1.5 of the maximum spacing AND the minimal distance is not bigger
393  than 0.5 of the minimum spacing then the current contour is an intersection contour
394  */
395 
396  for (polygonArray->InitTraversal(); polygonArray->GetNextCell(anotherInputPolygonSize, anotherInputPolygonIDs);)
397  {
398  // Choosing three plane points to calculate the plane vectors
399  double p1[3];
400  double p2[3];
401  double p3[3];
402 
403  // The plane vectors
404  double v1[3];
405  double v2[3] = {0};
406  // The plane normal
407  double normal[3];
408 
409  // Create first Vector
410  poly->GetPoint(anotherInputPolygonIDs[0], p1);
411  poly->GetPoint(anotherInputPolygonIDs[1], p2);
412 
413  v1[0] = p2[0] - p1[0];
414  v1[1] = p2[1] - p1[1];
415  v1[2] = p2[2] - p1[2];
416 
417  // Find 3rd point for 2nd vector (The angle between the two plane vectors should be bigger than 30 degrees)
418 
419  double maxDistance(0);
420  double minDistance(10000);
421 
422  for (vtkIdType j = 2; j < anotherInputPolygonSize; j++)
423  {
424  poly->GetPoint(anotherInputPolygonIDs[j], p3);
425 
426  v2[0] = p3[0] - p1[0];
427  v2[1] = p3[1] - p1[1];
428  v2[2] = p3[2] - p1[2];
429 
430  // Calculate the angle between the two vector for the current point
431  double dotV1V2 = vtkMath::Dot(v1, v2);
432  double absV1 = sqrt(vtkMath::Dot(v1, v1));
433  double absV2 = sqrt(vtkMath::Dot(v2, v2));
434  double cosV1V2 = dotV1V2 / (absV1 * absV2);
435 
436  double arccos = acos(cosV1V2);
437  double degree = vtkMath::DegreesFromRadians(arccos);
438 
439  // If angle is bigger than 30 degrees break
440  if (degree > 30)
441  break;
442 
443  } // for (to find 3rd point)
444 
445  // Calculate normal of the plane by taking the cross product of the two vectors
446  vtkMath::Cross(v1, v2, normal);
447  vtkMath::Normalize(normal);
448 
449  // Determine position of the plane
450  double lambda = vtkMath::Dot(normal, p1);
451 
452  /*
453  Calculate the distance to the plane for each point of the current polygon
454  If the distance is zero then save the currentPoint as intersection point
455  */
456  for (vtkIdType k = 0; k < currentCellSize; k++)
457  {
458  double currentPoint[3];
459  currentPoints->GetPoint(currentCell[k], currentPoint);
460 
461  double tempPoint[3];
462  tempPoint[0] = normal[0] * currentPoint[0];
463  tempPoint[1] = normal[1] * currentPoint[1];
464  tempPoint[2] = normal[2] * currentPoint[2];
465 
466  double temp = tempPoint[0] + tempPoint[1] + tempPoint[2] - lambda;
467  double distance = fabs(temp);
468 
469  if (distance > maxDistance)
470  {
471  maxDistance = distance;
472  }
473  if (distance < minDistance)
474  {
475  minDistance = distance;
476  }
477  } // for (to calculate distance and intersections with currentPolygon)
478 
479  if (maxDistance < 1.5 * m_MaxSpacing && minDistance < 0.5 * m_MinSpacing)
480  {
481  return false;
482  }
483 
484  // Because we are considering the plane defined by the acual input polygon only one iteration is sufficient
485  // We do not need to consider each cell of the plane
486  break;
487  } // for (to traverse through all cells of actualInputPolyData)
488 
489  } // for (to iterate through all inputs)
490 
491  return true;
492 }
493 
495 {
496  Superclass::GenerateOutputInformation();
497 }
498 
500 {
501  for (unsigned int i = 0; i < this->GetNumberOfIndexedInputs(); i++)
502  {
503  this->PopBackInput();
504  }
505  this->SetNumberOfIndexedInputs(0);
506  this->SetNumberOfIndexedOutputs(0);
507 
508  // BUG XXXXX Fix
510  this->SetNthOutput(0, output.GetPointer());
511 
512  m_NumberOfPointsAfterReduction = 0;
513 }
514 
516 {
517  this->m_UseProgressBar = status;
518 }
519 
521 {
522  this->m_ProgressStepSize = stepSize;
523 }
void Progress(unsigned int steps=1)
Sets the current amount of progress to current progress + steps.
Class for storing surfaces (vtkPolyData).
Definition: mitkSurface.h:32
virtual vtkPolyData * GetVtkPolyData(unsigned int t=0) const
void SetUseProgressBar(bool)
Set whether the mitkProgressBar should be used.
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...
virtual void SetInput(const mitk::Surface *surface)
void Normalize(itk::Image< TPixel, VImageDimension > *itkImage, mitk::Image::Pointer im2, mitk::Image::Pointer mask1, std::string output)
Definition: CLBrainMask.cpp:40
virtual void GenerateOutputInformation() override
static Pointer New()
static itkEventMacro(BoundingShapeInteractionEvent, itk::AnyEvent) class MITKBOUNDINGSHAPE_EXPORT BoundingShapeInteractor Pointer New()
Basic interaction methods for mitk::GeometryData.