Commit abff3e73 authored by Thibault Lejemble's avatar Thibault Lejemble

progress on mlsSphereFitDer methods

TODO:

- some cleaning

- compute d2w

- use nume, deno, dNume, dDeno from base class

- handle spaceDer/scaleDer in dPotential() and dNormal()
parent ece1dcce
......@@ -5,15 +5,9 @@
#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
{
......@@ -27,6 +21,19 @@ protected:
PROVIDES_NORMAL_SPACE_DERIVATIVE
};
enum
{
Dim = DataPoint::Dim //!< Dimension of the ambient space
};
//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)
enum
{
DerDim = DER_NB_DERIVATIVES(Base::Type,Dim) //!< Number of dimensions used for the differentiation
};
public:
typedef typename Base::Scalar Scalar; /*!< \brief Inherited scalar type*/
typedef typename Base::VectorType VectorType; /*!< \brief Inherited vector type*/
......@@ -34,19 +41,32 @@ public:
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;
/*!
\brief Static squared matrix of scalars with a size adapted to the
differentiation type
The size is adapted to the differentiation type (scale and/or space)
only. Such Matrix is used to represent Hessian matrix of multivariate
single-valued functions such as \f$ u_c \f$ and \f$ u_q \f$ of an
AlgebraicSphere obtained from a fitting procedure.
*/
typedef Eigen::Matrix< Scalar, DerDim, DerDim > Matrix;
/*!
\brief Static matrix of scalars with a size adapted to the
differentiation type and the dimension of the ambiant space.
typedef Eigen::Matrix< Scalar,
DER_NB_DERIVATIVES(Base::Type,DataPoint::Dim),
DataPoint::Dim*DER_NB_DERIVATIVES(Base::Type,DataPoint::Dim)> MatrixArray;
The size is adapted to the differentiation type (scale and/or space) and
the dimension of the ambiant space. Such Matrix is used to represent
Hessian matrix of multivariate multi-valued functions such as
\f$ u_l \f$ of an AlgebraicSphere obtained from a fitting procedure.
//TODO(thib) add explanation on how to use this matrix (using block...)
*/
typedef Eigen::Matrix< Scalar, DerDim, Dim*DerDim > MatrixArray;
protected:
// computation data
......@@ -54,6 +74,9 @@ protected:
m_d2SumDotPP, /*!< \brief Sum of the squared relative positions with twice differenciated weights */
m_d2SumW; /*!< \brief Sum of queries weight with twice differenciated weights */
MatrixArray m_d2SumP,
m_d2SumN;
public:
// results
Matrix m_d2Uc, /*!< \brief Second derivative of the hyper-sphere constant term */
......
......@@ -32,8 +32,13 @@ MlsSphereFitDer<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
m_d2SumDotPN += d2w * _nei.normal().dot(q);
m_d2SumDotPP += d2w * q.squaredNorm();
m_d2SumW += d2w;
}
for(int i=0; i<Dim; ++i)
{
m_d2SumP.template block<DerDim,DerDim>(0,i*Dim) += d2w * q[i];
m_d2SumN.template block<DerDim,DerDim>(0,i*Dim) += d2w * _nei.normal()[i];
}
}
return bResult;
}
......@@ -66,10 +71,27 @@ MlsSphereFitDer<DataPoint, _WFunctor, T>::finalize()
- 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 sumdSumPdSumN = Matrix::Zero(); //TODO(thib) is this the transpose of eachother ?
Matrix sumdSumNdSumP = Matrix::Zero(); //TODO(thib) is this the transpose of eachother ?
Matrix sumd2SumPdSumN = Matrix::Zero();
Matrix sumd2SumNdSumP = Matrix::Zero();
Matrix sumdSumPdSumP = Matrix::Zero();
Matrix sumd2SumPdSumP = Matrix::Zero();
for(int i=0; i<Dim; ++i)
{
sumdSumPdSumN += Base::m_dSumN.row(i).transpose()*Base::m_dSumP.row(i); //TODO(thib) is this the transpose of eachother ?
sumdSumNdSumP += Base::m_dSumP.row(i).transpose()*Base::m_dSumN.row(i); //TODO(thib) is this the transpose of eachother ?
sumd2SumPdSumN += m_d2SumP.template block<DerDim,DerDim>(0,i*Dim)*Base::m_sumN(i);
sumd2SumNdSumP += m_d2SumN.template block<DerDim,DerDim>(0,i*Dim)*Base::m_sumP(i);
sumdSumPdSumP += Base::m_dSumP.row(i).transpose()*Base::m_dSumP.row(i); //TODO(thib) simplification ?
sumd2SumPdSumP += m_d2SumP.template block<DerDim,DerDim>(0,i*Dim)*Base::m_sumP(i);
}
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_sumW*Base::m_sumW*( Base::m_sumW*(sumdSumPdSumN+sumdSumNdSumP+sumd2SumPdSumN+sumd2SumNdSumP)
+ 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)
......@@ -77,8 +99,8 @@ MlsSphereFitDer<DataPoint, _WFunctor, T>::finalize()
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_sumW*Base::m_sumW*( Scalar(2.)*Base::m_sumW*(sumdSumPdSumP+sumd2SumPdSumP)
+ 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
......@@ -86,19 +108,49 @@ MlsSphereFitDer<DataPoint, _WFunctor, T>::finalize()
Scalar deno2 = deno*deno;
m_d2Uq = Scalar(.5)/(deno2*deno2)*(deno2*(dDeno.transpose()*dNume + deno*d2Nume
- dNume.transpose()*dDeno - nume*d2Deno)
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
for(int i=0; i<Dim; ++i)
{
m_d2Ul.template block<DerDim,DerDim>(0,i*Dim) = invSumW*(
m_d2SumN.template block<DerDim,DerDim>(0,i*Dim)
- Scalar(2.)*( m_d2Uq*Base::m_sumP[i]
+ Base::m_dSumP.row(i).transpose()*Base::m_dUq
+ Base::m_uq*m_d2SumP.template block<DerDim,DerDim>(0,i*Dim)
+ Base::m_dUq.transpose()*Base::m_dSumP.row(i))
- Base::m_ul[i]*m_d2SumW
- Base::m_dUl.row(i).transpose()*Base::m_dSumW
- Base::m_dSumW.transpose()*Base::m_dUl.row(i));
}
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
Matrix sumdUldSumP = Matrix::Zero();
Matrix sumUld2SumP = Matrix::Zero();
Matrix sumd2UlsumP = Matrix::Zero();
Matrix sumdSumPdUl = Matrix::Zero();
for(int i=0; i<Dim; ++i)
{
sumdUldSumP += Base::m_dUl.row(i).transpose()*Base::m_dSumP.row(i);
sumUld2SumP += Base::m_ul[i]*m_d2SumP.template block<DerDim,DerDim>(0,i*Dim);
sumd2UlsumP += m_d2Ul.template block<DerDim,DerDim>(0,i*Dim)*Base::m_sumP[i];
sumdSumPdUl += Base::m_dSumP.row(i).transpose()*Base::m_dUl.row(i);
}
m_d2Uc = -invSumW*(
sumdUldSumP
+ sumUld2SumP
+ sumd2UlsumP
+ sumdSumPdUl
+ 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);
}
......
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