Medical Imaging Interaction Toolkit  2016.11.0
Medical Imaging Interaction Toolkit
All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
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()