ITK 6.0.0
Insight Toolkit
 
Loading...
Searching...
No Matches
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 * https://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
28namespace itk
29{
31template <typename TPoint>
33{
34public:
36 using PointType = TPoint;
37 using CoordType = typename PointType::CoordinateType;
38
39 static constexpr unsigned int PointDimension = PointType::PointDimension;
40 static constexpr unsigned int NumberOfCoefficients = PointDimension * (PointDimension + 1 / 2 + PointDimension + 1);
41
42 using VectorType = typename PointType::VectorType;
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 // *****************************************************************
57
68
70
71 [[nodiscard]] CoefficientVectorType
73 {
74 return this->m_Coefficients;
75 }
76
79 {
81 return m_A;
82 }
83
86 {
88 return m_B;
89 }
90
91 [[nodiscard]] unsigned int
92 GetRank() const
93 {
94 return m_Rank;
95 }
96
98 [[nodiscard]] inline CoordType
99 ComputeError(const PointType & iP) const
100 {
101 // ComputeAMatrixAndBVector();
102 vnl_svd<CoordType> svd(m_A, m_SVDAbsoluteThreshold);
103 svd.zero_out_relative(m_SVDRelativeThreshold);
104 const CoordType oError = inner_product(iP.GetVnlVector(), svd.recompose() * iP.GetVnlVector());
105
106 return this->m_Coefficients.back() - oError;
107 /*
108 CoordType oError( 0. );
109
110 std::vector< CoordType > pt( PointDimension + 1, 1. );
111
112 unsigned int dim1( 0 ), dim2, k( 0 );
113
114 while( dim1 < PointDimension )
115 {
116 pt[dim1] = iP[dim1];
117 ++dim1;
118 }
119
120 for( dim1 = 0; dim1 < PointDimension + 1; ++dim1 )
121 {
122 oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1];
123
124 for( dim2 = dim1 + 1; dim2 < PointDimension + 1; ++dim2 )
125 {
126 oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2];
127 }
128 }
129 oError += this->m_Coefficients[k++];
130
131 return oError;*/
132 }
133
135 inline CoordType
137 {
138 PointType optimal_location = ComputeOptimalLocation(iP);
139
140 return ComputeError(optimal_location);
141 }
142
145 {
147
148 vnl_svd<CoordType> svd(m_A, m_SVDAbsoluteThreshold);
149 svd.zero_out_relative(m_SVDRelativeThreshold);
150
151 m_Rank = svd.rank();
152
153 const auto y = (m_B.as_vector() - m_A * iP.GetVnlVector());
154
155 VNLVectorType displacement = svd.solve(y);
156
157 PointType oP;
158 for (unsigned int dim = 0; dim < PointDimension; ++dim)
159 {
160 oP[dim] = iP[dim] + displacement[dim];
161 }
162 return oP;
163 }
164
167 ComputeOptimalLocation(const unsigned int)
168 {}
169
170 void
172 const PointType & iP2,
173 const PointType & iP3,
174 const CoordType & iWeight = static_cast<CoordType>(1.))
175 {
176 const VectorType N = TriangleType::ComputeNormal(iP1, iP2, iP3);
177
178 AddPoint(iP1, N, iWeight);
179 }
180
181 void
182 AddPoint(const PointType & iP, const VectorType & iN, const CoordType & iWeight = static_cast<CoordType>(1.))
183 {
184 unsigned int k(0); /*one-line-declaration*/
185 const CoordType d = -iN * iP.GetVectorFromOrigin();
186 for (unsigned int dim1 = 0; dim1 < PointDimension; ++dim1)
187 {
188 for (unsigned int dim2 = dim1; dim2 < PointDimension; ++dim2)
189 {
190 this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
191 }
192 this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
193 }
194
195 this->m_Coefficients[k++] += iWeight * d * d;
196 }
197
198 // ***********************************************************************
199 // operators
200 Self &
201 operator=(const Self & iRight)
202 {
203 if (this != &iRight)
204 {
205 this->m_Coefficients = iRight.m_Coefficients;
206 }
207 return *this;
208 }
209
210 Self
211 operator+(const Self & iRight) const
212 {
213 return Self(this->m_Coefficients + iRight.m_Coefficients);
214 }
215
216 Self &
217 operator+=(const Self & iRight)
218 {
219 this->m_Coefficients += iRight.m_Coefficients;
220 return *this;
221 }
222
223 Self
224 operator-(const Self & iRight) const
225 {
226 return Self(this->m_Coefficients - iRight.m_Coefficients);
227 }
228
229 Self &
230 operator-=(const Self & iRight)
231 {
232 this->m_Coefficients -= iRight.m_Coefficients;
233 return *this;
234 }
235
236 Self
237 operator*(const CoordType & iV) const
238 {
239 Self oElement(this->m_Coefficients * iV);
240
241 return oElement;
242 }
243
244 Self &
246 {
247 this->m_Coefficients *= iV;
248 return *this;
249 }
250
251protected:
255 unsigned int m_Rank;
258 // bool m_MatrixFilled;
259
260 void
262 {
263 unsigned int k(0);
264 for (unsigned int dim1 = 0; dim1 < PointDimension; ++dim1)
265 {
266 for (unsigned int dim2 = dim1; dim2 < PointDimension; ++dim2)
267 {
268 m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++];
269 }
270 m_B[dim1] = -m_Coefficients[k++];
271 }
272 // m_MatrixFilled = true;
273 }
274};
275} // namespace itk
276#endif
vnl_vector_fixed< CoordType, Self::PointDimension > VNLVectorType
PointType ComputeOptimalLocation(const unsigned int)
TODO to be implemented!!!
CoordType ComputeErrorAtOptimalLocation(const PointType &iP)
TODO this method should be really optimized!!!
CoordType ComputeError(const PointType &iP) const
TODO this method should be really optimized!!!
vnl_vector_fixed< CoordType, Self::NumberOfCoefficients > CoefficientVectorType
QuadEdgeMeshDecimationQuadricElementHelper(const CoefficientVectorType &iCoefficients)
void AddPoint(const PointType &iP, const VectorType &iN, const CoordType &iWeight=static_cast< CoordType >(1.))
void AddTriangle(const PointType &iP1, const PointType &iP2, const PointType &iP3, const CoordType &iWeight=static_cast< CoordType >(1.))
A convenience class for computation of various triangle elements in 2D or 3D.
static VectorType ComputeNormal(const PointType &iA, const PointType &iB, const PointType &iC)
The "itk" namespace contains all Insight Segmentation and Registration Toolkit (ITK) classes....