Commit 8522a56d authored by Nicolas Mellado's avatar Nicolas Mellado

Add centeredBasis to CompactPlane

parent 56f27201
......@@ -64,8 +64,7 @@ public:
// computation data
Scalar m_sumW; /*!< \brief Sum of queries weight.*/
VectorType m_cog, /*!< \brief Gravity center of the neighborhood */
m_evalPos; /*!< \brief Center of the evaluation basis */
VectorType m_cog; /*!< \brief Gravity center of the neighborhood */
MatrixType m_cov; /*!< \brief Covariance matrix */
Solver m_solver; /*!<\brief Solver used to analyse the covariance matrix */
......@@ -122,12 +121,16 @@ public:
* \brief Express a point in ambiant space relatively to the tangent plane.
* Output vector is: [h, u, v]^T, where u, v are 2d coordinates on the plane,
* and h the height of the sample.
* \tparam ignoreTranslation must be set to true when passing vectors instead of points
*/
template <bool ignoreTranslation = false>
MULTIARCH inline VectorType worldToTangentPlane(const VectorType &_q) const;
/*!
* \brief Transform a point from the tangent plane [h, u, v]^T to ambiant space
* \tparam ignoreTranslation must be set to true when passing vectors instead of points
*/
template <bool ignoreTranslation = false>
MULTIARCH inline VectorType tangentPlaneToWorld(const VectorType &_q) const;
}; //class CovariancePlaneFit
......
......@@ -14,11 +14,10 @@ CovariancePlaneFit<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
// Setup primitive
Base::resetPrimitive();
// Base::basisCenter() = _evalPos;
Base::basisCenter() = _evalPos;
// Setup fitting internal values
m_sumW = Scalar(0.0);
m_evalPos = _evalPos;
m_cog = VectorType::Zero();
m_cov = MatrixType::Zero();
}
......@@ -27,7 +26,7 @@ template < class DataPoint, class _WFunctor, typename T>
bool
CovariancePlaneFit<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
{
VectorType q = _nei.pos() - m_evalPos;
VectorType q = _nei.pos() - Base::basisCenter();
// compute weight
Scalar w = m_w.w(q, _nei);
......@@ -57,25 +56,26 @@ CovariancePlaneFit<DataPoint, _WFunctor, T>::finalize ()
return Base::m_eCurrentState;
}
// Finalize the centroid (still expressed in local basis)
m_cog = m_cog/m_sumW;
// Center the covariance on the centroid:
m_cov -= m_cog * m_cog.transpose() / m_sumW;
m_cov += - Scalar(2) * m_cov/m_sumW + m_cog * m_cog.transpose();
// \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::m_eCurrentState = ( m_solver.info() == Eigen::Success ? STABLE : UNDEFINED );
Base::setPlane(m_solver.eigenvectors().col(0), m_cog);
Base::m_eCurrentState = ( m_solver.info() == Eigen::Success ? STABLE : UNDEFINED );
return Base::m_eCurrentState;
}
......@@ -88,22 +88,32 @@ CovariancePlaneFit<DataPoint, _WFunctor, T>::surfaceVariation () const
}
template < class DataPoint, class _WFunctor, typename T>
template <bool ignoreTranslation>
typename CovariancePlaneFit<DataPoint, _WFunctor, T>::VectorType
CovariancePlaneFit<DataPoint, _WFunctor, T>::worldToTangentPlane (const VectorType& _q) const
{
return m_solver.eigenvectors().transpose() * (_q - m_evalPos);
if (ignoreTranslation)
return m_solver.eigenvectors().transpose() * _q;
else {
// apply rotation and translation to get uv coordinates
return m_solver.eigenvectors().transpose() * (_q - Base::basisCenter());
}
}
template < class DataPoint, class _WFunctor, typename T>
template <bool ignoreTranslation>
typename CovariancePlaneFit<DataPoint, _WFunctor, T>::VectorType
CovariancePlaneFit<DataPoint, _WFunctor, T>::tangentPlaneToWorld (const VectorType& _lq) const
{
return m_solver.eigenvectors().transpose().inverse() * _lq + m_evalPos;
if (ignoreTranslation)
return m_solver.eigenvectors().transpose().inverse() * _lq;
else {
return m_solver.eigenvectors().transpose().inverse() * _lq + Base::basisCenter();
}
}
namespace internal
{
......@@ -133,7 +143,7 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::addNeighbor(const DataPoint
ScalarArray dw;
// centered basis
VectorType q = _nei.pos()-Base::m_evalPos;
VectorType q = _nei.pos()-Base::basisCenter();
// compute weight
if (Type & FitScaleDer)
......@@ -174,17 +184,10 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::finalize()
for(int k=0; k<NbDerivatives; ++k)
{
// Finalize the computation of dCov.
// Note that at this stage m_dCog = sum_i dw_i * (p_i-c).
// Since the covariance matrix is translation invariant,
// this step natuarally cancels the centered basis.
VectorType shifted_cog = Base::m_cog - Base::m_evalPos;
m_dCov[k] = m_dCov[k]
- 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;
- Base::m_cog * m_dCog.col(k).transpose()
- m_dCog.col(k) * Base::m_cog.transpose()
+ m_dSumW[k] * Base::m_cog * Base::m_cog.transpose();
// apply normalization by sumW:
m_dCog.col(k) = (m_dCog.col(k) - m_dSumW(k) * Base::m_cog) / Base::m_sumW;
......@@ -206,7 +209,10 @@ CovariancePlaneDer<DataPoint, _WFunctor, T, Type>::finalize()
VectorType dDiff = -m_dCog.col(k);
if(k>0 || !isScaleDer())
dDiff(isScaleDer() ? k-1 : k) += 1;
m_dDist(k) = m_dNormal.col(k).dot(Base::m_evalPos-Base::m_cog) + normal.dot(dDiff);
m_dDist(k) = m_dNormal.col(k).dot(Base::m_cog) + normal.dot(dDiff);
// \fixme we shouldn't need this normalization, however currently the derivatives are overestimated by a factor 2
m_dNormal /= Scalar(2.);
}
}
......
......@@ -52,8 +52,7 @@ public:
// computation data
Scalar m_sumW; /*!< \brief Sum of queries weight.*/
VectorType m_sumN, /*!< \brief Sum of the normal vectors */
m_sumP, /*!< \brief Sum of the relative positions */
m_evalPos; /*!< \brief Center of the evaluation basis */
m_sumP; /*!< \brief Sum of the relative positions */
WFunctor m_w; /*!< \brief Weight function (must inherits BaseWeightFunc) */
......
......@@ -13,7 +13,7 @@ MeanPlaneFit<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
// Setup primitive
Base::resetPrimitive();
m_evalPos = _evalPos;
Base::basisCenter() = _evalPos;
// Setup fitting internal values
m_sumP = VectorType::Zero();
......@@ -25,7 +25,7 @@ template < class DataPoint, class _WFunctor, typename T>
bool
MeanPlaneFit<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
{
VectorType q = _nei.pos() - m_evalPos;
VectorType q = _nei.pos() - Base::basisCenter();
// compute weight
Scalar w = m_w.w(q, _nei);
......@@ -55,7 +55,7 @@ MeanPlaneFit<DataPoint, _WFunctor, T>::finalize ()
return Base::m_eCurrentState;
}
Base::setPlane(m_sumN / m_sumW, m_sumP / m_sumW + m_evalPos);
Base::setPlane(m_sumN / m_sumW, m_sumP / m_sumW);
return Base::m_eCurrentState = STABLE;
}
......
......@@ -12,6 +12,7 @@
#include "primitive.h" // PrimitiveBase
#include <Eigen/Geometry>
#include <iostream>
namespace Grenaille
{
......@@ -65,12 +66,19 @@ public:
/*! \brief Weight Function */
typedef _WFunctor WFunctor;
private:
//! \brief Evaluation position (needed for centered basis)
VectorType m_p;
public:
/*! \brief Default constructor */
MULTIARCH inline CompactPlane()
: Base(), EigenBase()
{
m_p = VectorType::Zero();
resetPrimitive();
}
......@@ -79,7 +87,7 @@ public:
CompactPlane<DataPoint, WFunctor, T>& compactPlane()
{ return * static_cast<CompactPlane<DataPoint, WFunctor, T>*>(this); }
/*! \brief Set the scalar field values to 0 and reset the isNormalized()
/*! \brief Set the scalar field values to 0
status */
MULTIARCH inline void resetPrimitive()
{
......@@ -97,6 +105,11 @@ public:
return ! ((*this) == other);
}
/*! \brief Reading access to the basis center (evaluation position) */
MULTIARCH inline const VectorType& basisCenter () const { return m_p; }
/*! \brief Writing access to the (evaluation position) */
MULTIARCH inline VectorType& basisCenter () { return m_p; }
/* \brief Init the plane from a direction and a position
\param _dir Orientation of the plane, does not need to be normalized
\param _pos Position of the plane
......@@ -118,14 +131,16 @@ public:
MULTIARCH inline Scalar potential (const VectorType& _q) const
{
// The potential is the distance from the point to the plane
return EigenBase::signedDistance(_q);
return EigenBase::signedDistance(_q - m_p);
}
//! \brief Project a point on the plane
MULTIARCH inline VectorType project (const VectorType& _q) const
{
// Project on the normal vector and add the offset value
return EigenBase::projection(_q);
return EigenBase::projection(_q - m_p) + m_p;
}
//! \brief Scalar field gradient direction at the evaluation point
MULTIARCH inline VectorType primitiveGradient () const
{
......
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