Commit 3efe8402 authored by Nicolas Mellado's avatar Nicolas Mellado

Merge branch 'curvature-wip' into 'master'

Add curvature estimators, fix several warning and improve testing.

See merge request !1
parents 8eb5face cb680c91
/*
Copyright (C) 2014 Nicolas Mellado <nmellado0@gmail.com>
Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/
#ifndef _GRENAILLE_COVARIANCE_PLANE_FIT_
#define _GRENAILLE_COVARIANCE_PLANE_FIT_
#include <Eigen/Eigenvalues>
#include <Eigen/Eigenvalues>
#include "enums.h"
namespace Grenaille
......@@ -19,17 +19,13 @@ namespace Grenaille
/*!
\brief Plane fitting procedure using only points position
\note This procedure requires two passes to fit a plane
This class can also computes the surface variation measure introduced in
\cite Pauly:2002:PSSimplification. The solver used to analyse the covariance
\cite Pauly:2002:PSSimplification. The solver used to analyse the covariance
matrix is stored for further use.
\inherit Concept::FittingProcedureConcept
\warning This class is currently untested and should not be used !
\see CompactPlane
*/
template < class DataPoint, class _WFunctor, typename T >
......@@ -47,13 +43,13 @@ protected:
public:
/*! \brief Scalar type inherited from DataPoint*/
typedef typename Base::Scalar Scalar;
typedef typename Base::Scalar Scalar;
/*! \brief Vector type inherited from DataPoint*/
typedef typename Base::VectorType VectorType;
/*! \brief Vector type inherited from DataPoint*/
typedef typename Base::MatrixType MatrixType;
/*! \brief Weight Function*/
typedef _WFunctor WFunctor;
typedef _WFunctor WFunctor;
/*! \brief Solver used to analyse the covariance matrix*/
typedef Eigen::SelfAdjointEigenSolver<MatrixType> Solver;
......@@ -95,23 +91,23 @@ public:
/**************************************************************************/
/* Results */
/**************************************************************************/
using Base::potential;
/*! \brief Value of the scalar field at the evaluation point */
MULTIARCH inline Scalar potential() const { return Base::potential(m_evalPos); }
/*! \brief Value of the normal of the primitive at the evaluation point */
MULTIARCH inline VectorType normal() const { return Base::m_p.template head<DataPoint::Dim>(); }
/*! \brief Reading access to the Solver used to analyse the covariance
/*! \brief Reading access to the Solver used to analyse the covariance
matrix */
MULTIARCH inline const Solver& solver() const { return m_solver; }
/*! \brief Implements \cite Pauly:2002:PSSimplification surface variation.
It computes the ratio \f$ d \frac{\lambda_0}{\sum_i \lambda_i} \f$ with \c d the dimension of the ambient space.
\return 0 for invalid fits
*/
MULTIARCH inline Scalar surfaceVariation() const;
......@@ -122,20 +118,20 @@ namespace internal {
using ::Grenaille::internal::FitSpaceDer;
using ::Grenaille::internal::FitScaleDer;
/*!
/*!
\brief Internal generic class computing the derivatives of covariance plane fits
\inherit Concept::FittingExtensionConcept
The differentiation can be done automatically in scale and/or space, by
combining the enum values FitScaleDer and FitSpaceDer in the template
combining the enum values FitScaleDer and FitSpaceDer in the template
parameter Type.
The differenciated values are stored in static arrays. The size of the
arrays is computed with respect to the derivation type (scale and/or space)
and the number of the dimension of the ambiant space.
By convention, the scale derivatives are stored at index 0 when Type
and the number of the dimension of the ambiant space.
By convention, the scale derivatives are stored at index 0 when Type
contains at least FitScaleDer. The size of these arrays can be known using
derDimension(), and the differentiation type by isScaleDer() and
derDimension(), and the differentiation type by isScaleDer() and
isSpaceDer().
*/
template < class DataPoint, class _WFunctor, typename T, int Type>
......@@ -143,7 +139,7 @@ class CovariancePlaneDer : public T
{
private:
typedef T Base; /*!< \brief Generic base type */
protected:
enum
......@@ -152,7 +148,7 @@ protected:
PROVIDES_COVARIANCE_PLANE_DERIVATIVE, /*!< \brief Provides derivatives for hyper-planes */
PROVIDES_NORMAL_DERIVATIVE
};
static const int NbDerivatives = ((Type & FitScaleDer) ? 1 : 0 ) + ((Type & FitSpaceDer) ? DataPoint::Dim : 0);
static const int DerStorageOrder = (Type & FitSpaceDer) ? Eigen::RowMajor : Eigen::ColMajor;
......@@ -171,7 +167,7 @@ private:
// computation data
ScalarArray m_dSumW; /*!< \brief Sum of weight derivatives */
MatrixType m_dCov[NbDerivatives];
VectorArray m_dCog; /*!< \brief Derivatives of the centroid */
VectorArray m_dNormal; /*!< \brief Derivatives of the hyper-plane normal */
ScalarArray m_dDist; /*!< \brief Derivatives of the MLS scalar field */
......@@ -196,10 +192,10 @@ public:
/**************************************************************************/
/* Use results */
/**************************************************************************/
/*! \brief Returns the derivatives of the scalar field at the evaluation point */
MULTIARCH inline ScalarArray dPotential() const { return m_dDist; }
/*! \brief Returns the derivatives of the primitive normal */
MULTIARCH inline VectorArray dNormal() const { return m_dNormal; }
......@@ -212,15 +208,15 @@ public:
}; //class CovariancePlaneDer
}// namespace internal
}// namespace internal
/*!
\brief Differentiation in scale of the CovariancePlaneFit
\inherit Concept::FittingExtensionConcept
Requierement:
Requierement:
\verbatim PROVIDES_COVARIANCE_PLANE \endverbatim
Provide:
Provide:
\verbatim PROVIDES_COVARIANCE_PLANE_SCALE_DERIVATIVE \endverbatim
*/
template < class DataPoint, class _WFunctor, typename T>
......@@ -237,9 +233,9 @@ protected:
\brief Spatial differentiation of the CovariancePlaneFit
\inherit Concept::FittingExtensionConcept
Requierement:
Requierement:
\verbatim PROVIDES_COVARIANCE_PLANE \endverbatim
Provide:
Provide:
\verbatim PROVIDES_COVARIANCE_PLANE_SPACE_DERIVATIVE \endverbatim
*/
template < class DataPoint, class _WFunctor, typename T>
......@@ -256,10 +252,10 @@ protected:
\brief Differentiation both in scale and space of the CovariancePlaneFit
\inherit Concept::FittingExtensionConcept
Requierement:
Requierement:
\verbatim PROVIDES_COVARIANCE_PLANE \endverbatim
Provide:
\verbatim PROVIDES_COVARIANCE_PLANE_SCALE_DERIVATIVE
Provide:
\verbatim PROVIDES_COVARIANCE_PLANE_SCALE_DERIVATIVE
PROVIDES_COVARIANCE_PLANE_SPACE_DERIVATIVE
\endverbatim
*/
......
/*
Copyright (C) 2014 Nicolas Mellado <nmellado0@gmail.com>
Copyright (C) 2015 Gael Guennebaud <gael.guennebaud@inria.fr>
This Source Code Form is subject to the terms of the Mozilla Public
License, v. 2.0. If a copy of the MPL was not distributed with this
file, You can obtain one at http://mozilla.org/MPL/2.0/.
file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/
template < class DataPoint, class _WFunctor, typename T>
void
void
CovariancePlaneFit<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
// Setup primitive
......@@ -24,7 +24,7 @@ CovariancePlaneFit<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
}
template < class DataPoint, class _WFunctor, typename T>
bool
bool
CovariancePlaneFit<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
{
VectorType q = _nei.pos() - m_evalPos;
......@@ -36,7 +36,7 @@ CovariancePlaneFit<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
m_cog += w * q;
m_sumW += w;
m_cov += w * q * q.transpose();
++(Base::m_nbNeighbors);
return true;
}
......@@ -48,8 +48,6 @@ template < class DataPoint, class _WFunctor, typename T>
FIT_RESULT
CovariancePlaneFit<DataPoint, _WFunctor, T>::finalize ()
{
Scalar epsilon = Eigen::NumTraits<Scalar>::dummy_precision();
// handle specific configurations
// With less than 3 neighbors the fitting is undefined
if(m_sumW == Scalar(0.) || Base::m_nbNeighbors < 3)
......@@ -61,15 +59,20 @@ CovariancePlaneFit<DataPoint, _WFunctor, T>::finalize ()
// Center the covariance on the centroid:
m_cov -= m_cog * m_cog.transpose() / m_sumW;
// \note The covariance matrix should be here normalized by m_sumW.
// As it does not affect the eigen decomposition, we skip this normalization
// to save computation.
// Finalize the centroid
m_cog = m_cog/m_sumW + m_evalPos;
#ifdef __CUDACC__
m_solver.computeDirect(m_cov);
#else
m_solver.compute(m_cov);
#endif
Base::setPlane(m_solver.eigenvectors().col(0), m_cog);
// \todo Use the output of the solver to check stability
......@@ -93,7 +96,7 @@ namespace internal
{
template < class DataPoint, class _WFunctor, typename T, int Type>
void
void
CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::init(const VectorType& _evalPos)
{
Base::init(_evalPos);
......@@ -106,7 +109,7 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::init(const VectorType& _evalP
template < class DataPoint, class _WFunctor, typename T, int Type>
bool
bool
CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::addNeighbor(const DataPoint &_nei)
{
bool bResult = Base::addNeighbor(_nei);
......@@ -140,7 +143,7 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::addNeighbor(const DataPoint
template < class DataPoint, class _WFunctor, typename T, int Type>
FIT_RESULT
FIT_RESULT
CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::finalize()
{
MULTIARCH_STD_MATH(sqrt);
......@@ -155,7 +158,7 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::finalize()
Eigen::Matrix<Scalar,2,1> shifted_eivals = Base::m_solver.eigenvalues().template tail<2>().array() - Base::m_solver.eigenvalues()(0);
if(shifted_eivals(0) < consider_as_zero || shifted_eivals(0) < epsilon * shifted_eivals(1)) shifted_eivals(0) = 0;
if(shifted_eivals(1) < consider_as_zero) shifted_eivals(1) = 0;
for(int k=0; k<NbDerivatives; ++k)
{
// Finalize the computation of dCov.
......@@ -167,10 +170,10 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::finalize()
- shifted_cog * m_dCog.col(k).transpose()
- m_dCog.col(k) * shifted_cog.transpose()
+ m_dSumW[k] * shifted_cog * shifted_cog.transpose();
// cancel centered basis of dCog:
m_dCog.col(k) += m_dSumW[k] * Base::m_evalPos;
// apply normalization by sumW:
m_dCog.col(k) = (m_dCog.col(k) - m_dSumW(k) * Base::m_cog) / Base::m_sumW;
......@@ -187,7 +190,7 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::finalize()
if(shifted_eivals(0)>0) z(0) /= shifted_eivals(0);
if(shifted_eivals(1)>0) z(1) /= shifted_eivals(1);
m_dNormal.col(k) = Base::m_solver.eigenvectors().template rightCols<2>() * z;
VectorType dDiff = -m_dCog.col(k);
if(k>0 || !isScaleDer())
dDiff(isScaleDer() ? k-1 : k) += 1;
......
......@@ -5,6 +5,7 @@
*/
#include <Eigen/Eigenvalues>
#include <Eigen/Core>
#ifndef _GRENAILLE_CURVATURE_
#define _GRENAILLE_CURVATURE_
......@@ -45,6 +46,10 @@ public:
typedef typename Base::VectorType VectorType; /*!< \brief Inherited vector type*/
typedef typename DataPoint::MatrixType MatrixType; /*!< \brief Matrix type inherited from DataPoint*/
private:
typedef Eigen::Matrix<Scalar,3,2> Mat32; /*!< \brief Matrix type for tangent plane basis */
typedef Eigen::Matrix<Scalar,2,2> Mat22; /*!< \brief Matrix type for shape operator */
private:
Scalar m_k1, m_k2;
VectorType m_v1, m_v2;
......@@ -87,6 +92,20 @@ public:
//! \brief Returns an estimate of the Gaussian curvature
MULTIARCH inline Scalar GaussianCurvature() const { return m_k1 * m_k2;}
//! \brief Compute principal curvature directions
//!
//! The tangent plane can be calculated from the normal vector or from its
//! derivatives, depending of the useNormal parameter
//!
//! The finalize() method calls this function with useNormal=false by
//! default.
//!
MULTIARCH inline FIT_RESULT computeCurvature(bool useNormal = false);
protected:
//! \brief Compute a tangent plane basis
MULTIARCH inline Mat32 tangentPlane(bool useNormal = false) const;
};
#include "curvature.hpp"
......
......@@ -6,31 +6,111 @@
file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/
template < class DataPoint, class _WFunctor, typename T>
FIT_RESULT
CurvatureEstimator<DataPoint, _WFunctor, T>::finalize()
{
FIT_RESULT bResult = Base::finalize();
if(bResult != UNDEFINED)
{
return computeCurvature(false);
}
return bResult;
}
template < class DataPoint, class _WFunctor, typename T>
FIT_RESULT CurvatureEstimator<DataPoint, _WFunctor, T>::computeCurvature(bool useNormal)
{
// Get the object space Weingarten map dN
MatrixType dN = Base::dNormal().template middleCols<DataPoint::Dim>(Base::isScaleDer() ? 1: 0);
// Compute tangent-space basis
Mat32 B = tangentPlane(useNormal);
// Compute the 2x2 matrix representing the shape operator by transforming dN to the basis B.
// Recall that dN is a bilinear form, it thus transforms as follows:
Mat22 S = B.transpose() * dN * B;
// Recall that at this stage, the shape operator represented by S describes the normal curvature K_n(u) in the direction u \in R^2 as follows:
// K_n(u) = u^T S u
// The principal curvatures are fully defined by the values and the directions of the extrema of K_n.
//
// If the normal field N(x) comes from the gradient of a scalar field, then N(x) is curl-free, and dN and S are symmetric matrices.
// In this case, the extrema of the previous quadratic form are directly obtained by the the eigenvalue decomposition of S.
// However, if N(x) is only an approximation of the normal field of a surface, then N(x) is not necessarily curl-free, and in this case S is not symmetric.
// In this case, we first have to find an equivalent symmetric matrix S' such that:
// K_n(u) = u^T S' u,
// for any u \in R^2.
// It is easy to see that such a S' is simply obtained as:
// S' = (S + S^T)/2
S(0,1) = S(1,0) = (S(0,1) + S(1,0))/Scalar(2);
Eigen::SelfAdjointEigenSolver<Mat22> eig2;
eig2.computeDirect(S);
if (eig2.info() != Eigen::Success){
return UNDEFINED;
}
m_k1 = eig2.eigenvalues()(0);
m_k2 = eig2.eigenvalues()(1);
m_v1 = B * eig2.eigenvectors().col(0);
m_v2 = B * eig2.eigenvectors().col(1);
if(abs(m_k1)<abs(m_k2))
{
#ifdef __CUDACC__
Scalar tmpk = m_k1;
m_k1 = m_k2;
m_k2 = tmpk;
VectorType tmpv = m_v1;
m_v1 = m_v2;
m_v2 = tmpv;
#else
std::swap(m_k1, m_k2);
std::swap(m_v1, m_v2);
#endif
}
return this->getCurrentState();
}
template < class DataPoint, class _WFunctor, typename T>
typename CurvatureEstimator<DataPoint, _WFunctor, T>::Mat32
CurvatureEstimator<DataPoint, _WFunctor, T>::tangentPlane(bool useNormal) const
{
typedef typename VectorType::Index Index;
typedef Eigen::Matrix<Scalar,3,2> Mat32;
typedef Eigen::Matrix<Scalar,2,2> Mat22;
MULTIARCH_STD_MATH(sqrt);
MULTIARCH_STD_MATH(abs);
FIT_RESULT bResult = Base::finalize();
Mat32 B;
Index i0=Index(-1), i1=Index(-1), i2=Index(-1);
if(bResult != UNDEFINED)
if(useNormal)
{
VectorType n = Base::normal();
n.minCoeff(&i0);
i1 = (i0+1)%3;
i2 = (i0+2)%3;
B.col(0)[i0] = 0;
B.col(0)[i1] = n[i2];
B.col(0)[i2] = -n[i1];
B.col(0).normalize();
B.col(1) = B.col(0).cross(n);
}
else
{
// Get the object space Weingarten map dN
MatrixType dN = Base::dNormal().template middleCols<DataPoint::Dim>(Base::isScaleDer() ? 1: 0);
// Compute tangent-space basis from dN
// 1 - pick the column with maximal norm as the first tangent vector,
Index i0, i1, i2;
Scalar sqNorm = dN.colwise().squaredNorm().maxCoeff(&i0);
Mat32 B;
B.col(0) = dN.col(i0) / sqrt(sqNorm);
// 2 - orthogonalize the other column vectors, and pick the most reliable one
i1 = (i0+1)%3;
......@@ -41,53 +121,8 @@ CurvatureEstimator<DataPoint, _WFunctor, T>::finalize()
Scalar v2norm2 = v2.squaredNorm();
if(v1norm2 > v2norm2) B.col(1) = v1 / sqrt(v1norm2);
else B.col(1) = v2 / sqrt(v2norm2);
// Compute the 2x2 matrix representing the shape operator by transforming dN to the basis B.
// Recall that dN is a bilinear form, it thus transforms as follows:
Mat22 S = B.transpose() * dN * B;
// Recall that at this stage, the shape operator represented by S describes the normal curvature K_n(u) in the direction u \in R^2 as follows:
// K_n(u) = u^T S u
// The principal curvatures are fully defined by the values and the directions of the extrema of K_n.
//
// If the normal field N(x) comes from the gradient of a scalar field, then N(x) is curl-free, and dN and S are symmetric matrices.
// In this case, the extrema of the previous quadratic form are directly obtained by the the eigenvalue decomposition of S.
// However, if N(x) is only an approximation of the normal field of a surface, then N(x) is not necessarily curl-free, and in this case S is not symmetric.
// In this case, we first have to find an equivalent symmetric matrix S' such that:
// K_n(u) = u^T S' u,
// for any u \in R^2.
// It is easy to see that such a S' is simply obtained as:
// S' = (S + S^T)/2
S(0,1) = S(1,0) = (S(0,1) + S(1,0))/Scalar(2);
Eigen::SelfAdjointEigenSolver<Mat22> eig2;
eig2.computeDirect(S);
if (eig2.info() != Eigen::Success){
return UNDEFINED;
}
m_k1 = eig2.eigenvalues()(0);
m_k2 = eig2.eigenvalues()(1);
m_v1 = B * eig2.eigenvectors().col(0);
m_v2 = B * eig2.eigenvectors().col(1);
if(abs(m_k1)<abs(m_k2))
{
#ifdef __CUDACC__
Scalar tmpk = m_k1;
m_k1 = m_k2;
m_k2 = tmpk;
VectorType tmpv = m_v1;
m_v1 = m_v2;
m_v2 = tmpv;
#else
std::swap(m_k1, m_k2);
std::swap(m_v1, m_v2);
#endif
}
}
return bResult;
return B;
}

#ifndef _GRENAILLE_CURVATUREESTIMATION_
#define _GRENAILLE_CURVATUREESTIMATION_
namespace Grenaille
{
/** \addtogroup experimental
* @{
*/
template < class DataPoint, class _WFunctor, typename T>
class BaseCurvatureEstimator : public T
{
private:
typedef T Base;
protected:
enum
{
PROVIDES_PRINCIPALE_CURVATURES
};
public:
typedef typename Base::Scalar Scalar; /*!< \brief Inherited scalar type*/
typedef typename Base::VectorType VectorType; /*!< \brief Inherited vector type*/
typedef typename DataPoint::MatrixType MatrixType; /*!< \brief Matrix type inherited from DataPoint*/
protected:
Scalar m_k1, m_k2;
VectorType m_v1, m_v2;
public:
/*! \brief Default constructor */
MULTIARCH inline BaseCurvatureEstimator() : m_k1(0), m_k2(0) {}
public:
/**************************************************************************/
/* Initialization */
/**************************************************************************/
/*! \copydoc Concept::FittingProcedureConcept::init() */
MULTIARCH inline void init (const VectorType& _evalPos);
/**************************************************************************/
/* Use results */
/**************************************************************************/
//! \brief Returns an estimate of the first principal curvature value
//!
//! It is the greatest curvature in <b>absolute value</b>.
MULTIARCH inline Scalar k1() const { return m_k1; }
//! \brief Returns an estimate of the second principal curvature value
//!
//! It is the smallest curvature in <b>absolute value</b>.
MULTIARCH inline Scalar k2() const { return m_k2; }
//! \brief Returns an estimate of the first principal curvature direction
//!
//! It is the greatest curvature in <b>absolute value</b>.
MULTIARCH inline VectorType k1Direction() const { return m_v1; }
//! \brief Returns an estimate of the second principal curvature direction
//!
//! It is the smallest curvature in <b>absolute value</b>.
MULTIARCH inline VectorType k2Direction() const { return m_v2; }
//! \brief Returns an estimate of the mean curvature
MULTIARCH inline Scalar kMean() const { return (m_k1 + m_k2)/2.;}
//! \brief Returns an estimate of the Gaussian curvature
MULTIARCH inline Scalar GaussianCurvature() const { return m_k1 * m_k2;}
};
/*!
* \brief Extension to compute curvature values based on a covariance analysis
* of normal vectors of neighbors.
*
* A 3D covariance matrix is computed from the normals of the neighbors and the
* two principal curvature values and directions are given by the two extreme
* eigenvalues and associated eigenvectors of the covariance matrix
* \cite Liang:1990:RRSS.
*
* \warning This class is valid only in 3D.
*/
template < class DataPoint, class _WFunctor, typename T>
class NormalCovarianceCurvature : public BaseCurvatureEstimator<DataPoint,_WFunctor,T>
{
private:
typedef BaseCurvatureEstimator<DataPoint,_WFunctor,T> Base;
//TODO(thib) check the curvature values that might be wrong
//TODO(thib) use weighting function
//TODO(thib) which eigenvectors should be selected ? extreme of maximal ?
public:
typedef typename Base::Scalar Scalar; /*!< \brief Inherited scalar type*/
typedef typename Base::VectorType VectorType; /*!< \brief Inherited vector type*/
typedef typename DataPoint::MatrixType MatrixType; /*!< \brief Matrix type inherited from DataPoint*/
/*! \brief Solver used to analyse the covariance matrix*/
typedef Eigen::SelfAdjointEigenSolver<MatrixType> Solver;
protected:
MatrixType m_cov; /*!< \brief Covariance matrix */
VectorType m_cog; /*!< \brief Gravity center */
Solver m_solver; /*!< \brief Solver used to analyse the covariance matrix */
public:
/*! \brief Default constructor */
MULTIARCH inline NormalCovarianceCurvature() : Base() {}
public:
/**************************************************************************/
/* Initialization */
/**************************************************************************/
/*! \copydoc Concept::FittingProcedureConcept::init() */
MULTIARCH inline void init (const VectorType& _evalPos);
/**************************************************************************/
/* Processing */
/**************************************************************************/
/*! \copydoc Concept::FittingProcedureConcept::addNeighbor() */
MULTIARCH inline bool addNeighbor(const DataPoint &_nei);
/*! \copydoc Concept::FittingProcedureConcept::finalize() */
MULTIARCH inline FIT_RESULT finalize();
};
/*!
* \brief Extension to compute curvature values based on a covariance analysis
* of normal vectors of neighbors projected onto the tangent plane.
*
* A 2D covariance matrix is computed from the projections of normals of the
* neighbors and the two principal curvature values and directions are given by
* the eigenvalues and associated eigenvectors of the covariance matrix
* \cite Berkmann:1994:CSG.
*
* \note This procedure requires two passes, the first one for plane fitting,
* the second one for covariance analysis.
* \warning This class is valid only in 3D.
*/
template < class DataPoint, class _WFunctor, typename T>
class ProjectedNormalCovarianceCurvature : public BaseCurvatureEstimator<DataPoint,_WFunctor,T>
{
private:
typedef BaseCurvatureEstimator<DataPoint,_WFunctor,T> Base;
//TODO(thib) check the curvature values that might be wrong
//TODO(thib) use weighting function
protected:
enum
{
Check = Base::PROVIDES_PLANE
};
enum PASS : int
{
FIRST_PASS = 0,
SECOND_PASS,
PASS_COUNT
};
public:
typedef typename Base::Scalar Scalar; /*!< \brief Inherited scalar type*/
typedef typename Base::VectorType VectorType; /*!< \brief Inherited vector type*/
typedef typename DataPoint::MatrixType MatrixType; /*!< \brief Matrix type inherited from DataPoint*/
//TODO(thib) use of Eigen::RowAtCompileTime-1 ?
typedef Eigen::Matrix<Scalar,2,2> Mat22;
typedef Eigen::Matrix<Scalar,3,2> Mat32;
typedef Eigen::Matrix<Scalar,2,1> Vector2;
typedef typename VectorType::Index Index;
/*! \brief Solver used to analyse the covariance matrix*/
typedef Eigen::SelfAdjointEigenSolver<Mat22> Solver;
protected:
Vector2 m_cog; /*!< \brief Gravity center */
Mat22 m_cov; /*!< \brief Covariance matrix */
Solver m_solver; /*!< \brief Solver used to analyse the covariance matrix */
PASS m_pass; /*!< \brief Current pass */
Mat32 m_tframe; /*!< \brief Tangent frame */