Medical Imaging Interaction Toolkit  2018.4.99-4c24e3cb
Medical Imaging Interaction Toolkit
mitkPointCloudScoringFilter.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 <cmath>
16 
17 #include <vtkDoubleArray.h>
18 #include <vtkKdTree.h>
19 #include <vtkPointData.h>
20 #include <vtkPoints.h>
21 #include <vtkPolyVertex.h>
22 #include <vtkSmartPointer.h>
23 #include <vtkUnstructuredGrid.h>
24 
26 {
27  this->SetNumberOfIndexedOutputs(1);
28 }
29 
31 {
32 }
33 
35 {
36  if (UnstructuredGridToUnstructuredGridFilter::GetNumberOfInputs() != 2)
37  {
38  MITK_ERROR << "Not enough input arguments for PointCloudScoringFilter" << std::endl;
39  return;
40  }
41 
42  DataObjectPointerArray inputs = UnstructuredGridToUnstructuredGridFilter::GetInputs();
43  mitk::UnstructuredGrid::Pointer edgeGrid = dynamic_cast<mitk::UnstructuredGrid *>(inputs.at(0).GetPointer());
44  mitk::UnstructuredGrid::Pointer segmGrid = dynamic_cast<mitk::UnstructuredGrid *>(inputs.at(1).GetPointer());
45 
46  if (edgeGrid->IsEmpty() || segmGrid->IsEmpty())
47  {
48  if (edgeGrid->IsEmpty())
49  MITK_ERROR << "edgeGrid is empty" << std::endl;
50  if (segmGrid->IsEmpty())
51  MITK_ERROR << "segmGrid is empty" << std::endl;
52  }
53 
54  if (m_FilteredScores.size() > 0)
55  m_FilteredScores.clear();
56 
57  vtkSmartPointer<vtkUnstructuredGrid> edgevtkGrid = edgeGrid->GetVtkUnstructuredGrid();
58  vtkSmartPointer<vtkUnstructuredGrid> segmvtkGrid = segmGrid->GetVtkUnstructuredGrid();
59 
60  // KdTree from here
61  vtkSmartPointer<vtkPoints> kdPoints = vtkSmartPointer<vtkPoints>::New();
62  vtkSmartPointer<vtkKdTree> kdTree = vtkSmartPointer<vtkKdTree>::New();
63 
64  for (int i = 0; i < edgevtkGrid->GetNumberOfPoints(); i++)
65  {
66  kdPoints->InsertNextPoint(edgevtkGrid->GetPoint(i));
67  }
68 
69  kdTree->BuildLocatorFromPoints(kdPoints);
70  // KdTree until here
71 
72  vtkSmartPointer<vtkPoints> points = vtkSmartPointer<vtkPoints>::New();
73 
74  for (int i = 0; i < segmvtkGrid->GetNumberOfPoints(); i++)
75  {
76  points->InsertNextPoint(segmvtkGrid->GetPoint(i));
77  }
78 
79  std::vector<ScorePair> score;
80  std::vector<double> distances;
81 
82  double dist_glob = 0.0;
83  double dist = 0.0;
84 
85  for (int i = 0; i < points->GetNumberOfPoints(); i++)
86  {
87  double point[3];
88  points->GetPoint(i, point);
89  kdTree->FindClosestPoint(point[0], point[1], point[2], dist);
90  dist_glob += dist;
91  distances.push_back(dist);
92  score.push_back(std::make_pair(i, dist));
93  }
94 
95  double avg = dist_glob / points->GetNumberOfPoints();
96 
97  double tmpVar = 0.0;
98  double highest = 0.0;
99 
100  for (unsigned int i = 0; i < distances.size(); i++)
101  {
102  tmpVar = tmpVar + ((distances.at(i) - avg) * (distances.at(i) - avg));
103  if (distances.at(i) > highest)
104  highest = distances.at(i);
105  }
106 
107  // CUBIC MEAN
108  double cubicAll = 0.0;
109  for (unsigned i = 0; i < score.size(); i++)
110  {
111  cubicAll = cubicAll + score.at(i).second * score.at(i).second * score.at(i).second;
112  }
113  double root2 = cubicAll / static_cast<double>(score.size());
114  double cubic = pow(root2, 1.0 / 3.0);
115  // CUBIC END
116 
117  double metricValue = cubic;
118 
119  for (unsigned int i = 0; i < score.size(); i++)
120  {
121  if (score.at(i).second > metricValue)
122  {
123  m_FilteredScores.push_back(std::make_pair(score.at(i).first, score.at(i).second));
124  }
125  }
126 
127  m_NumberOfOutpPoints = m_FilteredScores.size();
128 
129  vtkSmartPointer<vtkPoints> filteredPoints = vtkSmartPointer<vtkPoints>::New();
130 
131  // storing the distances in the uGrid PointData
132  vtkSmartPointer<vtkDoubleArray> pointDataDistances = vtkSmartPointer<vtkDoubleArray>::New();
133  pointDataDistances->SetNumberOfComponents(1);
134  pointDataDistances->SetNumberOfTuples(m_FilteredScores.size());
135  pointDataDistances->SetName("Distances");
136 
137  for (unsigned int i = 0; i < m_FilteredScores.size(); i++)
138  {
139  mitk::Point3D point;
140  point = segmvtkGrid->GetPoint(m_FilteredScores.at(i).first);
141  filteredPoints->InsertNextPoint(point[0], point[1], point[2]);
142  if (score.at(i).second > 0.001)
143  {
144  double dist[1] = {score.at(i).second};
145  pointDataDistances->InsertTuple(i, dist);
146  }
147  else
148  {
149  double dist[1] = {0.0};
150  pointDataDistances->InsertTuple(i, dist);
151  }
152  }
153 
154  unsigned int numPoints = filteredPoints->GetNumberOfPoints();
155 
156  vtkSmartPointer<vtkPolyVertex> verts = vtkSmartPointer<vtkPolyVertex>::New();
157 
158  verts->GetPointIds()->SetNumberOfIds(numPoints);
159  for (unsigned int i = 0; i < numPoints; i++)
160  {
161  verts->GetPointIds()->SetId(i, i);
162  }
163 
164  vtkSmartPointer<vtkUnstructuredGrid> uGrid = vtkSmartPointer<vtkUnstructuredGrid>::New();
165  uGrid->Allocate(1);
166 
167  uGrid->InsertNextCell(verts->GetCellType(), verts->GetPointIds());
168  uGrid->SetPoints(filteredPoints);
169  uGrid->GetPointData()->AddArray(pointDataDistances);
170 
172  outputGrid->SetVtkUnstructuredGrid(uGrid);
173  this->SetNthOutput(0, outputGrid);
174 
175  score.clear();
176  distances.clear();
177 }
178 
180 {
181  Superclass::GenerateOutputInformation();
182 }
#define MITK_ERROR
Definition: mitkLogMacros.h:20
static Pointer New()
Class for storing unstructured grids (vtkUnstructuredGrid)