Go to the documentation of this file.00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018 #ifndef __itkQuadEdgeMeshDecimationQuadricElementHelper_h
00019 #define __itkQuadEdgeMeshDecimationQuadricElementHelper_h
00020
00021 #include <itkPoint.h>
00022 #include <vnl/vnl_vector_fixed.h>
00023 #include <vnl/vnl_matrix.h>
00024 #include <vnl/algo/vnl_matrix_inverse.h>
00025
00026 #include "itkTriangleHelper.h"
00027 #include "itkNumericTraits.h"
00028
00029 namespace itk
00030 {
00032 template< class TPoint >
00033 class QuadEdgeMeshDecimationQuadricElementHelper
00034 {
00035 public:
00036
00037 typedef QuadEdgeMeshDecimationQuadricElementHelper Self;
00038 typedef TPoint PointType;
00039 typedef typename PointType::CoordRepType CoordType;
00040
00041 itkStaticConstMacro(PointDimension, unsigned int, PointType::PointDimension);
00042 itkStaticConstMacro(NumberOfCoefficients, unsigned int,
00043 PointDimension * ( PointDimension + 1 ) / 2 + PointDimension + 1);
00044
00045 typedef typename PointType::VectorType VectorType;
00046 typedef vnl_matrix< CoordType > VNLMatrixType;
00047 typedef vnl_vector_fixed< CoordType,
00048 itkGetStaticConstMacro(PointDimension) > VNLVectorType;
00049 typedef vnl_vector_fixed< CoordType,
00050 itkGetStaticConstMacro(NumberOfCoefficients) > CoefficientVectorType;
00051 typedef TriangleHelper< PointType > TriangleType;
00052
00053
00054 QuadEdgeMeshDecimationQuadricElementHelper():
00055 m_Coefficients( itk::NumericTraits< CoordType >::Zero ),
00056 m_A( PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero ),
00057 m_B( itk::NumericTraits< CoordType >::Zero ),
00058 m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-6 ) ),
00059 m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) )
00060 {
00061 this->m_Rank = PointDimension;
00062 }
00063
00064 QuadEdgeMeshDecimationQuadricElementHelper( const CoefficientVectorType& iCoefficients ):
00065 m_Coefficients( iCoefficients ),
00066 m_A( PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero ),
00067 m_B( itk::NumericTraits< CoordType >::Zero ),
00068 m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-3 ) ),
00069 m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) )
00070 {
00071 this->m_Rank = PointDimension;
00072 this->ComputeAMatrixAndBVector();
00073 }
00074
00075 ~QuadEdgeMeshDecimationQuadricElementHelper()
00076 {}
00077
00078 CoefficientVectorType GetCoefficients( ) const
00079 {
00080 return this->m_Coefficients;
00081 }
00082
00083 VNLMatrixType GetAMatrix()
00084 {
00085 this->ComputeAMatrixAndBVector();
00086 return m_A;
00087 }
00088
00089 VNLVectorType GetBVector()
00090 {
00091 ComputeAMatrixAndBVector();
00092 return m_B;
00093 }
00094
00095 unsigned int GetRank() const
00096 {
00097 return m_Rank;
00098 }
00099
00101 inline CoordType ComputeError( const PointType& iP ) const
00102 {
00103
00104 vnl_svd< CoordType > svd( m_A, m_SVDAbsoluteThreshold );
00105 svd.zero_out_relative( m_SVDRelativeThreshold );
00106 CoordType oError = inner_product( iP.GetVnlVector(), svd.recompose() * iP.GetVnlVector() );
00107
00108 return this->m_Coefficients[ this->m_Coefficients.size() - 1] - oError;
00109
00110
00111
00112
00113
00114
00115
00116
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129
00130
00131
00132
00133
00134 }
00135
00137 inline CoordType ComputeErrorAtOptimalLocation( const PointType& iP )
00138 {
00139 PointType optimal_location = ComputeOptimalLocation( iP );
00140 return ComputeError( optimal_location );
00141 }
00142
00143 PointType ComputeOptimalLocation( const PointType& iP )
00144 {
00145 ComputeAMatrixAndBVector();
00146
00147 vnl_svd< CoordType > svd( m_A, m_SVDAbsoluteThreshold );
00148 svd.zero_out_relative( m_SVDRelativeThreshold );
00149
00150 m_Rank = svd.rank();
00151
00152 VNLVectorType y = m_B.as_vector() - m_A * iP.GetVnlVector();
00153
00154 VNLVectorType displacement = svd.solve( y );
00155 PointType oP;
00156
00157 for( unsigned int dim = 0; dim < PointDimension; dim++ )
00158 {
00159 oP[dim] = iP[dim] + displacement[dim];
00160 }
00161
00162 return oP;
00163 }
00164
00166 PointType ComputeOptimalLocation(
00167 const unsigned int& iNumberOfEigenValues )
00168 {
00169 }
00170
00171
00172
00173
00174
00175
00176
00177
00178
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188 void AddTriangle( const PointType& iP1,
00189 const PointType& iP2,
00190 const PointType& iP3,
00191 const CoordType& iWeight = static_cast< CoordType >( 1. ) )
00192 {
00193 VectorType N = TriangleType::ComputeNormal( iP1, iP2, iP3 );
00194 AddPoint( iP1, N, iWeight );
00195 }
00196
00197 void AddPoint( const PointType& iP,
00198 const VectorType& iN,
00199 const CoordType& iWeight = static_cast< CoordType >( 1. ) )
00200 {
00201 unsigned int k( 0 ), dim1, dim2;
00202
00203 CoordType d = -iN * iP.GetVectorFromOrigin();
00204
00205 for( dim1 = 0; dim1 < PointDimension; ++dim1 )
00206 {
00207 for( dim2 = dim1; dim2 < PointDimension; ++dim2 )
00208 {
00209 this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
00210 }
00211 this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
00212 }
00213
00214 this->m_Coefficients[k++] += iWeight * d * d;
00215 }
00216
00217
00218
00219
00220
00221
00222 Self& operator=( const Self& iRight )
00223 {
00224 this->m_Coefficients = iRight.m_Coefficients;
00225 return *this;
00226 }
00227
00228 Self operator+( const Self & iRight ) const
00229 {
00230 return Self( this->m_Coefficients + iRight.m_Coefficients );
00231 }
00232
00233 Self& operator+=( const Self& iRight )
00234 {
00235 this->m_Coefficients += iRight.m_Coefficients;
00236 return *this;
00237 }
00238
00239 Self operator-( const Self& iRight ) const
00240 {
00241 return Self( this->m_Coefficients - iRight.m_Coefficients );
00242 }
00243
00244 Self& operator-=( const Self& iRight )
00245 {
00246 this->m_Coefficients -= iRight.m_Coefficients;
00247 return *this;
00248 }
00249
00250 Self operator*( const CoordType& iV ) const
00251 {
00252 Self oElement = Self( this->m_Coefficients * iV );
00253 return oElement;
00254 }
00255
00256 Self& operator*=( const CoordType& iV )
00257 {
00258 this->m_Coefficients *= iV;
00259 return *this;
00260 }
00261
00262 protected:
00263
00264 CoefficientVectorType m_Coefficients;
00265 VNLMatrixType m_A;
00266 VNLVectorType m_B;
00267 unsigned int m_Rank;
00268 CoordType m_SVDAbsoluteThreshold;
00269 CoordType m_SVDRelativeThreshold;
00270
00271
00272 void ComputeAMatrixAndBVector()
00273 {
00274 unsigned int k( 0 ), dim1, dim2;
00275
00276 for( dim1 = 0; dim1 < PointDimension; ++dim1 )
00277 {
00278 for( dim2 = dim1; dim2 < PointDimension; ++dim2 )
00279 {
00280 m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++];
00281 }
00282 m_B[dim1] = - m_Coefficients[k++];
00283 }
00284
00285 }
00286
00287 };
00288
00289 }
00290 #endif
00291