Commit 5d6e7bdf authored by Thibault Lejemble's avatar Thibault Lejemble

restore inheritance to a common class for curvature estimation

parent a5208620
......@@ -5,54 +5,69 @@
namespace Grenaille
{
#define CURVATURE_EXTENSION \
protected: \
enum \
{ \
PROVIDES_PRINCIPALE_CURVATURES \
}; \
\
protected: \
Scalar m_k1, m_k2; \
VectorType m_v1, m_v2; \
\
public: \
/**************************************************************************/\
/* 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;} \
private: \
/*TODO(thib) is it the right way?*/ \
MULTIARCH inline void reset() { \
m_k1 = Scalar(0); \
m_k2 = Scalar(0); \
m_v1 = VectorType::Zero(); \
m_v2 = VectorType::Zero(); \
}
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
......@@ -66,10 +81,10 @@ private:
* \warning This class is valid only in 3D.
*/
template < class DataPoint, class _WFunctor, typename T>
class NormalCovarianceCurvature : public T
class NormalCovarianceCurvature : public BaseCurvatureEstimator<DataPoint,_WFunctor,T>
{
private:
typedef T Base;
typedef BaseCurvatureEstimator<DataPoint,_WFunctor,T> Base;
//TODO(thib) check the curvature values that might be wrong
//TODO(thib) use weighting function
......@@ -82,8 +97,6 @@ public:
/*! \brief Solver used to analyse the covariance matrix*/
typedef Eigen::SelfAdjointEigenSolver<MatrixType> Solver;
CURVATURE_EXTENSION
protected:
MatrixType m_cov; /*!< \brief Covariance matrix */
VectorType m_cog; /*!< \brief Gravity center */
......@@ -125,10 +138,10 @@ public:
* \warning This class is valid only in 3D.
*/
template < class DataPoint, class _WFunctor, typename T>
class ProjectedNormalCovarianceCurvature : public T
class ProjectedNormalCovarianceCurvature : public BaseCurvatureEstimator<DataPoint,_WFunctor,T>
{
private:
typedef T Base;
typedef BaseCurvatureEstimator<DataPoint,_WFunctor,T> Base;
//TODO(thib) check the curvature values that might be wrong
//TODO(thib) use weighting function
......@@ -158,8 +171,6 @@ public:
/*! \brief Solver used to analyse the covariance matrix*/
typedef Eigen::SelfAdjointEigenSolver<Mat22> Solver;
CURVATURE_EXTENSION
protected:
Vector2 m_cog; /*!< \brief Gravity center */
Mat22 m_cov; /*!< \brief Covariance matrix */
......

template < class DataPoint, class _WFunctor, typename T>
void
BaseCurvatureEstimator<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
Base::init(_evalPos);
m_k1 = 0;
m_k2 = 0;
m_v1 = VectorType::Zero();
m_v2 = VectorType::Zero();
}
template < class DataPoint, class _WFunctor, typename T>
void
NormalCovarianceCurvature<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
Base::init(_evalPos);
reset();
m_cov = MatrixType::Zero();
m_cog = VectorType::Zero();
}
......@@ -45,19 +57,19 @@ NormalCovarianceCurvature<DataPoint, _WFunctor, T>::finalize ()
m_solver.computeDirect(m_cov);
m_k1 = m_solver.eigenvalues()(1);
m_k2 = m_solver.eigenvalues()(2);
Base::m_k1 = m_solver.eigenvalues()(1);
Base::m_k2 = m_solver.eigenvalues()(2);
m_v1 = m_solver.eigenvectors().col(1);
m_v2 = m_solver.eigenvectors().col(2);
Base::m_v1 = m_solver.eigenvectors().col(1);
Base::m_v2 = m_solver.eigenvectors().col(2);
//TODO(thib) which epsilon value should be chosen ?
// Scalar epsilon = Eigen::NumTraits<Scalar>::dummy_precision();
Scalar epsilon = Scalar(1e-3);
if(m_k1<epsilon && m_k2<epsilon)
if(Base::m_k1<epsilon && Base::m_k2<epsilon)
{
m_k1 = Scalar(0);
m_k2 = Scalar(0);
Base::m_k1 = Scalar(0);
Base::m_k2 = Scalar(0);
// set principal directions from normals center of gravity
VectorType n = m_cog.normalized();
......@@ -65,12 +77,12 @@ NormalCovarianceCurvature<DataPoint, _WFunctor, T>::finalize ()
n.rowwise().squaredNorm().minCoeff(&i0);
i1 = (i0+1)%3;
i2 = (i0+2)%3;
m_v1[i0] = 0;
m_v1[i1] = n[i2];
m_v1[i2] = -n[i1];
m_v2[i0] = n[i1]*n[i1] + n[i2]*n[i2];
m_v2[i1] = -n[i1]*n[i0];
m_v2[i2] = -n[i2]*n[i0];
Base::m_v1[i0] = 0;
Base::m_v1[i1] = n[i2];
Base::m_v1[i2] = -n[i1];
Base::m_v2[i0] = n[i1]*n[i1] + n[i2]*n[i2];
Base::m_v2[i1] = -n[i1]*n[i0];
Base::m_v2[i2] = -n[i2]*n[i0];
}
}
return res;
......@@ -82,7 +94,6 @@ ProjectedNormalCovarianceCurvature<DataPoint, _WFunctor, T>::init(const VectorTy
{
Base::init(_evalPos);
reset();
m_cog = Vector2::Zero();
m_cov = Mat22::Zero();
m_pass = FIRST_PASS;
......@@ -159,24 +170,24 @@ ProjectedNormalCovarianceCurvature<DataPoint, _WFunctor, T>::finalize ()
m_solver.computeDirect(m_cov);
m_k1 = m_solver.eigenvalues()(0);
m_k2 = m_solver.eigenvalues()(1);
Base::m_k1 = m_solver.eigenvalues()(0);
Base::m_k2 = m_solver.eigenvalues()(1);
// transform from local plane coordinates to world coordinates
m_v1 = m_tframe * m_solver.eigenvectors().col(0);
m_v2 = m_tframe * m_solver.eigenvectors().col(1);
Base::m_v1 = m_tframe * m_solver.eigenvectors().col(0);
Base::m_v2 = m_tframe * m_solver.eigenvectors().col(1);
//TODO(thib) which epsilon value should be chosen ?
// Scalar epsilon = Eigen::NumTraits<Scalar>::dummy_precision();
Scalar epsilon = Scalar(1e-3);
if(m_k1<epsilon && m_k2<epsilon)
if(Base::m_k1<epsilon && Base::m_k2<epsilon)
{
m_k1 = Scalar(0);
m_k2 = Scalar(0);
Base::m_k1 = Scalar(0);
Base::m_k2 = Scalar(0);
// set principal directions from fitted plane
m_v2 = m_tframe.col(0);
m_v1 = m_tframe.col(1);
Base::m_v2 = m_tframe.col(0);
Base::m_v1 = m_tframe.col(1);
}
return STABLE;
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment