Commit 87112240 authored by Thibault Lejemble's avatar Thibault Lejemble

improve curvature estimation

add optional parameter to compute the tangent plane either from dN, or from N
parent 4ab8cc13
......@@ -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;
}
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