ITK  4.1.0
Insight Segmentation and Registration Toolkit
itkQuadEdgeMeshDecimationQuadricElementHelper.h
Go to the documentation of this file.
00001 /*=========================================================================
00002  *
00003  *  Copyright Insight Software Consortium
00004  *
00005  *  Licensed under the Apache License, Version 2.0 (the "License");
00006  *  you may not use this file except in compliance with the License.
00007  *  You may obtain a copy of the License at
00008  *
00009  *         http://www.apache.org/licenses/LICENSE-2.0.txt
00010  *
00011  *  Unless required by applicable law or agreed to in writing, software
00012  *  distributed under the License is distributed on an "AS IS" BASIS,
00013  *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
00014  *  See the License for the specific language governing permissions and
00015  *  limitations under the License.
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 
00028 namespace itk
00029 {
00031 template< class TPoint >
00032 class ITK_EXPORT QuadEdgeMeshDecimationQuadricElementHelper
00033 {
00034 public:
00035 
00036   typedef QuadEdgeMeshDecimationQuadricElementHelper Self;
00037   typedef TPoint                                     PointType;
00038   typedef typename PointType::CoordRepType           CoordType;
00039 
00040   itkStaticConstMacro(PointDimension, unsigned int, PointType::PointDimension);
00041   itkStaticConstMacro(NumberOfCoefficients, unsigned int,
00042                       PointDimension * ( PointDimension + 1 ) / 2 + PointDimension + 1);
00043 
00044   typedef typename PointType::VectorType VectorType;
00045   typedef vnl_matrix< CoordType >        VNLMatrixType;
00046   typedef vnl_vector_fixed< CoordType,
00047                             itkGetStaticConstMacro(PointDimension) >            VNLVectorType;
00048   typedef vnl_vector_fixed< CoordType,
00049                             itkGetStaticConstMacro(NumberOfCoefficients) >      CoefficientVectorType;
00050   typedef TriangleHelper< PointType > TriangleType;
00051 
00052   // *****************************************************************
00053   QuadEdgeMeshDecimationQuadricElementHelper():
00054     m_Coefficients(itk::NumericTraits< CoordType >::Zero),
00055     m_A(PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero),
00056     m_B(itk::NumericTraits< CoordType >::Zero),
00057     m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-6 ) ),
00058     m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) )
00059   {
00060     this->m_Rank = PointDimension;
00061   }
00062 
00063   QuadEdgeMeshDecimationQuadricElementHelper(const CoefficientVectorType & iCoefficients):
00064     m_Coefficients(iCoefficients),
00065     m_A(PointDimension, PointDimension, itk::NumericTraits< CoordType >::Zero),
00066     m_B(itk::NumericTraits< CoordType >::Zero),
00067     m_SVDAbsoluteThreshold( static_cast< CoordType >( 1e-3 ) ),
00068     m_SVDRelativeThreshold( static_cast< CoordType >( 1e-3 ) )
00069   {
00070     this->m_Rank = PointDimension;
00071     this->ComputeAMatrixAndBVector();
00072   }
00073 
00074   ~QuadEdgeMeshDecimationQuadricElementHelper()
00075   {}
00076 
00077   CoefficientVectorType GetCoefficients() const
00078   {
00079     return this->m_Coefficients;
00080   }
00081 
00082   VNLMatrixType GetAMatrix()
00083   {
00084     this->ComputeAMatrixAndBVector();
00085     return m_A;
00086   }
00087 
00088   VNLVectorType GetBVector()
00089   {
00090     ComputeAMatrixAndBVector();
00091     return m_B;
00092   }
00093 
00094   unsigned int GetRank() const
00095   {
00096     return m_Rank;
00097   }
00098 
00100   inline CoordType ComputeError(const PointType & iP) const
00101   {
00102     //     ComputeAMatrixAndBVector();
00103     vnl_svd< CoordType > svd(m_A, m_SVDAbsoluteThreshold);
00104     svd.zero_out_relative(m_SVDRelativeThreshold);
00105     CoordType oError = inner_product( iP.GetVnlVector(), svd.recompose() * iP.GetVnlVector() );
00106 
00107     return this->m_Coefficients[this->m_Coefficients.size() - 1] - oError;
00108     /*
00109     CoordType oError( 0. );
00110 
00111     std::vector< CoordType > pt( PointDimension + 1, 1. );
00112 
00113     unsigned int dim1( 0 ), dim2, k( 0 );
00114 
00115     while( dim1 < PointDimension )
00116       {
00117       pt[dim1] = iP[dim1];
00118       ++dim1;
00119       }
00120 
00121     for( dim1 = 0; dim1 < PointDimension + 1; ++dim1 )
00122       {
00123       oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1];
00124 
00125       for( dim2 = dim1 + 1; dim2 < PointDimension + 1; ++dim2 )
00126         {
00127         oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2];
00128         }
00129       }
00130     oError += this->m_Coefficients[k++];
00131 
00132     return oError;*/
00133   }
00134 
00136   inline CoordType ComputeErrorAtOptimalLocation(const PointType & iP)
00137   {
00138     PointType optimal_location = ComputeOptimalLocation(iP);
00139 
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 & )
00168   {}
00169 
00170   void AddTriangle( const PointType & iP1,
00171                     const PointType & iP2,
00172                     const PointType & iP3,
00173                     const CoordType & iWeight = static_cast< CoordType >( 1. ) )
00174   {
00175     VectorType N = TriangleType::ComputeNormal(iP1, iP2, iP3);
00176 
00177     AddPoint(iP1, N, iWeight);
00178   }
00179 
00180   void AddPoint( const PointType & iP,
00181                  const VectorType & iN,
00182                  const CoordType & iWeight = static_cast< CoordType >( 1. ) )
00183   {
00184     unsigned int k(0), dim1, dim2;
00185 
00186     CoordType d = -iN *iP.GetVectorFromOrigin();
00187 
00188     for ( dim1 = 0; dim1 < PointDimension; ++dim1 )
00189       {
00190       for ( dim2 = dim1; dim2 < PointDimension; ++dim2 )
00191         {
00192         this->m_Coefficients[k++] += iWeight * iN[dim1] * iN[dim2];
00193         }
00194       this->m_Coefficients[k++] += iWeight * iN[dim1] * d;
00195       }
00196 
00197     this->m_Coefficients[k++] += iWeight * d * d;
00198   }
00199 
00200   // ***********************************************************************
00201   // operators
00202   Self & operator=(const Self & iRight)
00203   {
00204     this->m_Coefficients = iRight.m_Coefficients;
00205     return *this;
00206   }
00207 
00208   Self operator+(const Self & iRight) const
00209   {
00210     return Self(this->m_Coefficients + iRight.m_Coefficients);
00211   }
00212 
00213   Self & operator+=(const Self & iRight)
00214   {
00215     this->m_Coefficients += iRight.m_Coefficients;
00216     return *this;
00217   }
00218 
00219   Self operator-(const Self & iRight) const
00220   {
00221     return Self(this->m_Coefficients - iRight.m_Coefficients);
00222   }
00223 
00224   Self & operator-=(const Self & iRight)
00225   {
00226     this->m_Coefficients -= iRight.m_Coefficients;
00227     return *this;
00228   }
00229 
00230   Self operator*(const CoordType & iV) const
00231   {
00232     Self oElement = Self(this->m_Coefficients * iV);
00233 
00234     return oElement;
00235   }
00236 
00237   Self & operator*=(const CoordType & iV)
00238   {
00239     this->m_Coefficients *= iV;
00240     return *this;
00241   }
00242 
00243 protected:
00244 
00245   CoefficientVectorType m_Coefficients;
00246   VNLMatrixType         m_A;
00247   VNLVectorType         m_B;
00248   unsigned int          m_Rank;
00249   CoordType             m_SVDAbsoluteThreshold;
00250   CoordType             m_SVDRelativeThreshold;
00251   //bool                        m_MatrixFilled;
00252 
00253   void ComputeAMatrixAndBVector()
00254   {
00255     unsigned int k(0), dim1, dim2;
00256 
00257     for ( dim1 = 0; dim1 < PointDimension; ++dim1 )
00258       {
00259       for ( 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     //m_MatrixFilled = true;
00266   }
00267 };
00268 }
00269 #endif
00270