Commit ece1dcce authored by Thibault Lejemble's avatar Thibault Lejemble

add new mlsSphereFitDer class

not yet implemented
parent 5d6e7bdf
#ifndef _GRENAILLE_MLS_SPHERE_FIT_DER_
#define _GRENAILLE_MLS_SPHERE_FIT_DER_
namespace Grenaille
{
//TODO(thib) get Type from *SphereDer class
//TODO(thib) or should it be put it in internal:: namespace ?
// and create other public class as OrientedSphereDer does
template < class DataPoint, class _WFunctor, typename T>
class MlsSphereFitDer : public T
{
private:
typedef T Base;
protected:
enum
{
Check = Base::PROVIDES_ALGEBRAIC_SPHERE_DERIVATIVE,
PROVIDES_NORMAL_SPACE_DERIVATIVE
};
public:
typedef typename Base::Scalar Scalar; /*!< \brief Inherited scalar type*/
typedef typename Base::VectorType VectorType; /*!< \brief Inherited vector type*/
typedef typename Base::VectorArray VectorArray; /*!< \brief Inherited vector array type */
typedef typename Base::ScalarArray ScalarArray; /*!< \brief Inherited scalar array type */
//TODO(thib) redundant with OrientedSphereDer macro, which cannot be used here due to internal namespace...
#define DER_NB_DERIVATIVES(TYPE,DIM) ((TYPE & internal::FitScaleDer) ? 1 : 0 ) + ((TYPE & internal::FitSpaceDer) ? DIM : 0)
//TODO(thib) use appropriate storage order
//#define GLS_DER_STORAGE_ORDER(TYPE) ((TYPE & FitSpaceDer) ? Eigen::RowMajor : Eigen::ColMajor )
typedef Eigen::Matrix< Scalar,
DER_NB_DERIVATIVES(Base::Type,DataPoint::Dim),
DER_NB_DERIVATIVES(Base::Type,DataPoint::Dim)> Matrix;
typedef Eigen::Matrix< Scalar,
DER_NB_DERIVATIVES(Base::Type,DataPoint::Dim),
DataPoint::Dim*DER_NB_DERIVATIVES(Base::Type,DataPoint::Dim)> MatrixArray;
protected:
// computation data
Matrix m_d2SumDotPN, /*!< \brief Sum of the dot product betwen relative positions and normals with twice differenciated weights */
m_d2SumDotPP, /*!< \brief Sum of the squared relative positions with twice differenciated weights */
m_d2SumW; /*!< \brief Sum of queries weight with twice differenciated weights */
public:
// results
Matrix m_d2Uc, /*!< \brief Second derivative of the hyper-sphere constant term */
m_d2Uq; /*!< \brief Second derivative of the hyper-sphere quadratic term */
MatrixArray m_d2Ul; /*!< \brief Second derivative of the hyper-sphere linear term */
public:
/************************************************************************/
/* Initialization */
/************************************************************************/
/*! \see Concept::FittingProcedureConcept::init() */
MULTIARCH void init(const VectorType &evalPos);
/************************************************************************/
/* Processing */
/************************************************************************/
/*! \see Concept::FittingProcedureConcept::addNeighbor() */
MULTIARCH bool addNeighbor(const DataPoint &nei);
/*! \see Concept::FittingProcedureConcept::finalize() */
MULTIARCH FIT_RESULT finalize();
/**************************************************************************/
/* Use results */
/**************************************************************************/
/*! \brief Returns the derivatives of the scalar field at the evaluation point */
MULTIARCH inline ScalarArray dPotential() const;
/*! \brief Returns the second derivatives of the scalar field at the evaluation point */
MULTIARCH inline VectorArray dNormal() const;
}; //class MlsSphereFitDer
#include "mlsSphereFitDer.hpp"
} //namespace Grenaille
#endif
template < class DataPoint, class _WFunctor, typename T>
void
MlsSphereFitDer<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
Base::init(_evalPos);
m_d2Uc = Matrix::Zero(),
m_d2Uq = Matrix::Zero();
m_d2Ul = MatrixArray::Zero();
m_d2SumDotPN = Matrix::Zero();
m_d2SumDotPP = Matrix::Zero();
m_d2SumW = Matrix::Zero();
}
template < class DataPoint, class _WFunctor, typename T>
bool
MlsSphereFitDer<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
{
bool bResult = Base::addNeighbor(_nei);
if(bResult)
{
// centered basis
VectorType q = _nei.pos() - Base::basisCenter();
Matrix d2w;
d2w;//TODO(thib) computes d2w wrt Scale/Space der type
m_d2SumDotPN += d2w * _nei.normal().dot(q);
m_d2SumDotPP += d2w * q.squaredNorm();
m_d2SumW += d2w;
}
return bResult;
}
//TODO(thib) delete this macro
template<typename T> struct ____The_Unkown_Type_Is____;
#define GLUE1(X,Y) X##Y
#define GLUE(X,Y) GLUE1(X,Y)
#define WhatIsTheTypeOf(expr) ____The_Unkown_Type_Is____<decltype(expr)> GLUE(_,__LINE__)
template < class DataPoint, class _WFunctor, typename T>
FIT_RESULT
MlsSphereFitDer<DataPoint, _WFunctor, T>::finalize()
{
Base::finalize();
if (this->isReady())
{
//TODO(thib) use nume, deno, dNume, dDeno from base class to avoid calculating them twice
Scalar invSumW = Scalar(1.)/Base::m_sumW;
Scalar nume = Base::m_sumDotPN - invSumW*Base::m_sumP.dot(Base::m_sumN);
Scalar deno = Base::m_sumDotPP - invSumW*Base::m_sumP.dot(Base::m_sumP);
ScalarArray dNume = Base::m_dSumDotPN
- invSumW*invSumW * ( Base::m_sumW * ( Base::m_sumN.transpose() * Base::m_dSumP + Base::m_sumP.transpose() * Base::m_dSumN )
- Base::m_dSumW*Base::m_sumP.dot(Base::m_sumN) );
ScalarArray dDeno = Base::m_dSumDotPP
- invSumW*invSumW*( Scalar(2.) * Base::m_sumW * Base::m_sumP.transpose() * Base::m_dSumP
- Base::m_dSumW*Base::m_sumP.dot(Base::m_sumP) );
Matrix d2Nume = m_d2SumDotPN
- invSumW*invSumW*invSumW*invSumW*(
Base::m_sumW*Base::m_sumW*(/*Base::m_sumW*(0)*/ //TODO(thib) replace this
/*+ */Base::m_dSumW.transpose()*(Base::m_sumN.transpose()*Base::m_dSumP + Base::m_sumP.transpose()*Base::m_dSumN)
- (Base::m_sumP.transpose()*Base::m_sumN)*m_d2SumW.transpose()
- (Base::m_dSumN.transpose()*Base::m_sumP + Base::m_dSumP.transpose()*Base::m_sumN)*Base::m_dSumW)
- Scalar(2.)*Base::m_sumW*Base::m_dSumW.transpose()*(Base::m_sumW*(Base::m_sumN.transpose()*Base::m_dSumP + Base::m_sumP.transpose()*Base::m_dSumN)
- (Base::m_sumP.transpose()*Base::m_sumN)*Base::m_dSumW));
Matrix d2Deno = m_d2SumDotPP
- invSumW*invSumW*invSumW*invSumW*(
Base::m_sumW*Base::m_sumW*(/*Scalar(2.)*Base::m_sumW*(0)*/ //TODO(thib) replace this
/*+ */Scalar(2.)*Base::m_dSumW.transpose()*(Base::m_sumP.transpose()*Base::m_dSumP)
- (Base::m_sumP.transpose()*Base::m_sumP)*m_d2SumW.transpose()
- Scalar(2.)*(Base::m_dSumP.transpose()*Base::m_sumP)*Base::m_dSumW)
- Scalar(2.)*Base::m_sumW*Base::m_dSumW.transpose()*(Scalar(2.)*Base::m_sumW*Base::m_sumP.transpose()*Base::m_dSumP
- (Base::m_sumP.transpose()*Base::m_sumP)*Base::m_dSumW));
Scalar deno2 = deno*deno;
m_d2Uq = Scalar(.5)/(deno2*deno2)*(deno2*(dDeno.transpose()*dNume + deno*d2Nume
- dNume.transpose()*dDeno - nume*d2Deno)
- Scalar(2.)*deno*dDeno.transpose()*(deno*dNume - nume*dDeno));
m_d2Ul; //TODO(thib) compute this
m_d2Uc = -invSumW*(/*0*/ // TODO(thib) replace this
//+ 0 // TODO(thib) replace this
//+ 0 // TODO(thib) replace this
//+ 0 // TODO(thib) replace this
/*+ */Base::m_dUq.transpose()*Base::m_dSumDotPP + Base::m_uq*m_d2SumDotPP
+ Base::m_dSumDotPP.transpose()*Base::m_dUq + m_d2Uq*Base::m_sumDotPP
+ Base::m_uc*m_d2SumW + Base::m_dUc.transpose()*Base::m_dSumW
- Base::m_dSumW.transpose()*Base::m_dUc);
}
return Base::m_eCurrentState;
}
template < class DataPoint, class _WFunctor, typename T>
typename MlsSphereFitDer<DataPoint, _WFunctor, T>::ScalarArray
MlsSphereFitDer<DataPoint, _WFunctor, T>::dPotential() const
{
//TODO(thib) handle spaceDer/scaleDer
return Base::m_dUc + Base::m_ul;
}
template < class DataPoint, class _WFunctor, typename T>
typename MlsSphereFitDer<DataPoint, _WFunctor, T>::VectorArray
MlsSphereFitDer<DataPoint, _WFunctor, T>::dNormal() const
{
//TODO(thib) handle spaceDer/scaleDer
//TODO(thib) this is just the hessian for now, need to normalize by the norm and then differenciate
return m_d2Uc + Base::m_dUl + Base::m_dUl.transpose() + Scalar(2.)*Base::m_uq*VectorArray::Ones();
}
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