Main Page   Groups   Namespace List   Class Hierarchy   Alphabetical List   Compound List   File List   Namespace Members   Compound Members   File Members   Concepts

itkQuadEdgeMeshDecimationQuadricElementHelper.h

Go to the documentation of this file.
00001 /*=========================================================================
00002 
00003   Program:   Insight Segmentation & Registration Toolkit
00004   Module:    $RCSfile: itkQuadEdgeMeshDecimationQuadricElementHelper.h,v $
00005   Language:  C++
00006   Date:      $Date: 2010-03-24 21:51:55 $
00007   Version:   $Revision: 1.8 $
00008 
00009   Copyright (c) Insight Software Consortium. All rights reserved.
00010   See ITKCopyright.txt or http://www.itk.org/HTML/Copyright.htm for details.
00011 
00012      This software is distributed WITHOUT ANY WARRANTY; without even
00013      the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
00014      PURPOSE.  See the above copyright notices for more information.
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     //     ComputeAMatrixAndBVector();
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     CoordType oError( 0. );
00111      
00112     std::vector< CoordType > pt( PointDimension + 1, 1. );
00113 
00114     unsigned int dim1( 0 ), dim2, k( 0 );
00115 
00116     while( dim1 < PointDimension )
00117       {
00118       pt[dim1] = iP[dim1];
00119       ++dim1;
00120       }
00121     
00122     for( dim1 = 0; dim1 < PointDimension + 1; ++dim1 )
00123       {
00124       oError += this->m_Coefficients[k++] * pt[dim1] * pt[dim1];
00125       
00126       for( dim2 = dim1 + 1; dim2 < PointDimension + 1; ++dim2 )
00127         {
00128         oError += 2. * this->m_Coefficients[k++] * pt[dim1] * pt[dim2];
00129         }
00130       }
00131     oError += this->m_Coefficients[k++];
00132 
00133     return oError;*/
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 //   PointType ComputeOptimalLocation( 
00172 //     const CoordType& iValue )
00173 //   {
00174 //     ComputeAMatrixAndBVector();
00175 //     vnl_svd< CoordType > svd( m_A );
00176 //     svd.zero.zero_out_relative( iValue );
00177 //     m_Rank = svd.rank();
00178 //
00179 //     VNLVectorType location = svd.solve( m_B );
00180 //     PointType oP;
00181 //
00182 //     for( unsigned int dim = 0; dim < PointDimension; dim++ )
00183 //       oP[dim] = location[dim];
00184 //
00185 //     return oP;
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   // operators
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   //bool                        m_MatrixFilled;
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     //m_MatrixFilled = true;
00285     }
00286   
00287 };
00288 
00289 }
00290 #endif
00291 

Generated at Fri Apr 16 19:21:04 2010 for ITK by doxygen 1.6.1 written by Dimitri van Heesch, © 1997-2000