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