ITK  5.2.0
Insight Toolkit
itkQuadEdgeMeshDecimationQuadricElementHelper.h
Go to the documentation of this file.
1 /*=========================================================================
2  *
3  * Copyright NumFOCUS
4  *
5  * Licensed under the Apache License, Version 2.0 (the "License");
6  * you may not use this file except in compliance with the License.
7  * You may obtain a copy of the License at
8  *
9  * http://www.apache.org/licenses/LICENSE-2.0.txt
10  *
11  * Unless required by applicable law or agreed to in writing, software
12  * distributed under the License is distributed on an "AS IS" BASIS,
13  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
14  * See the License for the specific language governing permissions and
15  * limitations under the License.
16  *
17  *=========================================================================*/
18 #ifndef itkQuadEdgeMeshDecimationQuadricElementHelper_h
19 #define itkQuadEdgeMeshDecimationQuadricElementHelper_h
20 
21 #include "itkPoint.h"
22 #include "vnl/vnl_vector_fixed.h"
23 #include "vnl/vnl_matrix.h"
24 #include "vnl/algo/vnl_matrix_inverse.h"
25 
26 #include "itkTriangleHelper.h"
27 
28 namespace itk
29 {
31 template <typename TPoint>
33 {
34 public:
36  using PointType = TPoint;
37  using CoordType = typename PointType::CoordRepType;
38 
39  static constexpr unsigned int PointDimension = PointType::PointDimension;
40  static constexpr unsigned int NumberOfCoefficients = PointDimension * (PointDimension + 1 / 2 + PointDimension + 1);
41 
43  using VNLMatrixType = vnl_matrix<CoordType>;
44  using VNLVectorType = vnl_vector_fixed<CoordType, Self::PointDimension>;
45  using CoefficientVectorType = vnl_vector_fixed<CoordType, Self::NumberOfCoefficients>;
47 
48  // *****************************************************************
50  : m_Coefficients(itk::NumericTraits<CoordType>::ZeroValue())
52  , m_B(itk::NumericTraits<CoordType>::ZeroValue())
53  , m_SVDAbsoluteThreshold(static_cast<CoordType>(1e-6))
54  , m_SVDRelativeThreshold(static_cast<CoordType>(1e-3))
55  {
56  this->m_Rank = PointDimension;
57  }
58 
60  : m_Coefficients(iCoefficients)
62  , m_B(itk::NumericTraits<CoordType>::ZeroValue())
63  , m_SVDAbsoluteThreshold(static_cast<CoordType>(1e-3))
64  , m_SVDRelativeThreshold(static_cast<CoordType>(1e-3))
65  {
66  this->m_Rank = PointDimension;
68  }
69 
71 
74  {
75  return this->m_Coefficients;
76  }
77 
80  {
82  return m_A;
83  }
84 
87  {
89  return m_B;
90  }
91 
92  unsigned int
93  GetRank() const
94  {
95  return m_Rank;
96  }
97 
99  inline CoordType
100  ComputeError(const PointType & iP) const
101  {
102  // ComputeAMatrixAndBVector();
103  vnl_svd<CoordType> svd(m_A, m_SVDAbsoluteThreshold);
104  svd.zero_out_relative(m_SVDRelativeThreshold);
105  CoordType oError = inner_product(iP.GetVnlVector(), svd.recompose() * iP.GetVnlVector());
106 
107  return this->m_Coefficients.back() - oError;
108  /*
109  CoordType oError( 0. );
110 
111  std::vector< CoordType > pt( PointDimension + 1, 1. );
112 
113  unsigned int dim1( 0 ), dim2, k( 0 );
114 
115  while( dim1 < PointDimension )
116  {
117  pt[dim1] = iP[dim1];
118  ++dim1;
119  }
120 
121  for( dim1 = 0; dim1 < PointDimension + 1; ++dim1 )
122  {
123  oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1];
124 
125  for( dim2 = dim1 + 1; dim2 < PointDimension + 1; ++dim2 )
126  {
127  oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2];
128  }
129  }
130  oError += this->m_Coefficients[k++];
131 
132  return oError;*/
133  }
134 
136  inline CoordType
138  {
139  PointType optimal_location = ComputeOptimalLocation(iP);
140 
141  return ComputeError(optimal_location);
142  }
143 
144  PointType
146  {
148 
149  vnl_svd<CoordType> svd(m_A, m_SVDAbsoluteThreshold);
150  svd.zero_out_relative(m_SVDRelativeThreshold);
151 
152  m_Rank = svd.rank();
153 
154  const auto y = (m_B.as_vector() - m_A * iP.GetVnlVector());
155 
156  VNLVectorType displacement = svd.solve(y);
157 
158  PointType oP;
159  for (unsigned int dim = 0; dim < PointDimension; dim++)
160  {
161  oP[dim] = iP[dim] + displacement[dim];
162  }
163  return oP;
164  }
165 
167  PointType
168  ComputeOptimalLocation(const unsigned int &)
169  {}
170 
171  void
172  AddTriangle(const PointType & iP1,
173  const PointType & iP2,
174  const PointType & iP3,
175  const CoordType & iWeight = static_cast<CoordType>(1.))
176  {
177  VectorType N = TriangleType::ComputeNormal(iP1, iP2, iP3);
178 
179  AddPoint(iP1, N, iWeight);
180  }
181 
182  void
183  AddPoint(const PointType & iP, const VectorType & iN, const CoordType & iWeight = static_cast<CoordType>(1.))
184  {
185  unsigned int k(0), dim1, dim2;
186 
187  CoordType d = -iN * iP.GetVectorFromOrigin();
188 
189  for (dim1 = 0; dim1 < PointDimension; ++dim1)
190  {
191  for (dim2 = dim1; dim2 < PointDimension; ++dim2)
192  {
193  this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
194  }
195  this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
196  }
197 
198  this->m_Coefficients[k++] += iWeight * d * d;
199  }
200 
201  // ***********************************************************************
202  // operators
203  Self &
204  operator=(const Self & iRight)
205  {
206  if (this != &iRight)
207  {
208  this->m_Coefficients = iRight.m_Coefficients;
209  }
210  return *this;
211  }
212 
213  Self
214  operator+(const Self & iRight) const
215  {
216  return Self(this->m_Coefficients + iRight.m_Coefficients);
217  }
218 
219  Self &
220  operator+=(const Self & iRight)
221  {
222  this->m_Coefficients += iRight.m_Coefficients;
223  return *this;
224  }
225 
226  Self
227  operator-(const Self & iRight) const
228  {
229  return Self(this->m_Coefficients - iRight.m_Coefficients);
230  }
231 
232  Self &
233  operator-=(const Self & iRight)
234  {
235  this->m_Coefficients -= iRight.m_Coefficients;
236  return *this;
237  }
238 
239  Self operator*(const CoordType & iV) const
240  {
241  Self oElement = Self(this->m_Coefficients * iV);
242 
243  return oElement;
244  }
245 
246  Self &
247  operator*=(const CoordType & iV)
248  {
249  this->m_Coefficients *= iV;
250  return *this;
251  }
252 
253 protected:
257  unsigned int m_Rank;
260  // bool m_MatrixFilled;
261 
262  void
264  {
265  unsigned int k(0), dim1, dim2;
266 
267  for (dim1 = 0; dim1 < PointDimension; ++dim1)
268  {
269  for (dim2 = dim1; dim2 < PointDimension; ++dim2)
270  {
271  m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++];
272  }
273  m_B[dim1] = -m_Coefficients[k++];
274  }
275  // m_MatrixFilled = true;
276  }
277 };
278 } // namespace itk
279 #endif
itk::QuadEdgeMeshDecimationQuadricElementHelper::operator*=
Self & operator*=(const CoordType &iV)
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:247
itk::QuadEdgeMeshDecimationQuadricElementHelper::m_Rank
unsigned int m_Rank
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:257
itk::QuadEdgeMeshDecimationQuadricElementHelper::m_SVDRelativeThreshold
CoordType m_SVDRelativeThreshold
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:259
itk::QuadEdgeMeshDecimationQuadricElementHelper::m_A
VNLMatrixType m_A
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:255
itk::QuadEdgeMeshDecimationQuadricElementHelper::GetCoefficients
CoefficientVectorType GetCoefficients() const
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:73
itk::QuadEdgeMeshDecimationQuadricElementHelper::ComputeOptimalLocation
PointType ComputeOptimalLocation(const PointType &iP)
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:145
itk::QuadEdgeMeshDecimationQuadricElementHelper::m_SVDAbsoluteThreshold
CoordType m_SVDAbsoluteThreshold
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:258
itk::QuadEdgeMeshDecimationQuadricElementHelper::m_B
VNLVectorType m_B
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:256
itkPoint.h
itk::QuadEdgeMeshDecimationQuadricElementHelper::operator-
Self operator-(const Self &iRight) const
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:227
itk::GTest::TypedefsAndConstructors::Dimension2::VectorType
ImageBaseType::SpacingType VectorType
Definition: itkGTestTypedefsAndConstructors.h:53
itk::TriangleHelper
A convenience class for computation of various triangle elements in 2D or 3D.
Definition: itkTriangleHelper.h:31
itk::QuadEdgeMeshDecimationQuadricElementHelper::operator+=
Self & operator+=(const Self &iRight)
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:220
itk::QuadEdgeMeshDecimationQuadricElementHelper::VectorType
typename PointType::VectorType VectorType
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:42
itk::QuadEdgeMeshDecimationQuadricElementHelper::PointType
TPoint PointType
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:36
itk::QuadEdgeMeshDecimationQuadricElementHelper::ComputeErrorAtOptimalLocation
CoordType ComputeErrorAtOptimalLocation(const PointType &iP)
TODO this method should be really optimized!!!
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:137
itk::QuadEdgeMeshDecimationQuadricElementHelper::m_Coefficients
CoefficientVectorType m_Coefficients
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:254
itk::QuadEdgeMeshDecimationQuadricElementHelper::GetAMatrix
VNLMatrixType GetAMatrix()
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:79
itk::QuadEdgeMeshDecimationQuadricElementHelper::PointDimension
static constexpr unsigned int PointDimension
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:39
itk::QuadEdgeMeshDecimationQuadricElementHelper::ComputeOptimalLocation
PointType ComputeOptimalLocation(const unsigned int &)
TODO to be implemented!!!
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:168
itk::QuadEdgeMeshDecimationQuadricElementHelper::CoefficientVectorType
vnl_vector_fixed< CoordType, Self::NumberOfCoefficients > CoefficientVectorType
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:45
itk::QuadEdgeMeshDecimationQuadricElementHelper::AddTriangle
void AddTriangle(const PointType &iP1, const PointType &iP2, const PointType &iP3, const CoordType &iWeight=static_cast< CoordType >(1.))
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:172
itk::QuadEdgeMeshDecimationQuadricElementHelper::operator+
Self operator+(const Self &iRight) const
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:214
itk::QuadEdgeMeshDecimationQuadricElementHelper::operator-=
Self & operator-=(const Self &iRight)
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:233
itk::QuadEdgeMeshDecimationQuadricElementHelper::ComputeError
CoordType ComputeError(const PointType &iP) const
TODO this method should be really optimized!!!
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:100
itkTriangleHelper.h
itk::QuadEdgeMeshDecimationQuadricElementHelper::QuadEdgeMeshDecimationQuadricElementHelper
QuadEdgeMeshDecimationQuadricElementHelper()
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:49
itk::QuadEdgeMeshDecimationQuadricElementHelper::operator*
Self operator*(const CoordType &iV) const
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:239
itk::QuadEdgeMeshDecimationQuadricElementHelper::Self
QuadEdgeMeshDecimationQuadricElementHelper Self
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:35
itk::NumericTraits
Define additional traits for native types such as int or float.
Definition: itkNumericTraits.h:58
itk::QuadEdgeMeshDecimationQuadricElementHelper::AddPoint
void AddPoint(const PointType &iP, const VectorType &iN, const CoordType &iWeight=static_cast< CoordType >(1.))
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:183
itk::QuadEdgeMeshDecimationQuadricElementHelper::GetBVector
VNLVectorType GetBVector()
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:86
itk::QuadEdgeMeshDecimationQuadricElementHelper::QuadEdgeMeshDecimationQuadricElementHelper
QuadEdgeMeshDecimationQuadricElementHelper(const CoefficientVectorType &iCoefficients)
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:59
itk::TriangleHelper::ComputeNormal
static VectorType ComputeNormal(const PointType &iA, const PointType &iB, const PointType &iC)
Compute Normal vector to the triangle formed by (iA,iB,iC)
itk::QuadEdgeMeshDecimationQuadricElementHelper::~QuadEdgeMeshDecimationQuadricElementHelper
~QuadEdgeMeshDecimationQuadricElementHelper()=default
itk
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....
Definition: itkAnnulusOperator.h:24
itk::QuadEdgeMeshDecimationQuadricElementHelper::VNLVectorType
vnl_vector_fixed< CoordType, Self::PointDimension > VNLVectorType
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:44
itk::QuadEdgeMeshDecimationQuadricElementHelper::ComputeAMatrixAndBVector
void ComputeAMatrixAndBVector()
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:263
itk::QuadEdgeMeshDecimationQuadricElementHelper
TODO explicit specification for VDimension=3!!!
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:32
itk::QuadEdgeMeshDecimationQuadricElementHelper::NumberOfCoefficients
static constexpr unsigned int NumberOfCoefficients
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:40
itk::QuadEdgeMeshDecimationQuadricElementHelper::operator=
Self & operator=(const Self &iRight)
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:204
itk::QuadEdgeMeshDecimationQuadricElementHelper::GetRank
unsigned int GetRank() const
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:93
itk::Math::e
static constexpr double e
Definition: itkMath.h:54
itk::QuadEdgeMeshDecimationQuadricElementHelper::VNLMatrixType
vnl_matrix< CoordType > VNLMatrixType
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:43
itk::QuadEdgeMeshDecimationQuadricElementHelper::CoordType
typename PointType::CoordRepType CoordType
Definition: itkQuadEdgeMeshDecimationQuadricElementHelper.h:37