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: 2008-10-01 18:36:29 $
00007   Version:   $Revision: 1.5 $
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 //    static const unsigned int PointDimension = PointType::PointDimension;
00046 //    static const unsigned int NumberOfCoefficients = 
00047 //      PointDimension * ( PointDimension + 1 ) / 2 + PointDimension + 1;
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   // operators
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 

Generated at Tue Sep 15 04:25:46 2009 for ITK by doxygen 1.5.8 written by Dimitri van Heesch, © 1997-2000