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
00046
00047
00048
00049 typedef typename PointType::VectorType VectorType;
00050 typedef vnl_matrix< CoordType > VNLMatrixType;
00051 typedef vnl_vector_fixed< CoordType,
00052 itkGetStaticConstMacro(PointDimension) > VNLVectorType;
00053 typedef vnl_vector_fixed< CoordType,
00054 itkGetStaticConstMacro(NumberOfCoefficients) > CoefficientVectorType;
00055 typedef TriangleHelper< PointType > TriangleType;
00056
00057
00058 QuadEdgeMeshDecimationQuadricElementHelper():
00059 m_Coefficients( itk::NumericTraits< CoordType >::Zero ),
00060 m_A( PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero ),
00061 m_B( itk::NumericTraits< CoordType >::Zero )
00062 {
00063 this->m_Rank = PointDimension;
00064 }
00065
00066 QuadEdgeMeshDecimationQuadricElementHelper( const CoefficientVectorType& iCoefficients ):
00067 m_Coefficients( iCoefficients ),
00068 m_A( PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero ),
00069 m_B( itk::NumericTraits< CoordType >::Zero )
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 )
00102 {
00103 CoordType oError( 0. );
00104
00105 std::vector< CoordType > pt( PointDimension + 1, 1. );
00106
00107 unsigned int dim1 = 0;
00108 unsigned int k =0;
00109
00110 while( dim1 < PointDimension )
00111 {
00112 pt[dim1] = iP[dim1];
00113 dim1++;
00114 }
00115
00116 for( dim1 = 0; dim1 < PointDimension + 1; dim1++ )
00117 {
00118 oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1];
00119
00120 for( unsigned int dim2 = dim1 + 1; dim2 < PointDimension + 1; dim2++ )
00121 {
00122 oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2];
00123 }
00124 }
00125 return oError;
00126 }
00127
00129 inline CoordType ComputeErrorAtOptimalLocation()
00130 {
00131 return ComputeError( ComputeOptimalLocation() );
00132 }
00133
00134 PointType ComputeOptimalLocation()
00135 {
00136 ComputeAMatrixAndBVector();
00137 vnl_svd< CoordType > svd( m_A );
00138 svd.zero_out_relative( 1e-10 );
00139 m_Rank = svd.rank();
00140
00141 VNLVectorType location = svd.solve( m_B.as_vector() );
00142 PointType oP;
00143
00144 for( unsigned int dim = 0; dim < PointDimension; dim++ )
00145 oP[dim] = location[dim];
00146
00147 return oP;
00148 }
00149
00151 PointType ComputeOptimalLocation(
00152 const unsigned int& iNumberOfEigenValues )
00153 {
00154 }
00155
00156 PointType ComputeOptimalLocation(
00157 const CoordType& iValue )
00158 {
00159 ComputeAMatrixAndBVector();
00160 vnl_svd< CoordType > svd( m_A );
00161 svd.zero.zero_out_relative( iValue );
00162 m_Rank = svd.rank();
00163
00164 VNLVectorType location = svd.solve( m_B );
00165 PointType oP;
00166
00167 for( unsigned int dim = 0; dim < PointDimension; dim++ )
00168 oP[dim] = location[dim];
00169
00170 return oP;
00171 }
00172
00173 void AddTriangle( const PointType& iP1,
00174 const PointType& iP2,
00175 const PointType& iP3,
00176 const CoordType& iWeight = static_cast< CoordType >( 1. ) )
00177 {
00178 AddPoint( iP1, TriangleType::ComputeNormal( iP1, iP2, iP3 ), iWeight );
00179 }
00180
00181 void AddPoint( const PointType& iP,
00182 const VectorType& iN,
00183 const CoordType& iWeight = static_cast< CoordType >( 1. ) )
00184 {
00185 unsigned int k = 0;
00186
00187 CoordType d = -iN * iP.GetVectorFromOrigin();
00188
00189 for( unsigned int dim1 = 0; dim1 < PointDimension; dim1++ )
00190 {
00191 for( unsigned int dim2 = dim1; dim2 < PointDimension; dim2++ )
00192 {
00193 this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
00194 }
00195 this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
00196 }
00197
00198 this->m_Coefficients[k++] += iWeight * d * d;
00199 }
00200
00201
00202
00203
00204
00205
00206 Self& operator=( const Self& iRight )
00207 {
00208 this->m_Coefficients = iRight.m_Coefficients;
00209 return *this;
00210 }
00211
00212 Self operator+( const Self & iRight ) const
00213 {
00214 return Self( this->m_Coefficients + iRight.m_Coefficients );
00215 }
00216
00217 Self& operator+=( const Self& iRight )
00218 {
00219 this->m_Coefficients += iRight.m_Coefficients;
00220 return *this;
00221 }
00222
00223 Self operator-( const Self& iRight ) const
00224 {
00225 return Self( this->m_Coefficients - iRight.m_Coefficients );
00226 }
00227
00228 Self& operator-=( const Self& iRight )
00229 {
00230 this->m_Coefficients -= iRight.m_Coefficients;
00231 return *this;
00232 }
00233
00234 Self operator*( const CoordType& iV ) const
00235 {
00236 Self oElement = Self( this->m_Coefficients * iV );
00237 return oElement;
00238 }
00239
00240 Self& operator*=( const CoordType& iV )
00241 {
00242 this->m_Coefficients *= iV;
00243 return *this;
00244 }
00245
00246 protected:
00247
00248 CoefficientVectorType m_Coefficients;
00249 VNLMatrixType m_A;
00250 VNLVectorType m_B;
00251 unsigned int m_Rank;
00252
00253 void ComputeAMatrixAndBVector()
00254 {
00255 unsigned int k = 0;
00256
00257 for( unsigned int dim1 = 0; dim1 < PointDimension; dim1++ )
00258 {
00259 for( unsigned int dim2 = dim1; dim2 < PointDimension; dim2++ )
00260 {
00261 m_A[dim1][dim2] = m_A[dim2][dim1] = m_Coefficients[k++];
00262 }
00263 m_B[dim1] = - m_Coefficients[k++];
00264 }
00265 }
00266
00267 };
00268
00269 }
00270 #endif
00271