18 #ifndef __itkQuadEdgeMeshDecimationQuadricElementHelper_h
19 #define __itkQuadEdgeMeshDecimationQuadricElementHelper_h
22 #include "vnl/vnl_vector_fixed.h"
23 #include "vnl/vnl_matrix.h"
24 #include "vnl/algo/vnl_matrix_inverse.h"
31 template<
class TPo
int >
40 itkStaticConstMacro(PointDimension,
unsigned int, PointType::PointDimension);
41 itkStaticConstMacro(NumberOfCoefficients,
unsigned int,
42 PointDimension * ( PointDimension + 1 ) / 2 + PointDimension + 1);
57 m_SVDAbsoluteThreshold( static_cast<
CoordType >( 1
e-6 ) ),
58 m_SVDRelativeThreshold( static_cast<
CoordType >( 1
e-3 ) )
60 this->m_Rank = PointDimension;
64 m_Coefficients(iCoefficients),
67 m_SVDAbsoluteThreshold( static_cast<
CoordType >( 1
e-3 ) ),
68 m_SVDRelativeThreshold( static_cast<
CoordType >( 1
e-3 ) )
70 this->m_Rank = PointDimension;
71 this->ComputeAMatrixAndBVector();
79 return this->m_Coefficients;
84 this->ComputeAMatrixAndBVector();
90 ComputeAMatrixAndBVector();
94 unsigned int GetRank()
const
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() );
107 return this->m_Coefficients[this->m_Coefficients.size() - 1] - oError;
138 PointType optimal_location = ComputeOptimalLocation(iP);
140 return ComputeError(optimal_location);
145 ComputeAMatrixAndBVector();
147 vnl_svd< CoordType > svd(m_A, m_SVDAbsoluteThreshold);
148 svd.zero_out_relative(m_SVDRelativeThreshold);
157 for (
unsigned int dim = 0; dim < PointDimension; dim++ )
159 oP[dim] = iP[dim] + displacement[dim];
167 const unsigned int & )
173 const CoordType & iWeight = static_cast< CoordType >( 1. ) )
175 VectorType N = TriangleType::ComputeNormal(iP1, iP2, iP3);
177 AddPoint(iP1, N, iWeight);
182 const CoordType & iWeight = static_cast< CoordType >( 1. ) )
184 unsigned int k(0), dim1, dim2;
186 CoordType d = -iN *iP.GetVectorFromOrigin();
188 for ( dim1 = 0; dim1 < PointDimension; ++dim1 )
190 for ( dim2 = dim1; dim2 < PointDimension; ++dim2 )
192 this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
194 this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
197 this->m_Coefficients[k++] += iWeight * d * d;
232 Self oElement =
Self(this->m_Coefficients * iV);
239 this->m_Coefficients *= iV;
253 void ComputeAMatrixAndBVector()
255 unsigned int k(0), dim1, dim2;
257 for ( dim1 = 0; dim1 < PointDimension; ++dim1 )
259 for ( dim2 = dim1; dim2 < PointDimension; ++dim2 )
261 m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++];
263 m_B[dim1] = -m_Coefficients[k++];