Commit 0695db91 authored by Nicolas Mellado's avatar Nicolas Mellado

Add Sphere fit without normal

parent 3efe8402
/*
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/.
*/
#ifndef _GRENAILLE_SPHERE_FIT_
#define _GRENAILLE_SPHERE_FIT_
#include "algebraicSphere.h"
namespace Grenaille
{
/*!
\brief Algebraic Sphere fitting procedure on point sets without normals
Method published in \cite Guennebaud:2007:APSS.
\inherit Concept::FittingProcedureConcept
\see AlgebraicSphere
*/
template < class DataPoint, class _WFunctor, typename T = void >
class SphereFit : public AlgebraicSphere<DataPoint, _WFunctor>
{
private:
typedef AlgebraicSphere<DataPoint, _WFunctor> Base;
public:
/*! \brief Scalar type inherited from DataPoint*/
typedef typename Base::Scalar Scalar;
/*! \brief Vector type inherited from DataPoint*/
typedef typename Base::VectorType VectorType;
/*! \brief Weight Function*/
typedef _WFunctor WFunctor;
protected:
typedef Eigen::Matrix<Scalar, DataPoint::Dim+2, 1> VectorA;
typedef Eigen::Matrix<Scalar, DataPoint::Dim+2, DataPoint::Dim+2> MatrixA;
// computation data
MatrixA m_matA; /*!< \brief Covariance matrix of [1, p, p^2] */
Scalar m_sumW; /*!< \brief Sum of queries weight */
WFunctor m_w; /*!< \brief Weight function (must inherits BaseWeightFunc) */
public:
/*! \brief Default constructor */
MULTIARCH inline SphereFit()
: Base(){}
/**************************************************************************/
/* Initialization */
/**************************************************************************/
/*! \copydoc Concept::FittingProcedureConcept::setWeightFunc() */
MULTIARCH inline void setWeightFunc (const WFunctor& _w) { m_w = _w; }
/*! \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();
/**************************************************************************/
/* Results */
/**************************************************************************/
using Base::potential;
/*! \brief Value of the scalar field at the evaluation point */
MULTIARCH inline Scalar potential() const { return Base::m_uc; }
/*! \brief Value of the normal of the primitive at the evaluation point */
MULTIARCH inline VectorType normal() const { return Base::m_ul.normalized(); }
}; //class SphereFit
#include "sphereFit.hpp"
} //namespace Grenaille
#endif
/*
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/.
*/
template < class DataPoint, class _WFunctor, typename T>
void
SphereFit<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
// Setup primitive
Base::resetPrimitive();
Base::basisCenter() = _evalPos;
// Setup fitting internal values
m_sumW = Scalar(0.0);
m_matA.setZero();
}
template < class DataPoint, class _WFunctor, typename T>
bool
SphereFit<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
{
// centered basis
VectorType q = _nei.pos() - Base::basisCenter();
// compute weight
Scalar w = m_w.w(q, _nei);
if (w > Scalar(0.))
{
VectorA a;
a << 1, q, q.squaredNorm();
m_matA += w * a * a.transpose();
m_sumW += w;
++(Base::m_nbNeighbors);
return true;
}
return false;
}
template < class DataPoint, class _WFunctor, typename T>
FIT_RESULT
SphereFit<DataPoint, _WFunctor, T>::finalize ()
{
// handle specific configurations
// With less than 3 neighbors the fitting is undefined
if(m_sumW == Scalar(0.) || Base::m_nbNeighbors < 3)
{
Base::m_ul.setZero();
Base::m_uc = Scalar(0.);
Base::m_uq = Scalar(0.);
Base::m_isNormalized = false;
Base::m_eCurrentState = UNDEFINED;
return Base::m_eCurrentState;
}
MatrixA matC;
matC.setIdentity();
matC.template topRightCorner<1,1>() << -2;
matC.template bottomLeftCorner<1,1>() << -2;
matC.template topLeftCorner<1,1>() << 0;
matC.template bottomRightCorner<1,1>() << 0;
MatrixA invCpratt;
invCpratt.setIdentity();
invCpratt.template topRightCorner<1,1>() << -0.5;
invCpratt.template bottomLeftCorner<1,1>() << -0.5;
invCpratt.template topLeftCorner<1,1>() << 0;
invCpratt.template bottomRightCorner<1,1>() << 0;
MatrixA M = invCpratt * m_matA;
Eigen::EigenSolver<MatrixA> eig(M);
VectorA eivals = eig.eigenvalues().real();
int minId = -1;
for(int i=0 ; i<DataPoint::Dim+2 ; ++i)
{
Scalar ev = eivals(i);
if((ev>0) && (minId==-1 || ev<eivals(minId)))
minId = i;
}
//mLambda = eivals(minId);
VectorA vecU = eig.eigenvectors().col(minId).real();
Base::m_uq = vecU[1+DataPoint::Dim];
Base::m_ul = vecU.template segment<DataPoint::Dim>(1);
Base::m_uc = vecU[0];
Base::m_isNormalized = false;
if(Base::m_nbNeighbors < 6)
{
Base::m_eCurrentState = UNSTABLE;
}
else
{
Base::m_eCurrentState = STABLE;
}
return Base::m_eCurrentState;
}
/*
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/.
*/
......@@ -34,6 +34,7 @@
// not supported on cuda
#ifndef __CUDACC__
# include "Grenaille/Core/unorientedSphereFit.h"
# include "Grenaille/Core/sphereFit.h"
#endif
......
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