Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
mitkPlaneProposer.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 "mitkPlaneProposer.h"
18 
19 #include <mitkGeometryData.h>
20 #include <mitkPlaneFit.h>
21 #include <mitkPointSet.h>
23 #include <mitkUnstructuredGrid.h>
24 #include <mitkVector.h>
25 
26 #include <vtkCenterOfMass.h>
27 #include <vtkDoubleArray.h>
28 #include <vtkMath.h>
29 #include <vtkPointData.h>
30 #include <vtkSmartPointer.h>
31 #include <vtkUnstructuredGrid.h>
32 
33 bool compare(std::pair<double, int> i, std::pair<double, int> j)
34 {
35  return (i.first > j.first);
36 }
37 
38 mitk::PlaneProposer::PlaneProposer() : m_UseDistances(false), m_UseLeastSquares(false), m_NumberOfClustersToUse(3)
39 {
40 }
41 
43 {
44 }
45 
46 void mitk::PlaneProposer::SetUnstructuredGrids(std::vector<mitk::UnstructuredGrid::Pointer> &grids)
47 {
48  if (grids.empty())
49  {
50  MITK_WARN << "Input grids are empty!";
51  return;
52  }
53  m_Grids = grids;
54 }
55 
57 {
58  m_UseDistances = status;
59 }
60 
62 {
63  m_NumberOfClustersToUse = num;
64 }
65 
67 {
68  m_SNC = snc;
69 }
70 
71 std::array<std::array<double, 3>, 3> mitk::PlaneProposer::GetCentroids()
72 {
73  return m_Centroids;
74 }
75 
77 {
78  return m_ProposedPlaneInfo;
79 }
80 
82 {
83  // find the avg distances of every cluster maybe a square mean for better looking at high values?
84  std::vector<std::pair<double, int>> avgDistances;
85 
86  // get the number of points and the id of every cluster
87  std::vector<std::pair<int, int>> sizeIDs;
88 
89  for (unsigned int i = 0; i < m_Grids.size(); i++)
90  {
91  mitk::UnstructuredGrid::Pointer cluster = m_Grids.at(i);
92  vtkSmartPointer<vtkDoubleArray> data =
93  dynamic_cast<vtkDoubleArray *>(cluster->GetVtkUnstructuredGrid()->GetPointData()->GetArray(0));
94  double avg = 0.0;
95  for (int j = 0; j < data->GetSize(); j++)
96  {
97  avg += data->GetValue(j);
98  }
99  avgDistances.push_back(std::make_pair(avg / static_cast<double>(data->GetSize()), i));
100  }
101 
102  // now sort the clusters with their distances
103  for (unsigned int i = 0; i < m_Grids.size(); i++)
104  {
105  sizeIDs.push_back(std::make_pair(m_Grids.at(i)->GetVtkUnstructuredGrid()->GetNumberOfPoints(), i));
106  }
107 
108  // sort the point IDs for finding the biggest clusters and clusters with the
109  // highest distance
110  // sizeIDs is sorted by number of points in each cluster
111  // avgDistances is sorted by avg distance to next edge
112  std::sort(sizeIDs.begin(), sizeIDs.end(), compare);
113  std::sort(avgDistances.begin(), avgDistances.end(), compare);
114 
115  if (m_Grids.size() < m_NumberOfClustersToUse || m_UseLeastSquares)
116  {
117  // We need at least three points to define a plane
118  m_NumberOfClustersToUse = m_Grids.size();
119  m_ProposedPlaneInfo = this->CreatePlaneByLeastSquares(sizeIDs, avgDistances);
120  }
121  else
122  {
123  m_ProposedPlaneInfo = this->CreatePlaneByCentroids(sizeIDs, avgDistances);
124  }
125 
126  if (m_SNC.IsNotNull())
127  {
128  m_SNC->ReorientSlices(m_SNC->GetCurrentPlaneGeometry()->GetCenter(), m_ProposedPlaneInfo.x, m_ProposedPlaneInfo.y);
129  m_SNC->SelectSliceByPoint(m_ProposedPlaneInfo.pointOnPlane);
130  }
131 }
132 
133 mitk::PlaneProposer::PlaneInfo mitk::PlaneProposer::CreatePlaneByCentroids(
134  const std::vector<std::pair<int, int>> &sizeIDs, const std::vector<std::pair<double, int>> &avgDistances)
135 {
136  for (unsigned int j = 0; j < m_NumberOfClustersToUse; j++) // iterate over the clusters
137  {
138  // Get the cluster with "id" from all clusters, "id" is saved in sizeIDs.second
139  // sizeIDs.first represent the number of points which is only needed for sorting the "id"s
140  vtkSmartPointer<vtkUnstructuredGrid> tmpGrid;
141 
142  if (m_UseDistances)
143  tmpGrid = m_Grids.at(avgDistances.at(j).second)->GetVtkUnstructuredGrid(); // highest distance
144  else
145  tmpGrid = m_Grids.at(sizeIDs.at(j).second)->GetVtkUnstructuredGrid(); // biggest cluster
146 
147  double tmpCenter[3];
148  vtkCenterOfMass::ComputeCenterOfMass(tmpGrid->GetPoints(), NULL, tmpCenter);
149  std::array<double, 3> center;
150  center[0] = tmpCenter[0];
151  center[1] = tmpCenter[1];
152  center[2] = tmpCenter[2];
153  m_Centroids[j] = center;
154  }
155 
156  double x[3];
157  x[0] = m_Centroids[1][0] - m_Centroids[0][0];
158  x[1] = m_Centroids[1][1] - m_Centroids[0][1];
159  x[2] = m_Centroids[1][2] - m_Centroids[0][2];
160  double y[3];
161  y[0] = m_Centroids[2][0] - m_Centroids[0][0];
162  y[1] = m_Centroids[2][1] - m_Centroids[0][1];
163  y[2] = m_Centroids[2][2] - m_Centroids[0][2];
164  double normal[3];
165  vtkMath::Cross(x, y, normal); // TODO do it manually!
167  planeInfo.normal = normal;
168  planeInfo.x = x;
169  planeInfo.y = y;
170  planeInfo.pointOnPlane[0] = m_Centroids[0][0];
171  planeInfo.pointOnPlane[1] = m_Centroids[0][1];
172  planeInfo.pointOnPlane[2] = m_Centroids[0][2];
173  return planeInfo;
174 }
175 
176 mitk::PlaneProposer::PlaneInfo mitk::PlaneProposer::CreatePlaneByLeastSquares(
177  const std::vector<std::pair<int, int>> &sizeIDs, const std::vector<std::pair<double, int>> &avgDistances)
178 {
179  // Generate a pointset from UnstructuredGrid for the PlaneFitFilter:
181 
182  int pointId = 0;
183 
184  for (unsigned int j = 0; j < m_NumberOfClustersToUse; j++) // iterate over the clusters
185  {
186  // Get the cluster with "id" from all clusters, "id" is saved in sizeIDs.second
187  // sizeIDs.first represent the number of points which is only needed for sorting the "id"s
188  vtkSmartPointer<vtkUnstructuredGrid> tmpGrid;
189 
190  if (m_UseDistances)
191  tmpGrid = m_Grids.at(avgDistances.at(j).second)->GetVtkUnstructuredGrid(); // highest distance
192  else
193  tmpGrid = m_Grids.at(sizeIDs.at(j).second)->GetVtkUnstructuredGrid(); // biggest cluster
194 
195  for (int i = 0; i < tmpGrid->GetNumberOfPoints(); i++)
196  {
197  mitk::Point3D point;
198  point[0] = tmpGrid->GetPoint(i)[0];
199  point[1] = tmpGrid->GetPoint(i)[1];
200  point[2] = tmpGrid->GetPoint(i)[2];
201 
202  pointset->InsertPoint(pointId, point);
203  pointId++;
204  }
205  }
207  planeFilter->SetInput(pointset);
208  planeFilter->Update();
209 
210  mitk::GeometryData::Pointer geoData = planeFilter->GetOutput();
211 
212  if (geoData.IsNull())
213  {
214  mitkThrow() << "GeometryData output from PlaneFit filter is null!";
215  }
217  mitk::PlaneGeometry::Pointer plane = dynamic_cast<mitk::PlaneGeometry *>(geoData->GetGeometry());
218  planeInfo.normal = plane->GetNormal();
219  planeInfo.pointOnPlane = plane->GetCenter();
220  return planeInfo;
221 }
static Pointer New()
void SetUseDistances(bool)
If true, the three clusters with the biggest mean distances are used for plane proposal Required the ...
#define MITK_WARN
Definition: mitkLogMacros.h:23
bool compare(std::pair< double, int > i, std::pair< double, int > j)
Encapsulates the geometrical information needed to descripe a PlaneInfo.
void CreatePlaneInfo()
Creates the actual plane proposal.
#define mitkThrow()
void SetNumberOfClustersToUse(unsigned int)
Sets the number of the clusters to be used for plane creation (default=3)
void SetSliceNavigationController(itk::SmartPointer< mitk::SliceNavigationController > &snc)
PlaneInfo GetProposedPlaneInfo()
static Pointer New()
Describes a two-dimensional, rectangular plane.
std::array< std::array< double, 3 >, 3 > GetCentroids()
void SetUnstructuredGrids(std::vector< itk::SmartPointer< mitk::UnstructuredGrid >> &grids)