Medical Imaging Interaction Toolkit  2018.4.99-3e3f1a6e
Medical Imaging Interaction Toolkit
mitkPlaneFit.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 "mitkPlaneFit.h"
14 
15 #include "mitkGeometryData.h"
16 #include "mitkPlaneGeometry.h"
18 
19 #include <vcl_iostream.h>
20 #include <vnl/algo/vnl_svd.h>
21 
22 mitk::PlaneFit::PlaneFit() : m_PointSet(nullptr)
23 {
24  m_TimeGeometry = mitk::ProportionalTimeGeometry::New();
25 }
26 
28 {
29 }
30 
32 {
33  mitk::PointSet::ConstPointer input = this->GetInput();
34  mitk::GeometryData::Pointer output = this->GetOutput();
35 
36  itkDebugMacro(<< "GenerateOutputInformation()");
37 
38  if (input.IsNull())
39  return;
40 
41  if (m_PointSet == nullptr)
42  {
43  return;
44  }
45 
46  bool update = false;
47  if (output->GetGeometry() == nullptr || output->GetTimeGeometry() == nullptr)
48  update = true;
49  if ((!update) && (output->GetTimeGeometry()->CountTimeSteps() != input->GetTimeGeometry()->CountTimeSteps()))
50  update = true;
51  if (update)
52  {
54 
56  dynamic_cast<ProportionalTimeGeometry *>(m_TimeGeometry.GetPointer());
57  timeGeometry->Initialize(planeGeometry, m_PointSet->GetPointSetSeriesSize());
58  // m_TimeGeometry->InitializeEvenlyTimed(
59  // planeGeometry, m_PointSet->GetPointSetSeriesSize() );
60 
61  TimeStepType timeStep;
62  for (timeStep = 0; (timeStep < m_PointSet->GetPointSetSeriesSize()) && (timeStep < m_Planes.size()); ++timeStep)
63  {
64  timeGeometry->SetTimeStepGeometry(m_Planes[timeStep], timeStep);
65  }
66 
67  output->SetTimeGeometry(m_TimeGeometry);
68  }
69 }
70 
72 {
73  unsigned int t;
74  for (t = 0; t < m_PointSet->GetPointSetSeriesSize(); ++t)
75  {
76  // check number of data points - less then 3points isn't enough
77  if (m_PointSet->GetSize(t) >= 3)
78  {
79  this->CalculateCentroid(t);
80 
81  this->ProcessPointSet(t);
82 
83  this->InitializePlane(t);
84  }
85  }
86 }
87 
89 {
90  // Process object is not const-correct so the const_cast is required here
91  this->ProcessObject::SetNthInput(0, const_cast<mitk::PointSet *>(pointSet));
92 
93  m_PointSet = pointSet;
94  unsigned int pointSetSize = pointSet->GetPointSetSeriesSize();
95 
96  m_Planes.resize(pointSetSize);
97  m_Centroids.resize(pointSetSize);
98  m_PlaneVectors.resize(pointSetSize);
99 
100  unsigned int t;
101  for (t = 0; t < pointSetSize; ++t)
102  {
103  m_Planes[t] = mitk::PlaneGeometry::New();
104  }
105 }
106 
108 {
109  if (this->GetNumberOfInputs() < 1)
110  {
111  return nullptr;
112  }
113 
114  return static_cast<const mitk::PointSet *>(this->ProcessObject::GetInput(0));
115 }
116 
118 {
119  if (m_PointSet == nullptr)
120  return;
121 
122  int ps_total = m_PointSet->GetSize(t);
123 
124  m_Centroids[t][0] = m_Centroids[t][1] = m_Centroids[t][2] = 0.0;
125 
126  for (int i = 0; i < ps_total; i++)
127  {
128  mitk::Point3D p3d = m_PointSet->GetPoint(i, t);
129  m_Centroids[t][0] += p3d[0];
130  m_Centroids[t][1] += p3d[1];
131  m_Centroids[t][2] += p3d[2];
132  }
133 
134  // calculation of centroid
135  m_Centroids[t][0] /= ps_total;
136  m_Centroids[t][1] /= ps_total;
137  m_Centroids[t][2] /= ps_total;
138 }
139 
141 {
142  if (m_PointSet == nullptr)
143  return;
144 
145  // int matrix with POINTS x (X,Y,Z)
146  vnl_matrix<mitk::ScalarType> dataM(m_PointSet->GetSize(t), 3);
147 
148  int ps_total = m_PointSet->GetSize(t);
149  for (int i = 0; i < ps_total; i++)
150  {
151  mitk::Point3D p3d = m_PointSet->GetPoint(i, t);
152  dataM[i][0] = p3d[0] - m_Centroids[t][0];
153  dataM[i][1] = p3d[1] - m_Centroids[t][1];
154  dataM[i][2] = p3d[2] - m_Centroids[t][2];
155  }
156  // process the SVD (singular value decomposition) from ITK
157  // the vector will be orderd descending
158  vnl_svd<mitk::ScalarType> svd(dataM, 0.0);
159 
160  // calculate the SVD of A
161  vnl_vector<mitk::ScalarType> v = svd.nullvector();
162 
163  // Avoid erratic normal sign switching when the plane changes minimally
164  // by negating the vector for negative x values.
165  if (v[0] < 0)
166  {
167  v = -v;
168  }
169 
170  m_PlaneVectors[t][0] = v[0];
171  m_PlaneVectors[t][1] = v[1];
172  m_PlaneVectors[t][2] = v[2];
173 }
174 
176 {
177  return m_Planes[t];
178 }
179 
181 {
182  return m_PlaneVectors[t];
183 }
184 
186 {
187  return m_Centroids[t];
188 }
189 
191 {
192  m_Planes[t]->InitializePlane(m_Centroids[t], m_PlaneVectors[t]);
193 }
void CalculateCentroid(int t=0)
void InitializePlane(int t=0)
virtual unsigned int GetPointSetSeriesSize() const
virtual void SetInput(const mitk::PointSet *ps)
void GenerateOutputInformation() override
void ProcessPointSet(int t=0)
OutputType * GetOutput()
void GenerateData() override
virtual int GetSize(unsigned int t=0) const
returns the current size of the point-list
Data structure which stores a set of points. Superclass of mitk::Mesh.
Definition: mitkPointSet.h:75
std::vcl_size_t TimeStepType
virtual mitk::PlaneGeometry::Pointer GetPlaneGeometry(int t=0)
~PlaneFit() override
static Pointer New()
virtual const mitk::Vector3D & GetPlaneNormal(int t=0) const
virtual const mitk::Point3D & GetCentroid(int t=0) const
PointType GetPoint(PointIdentifier id, int t=0) const
Get the point with ID id in world coordinates.
const mitk::PointSet * GetInput()