Commit 2017b32c authored by Thibault Lejemble's avatar Thibault Lejemble

tmp

parent d8ff5cf6
......@@ -4,11 +4,17 @@
file, You can obtain one at http://mozilla.org/MPL/2.0/.
*/
template < class DataPoint, class _WFunctor, typename T>
void
MlsSphereFitDer<DataPoint, _WFunctor, T>::init(const VectorType& _evalPos)
{
// Some operations in addNeighbor() cannot compile in Scale Differentiation
// mode only because Eigen find matrix size problem, even though these
// operations theoretically happen when Space Differentiation mode is used.
// TODO: find more flexible accessors, or use Constexpr If (c++17)
static_assert(Base::isSpaceDer(),"This extension cannot be used with scale \
differentiation only due to compile time checking errors.");
Base::init(_evalPos);
m_d2Uc = Matrix::Zero(),
......@@ -35,13 +41,21 @@ MlsSphereFitDer<DataPoint, _WFunctor, T>::addNeighbor(const DataPoint& _nei)
VectorType q = _nei.pos() - Base::basisCenter();
// compute weight derivatives
Scalar w = Base::m_w.w(q, _nei);
ScalarArray dw = ScalarArray::Zero();
Matrix d2w = Matrix::Zero();
if (Base::isScaleDer())
{
dw[0] = Base::m_w.scaledw(q, _nei);
d2w(0,0) = Base::m_w.scaled2w(q, _nei);
}
if (Base::isSpaceDer())
{
dw.template tail<Dim>() = Base::m_w.spacedw(q, _nei).transpose();
d2w.template bottomRightCorner<Dim,Dim>() = Base::m_w.spaced2w(q, _nei);
}
if (Base::isScaleDer() && Base::isSpaceDer())
{
......@@ -52,12 +66,41 @@ 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*DerDim) += d2w * q[i];
m_d2SumN.template block<DerDim,DerDim>(0,i*DerDim) += d2w * _nei.normal()[i];
}
if (Base::isSpaceDer())
{
m_d2SumDotPP.template bottomRightCorner<Dim,Dim>() -= Scalar(4.) * dw.template tail<Dim>().transpose() * q.transpose();
m_d2SumDotPP.template bottomRightCorner<Dim,Dim>().diagonal().array() += Scalar(2.) * w;
m_d2SumDotPN.template bottomRightCorner<Dim,Dim>() -= dw.template tail<Dim>().transpose() * _nei.normal().transpose() +
_nei.normal() * dw.template tail<Dim>();
for(int i=0; i<Dim; ++i)
{
m_d2SumP.template block<DerDim,DerDim>(0,i*DerDim).template bottomRightCorner<Dim,Dim>().col(i) -= dw.template tail<Dim>().transpose();
m_d2SumP.template block<DerDim,DerDim>(0,i*DerDim).template bottomRightCorner<Dim,Dim>().row(i) -= dw.template tail<Dim>();
}
}
if (Base::isSpaceDer() && Base::isScaleDer())
{
m_d2SumDotPP.template bottomLeftCorner<Dim,1>() -= Scalar(2.) * dw[0] * q;
m_d2SumDotPP.template topRightCorner<1,Dim>() = m_d2SumDotPP.template bottomLeftCorner<Dim,1>().transpose();
m_d2SumDotPN.template bottomLeftCorner<Dim,1>() -= dw[0] * _nei.normal();
m_d2SumDotPN.template topRightCorner<1,Dim>() = m_d2SumDotPN.template bottomLeftCorner<Dim,1>().transpose();
for(int i=0; i<Dim; ++i)
{
m_d2SumP.template block<DerDim,DerDim>(0,i*DerDim).template topRightCorner<1,Dim>()[i] -= dw[0];
m_d2SumP.template block<DerDim,DerDim>(0,i*DerDim).template bottomLeftCorner<Dim,1>() =
m_d2SumP.template block<DerDim,DerDim>(0,i*DerDim).template topRightCorner<1,Dim>();
}
}
}
return bResult;
}
......
......@@ -154,7 +154,7 @@ OrientedSphereDer<DataPoint, _WFunctor, T, Type>::addNeighbor(const DataPoint &
dw[0] = Base::m_w.scaledw(q, _nei);
if (Type & FitSpaceDer)
dw.template tail<int(DataPoint::Dim)>() = -Base::m_w.spacedw(q, _nei).transpose();
dw.template tail<int(DataPoint::Dim)>() = Base::m_w.spacedw(q, _nei).transpose();
// increment
m_dSumW += dw;
......@@ -163,6 +163,15 @@ OrientedSphereDer<DataPoint, _WFunctor, T, Type>::addNeighbor(const DataPoint &
m_dSumDotPN += dw * _nei.normal().dot(q);
m_dSumDotPP += dw * q.squaredNorm();
if (isSpaceDer())
{
Scalar w = Base::m_w.w(q, _nei);
m_dSumDotPN.template tail<DataPoint::Dim>() -= w * _nei.normal();
m_dSumDotPP.template tail<DataPoint::Dim>() -= Scalar(2.) * w * q;
m_dSumP.template rightCols<DataPoint::Dim>().diagonal().array() -= w;
}
return true;
}
......
......@@ -22,7 +22,7 @@ DistWeightFunc<DataPoint, WeightKernel>::spacedw( const VectorType& _q,
{
VectorType result = VectorType::Zero();
Scalar d = _q.norm();
if (d <= m_t && d != Scalar(0.)) result = (_q / (d * m_t)) * m_wk.df(d/m_t);
if (d <= m_t && d != Scalar(0.)) result = -(_q / (d * m_t)) * m_wk.df(d/m_t);
return result;
}
......@@ -71,7 +71,7 @@ DistWeightFunc<DataPoint, WeightKernel>::scaleSpaced2w( const VectorType& _q,
VectorType result = VectorType::Zero();
Scalar d = _q.norm();
if (d <= m_t && d != Scalar(0.))
result = -_q/(m_t*m_t)*(m_wk.df(d/m_t)/d + m_wk.ddf(d/m_t)/m_t);
result = _q/(m_t*m_t)*(m_wk.df(d/m_t)/d + m_wk.ddf(d/m_t)/m_t);
return result;
}
......
d_x(p-x) = +Id ... not -Id ???
=> ça voudrait dire que y'a des signes - presque partout ?
DistWeightFunc need const DataPoint& /*attributes*/
\ No newline at end of file
......@@ -26,3 +26,5 @@ add_multi_test(dnormal_sphere.cpp)
add_multi_test(scale_space_consistency.cpp)
add_multi_test(smooth_weight_kernel.cpp)
add_multi_test(dist_weight_func.cpp)
add_multi_test(mls_fit_der.cpp)
add_multi_test(testTmp.cpp)
/*
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 test/Grenaille/dnormal_paraboloid.cpp
\brief Test validity of dnormal
*/
#include "../common/testing.h"
#include "../common/testUtils.h"
#include <vector>
using namespace std;
using namespace Grenaille;
#include <fstream>
template<typename VectorArrayT>
void exportOBJ(const std::string& filename, const VectorArrayT& data)
{
std::ofstream file(filename);
if(!file.is_open())
{
std::cout << "Failed to open file [" << filename << "]" << std::endl;
return;
}
for(const auto& pt : data)
{
file << "v " << pt.pos().transpose() << "\n";
}
file.close();
std::cout << "### Data exported to [" << filename << "]" << std::endl;
}
template<typename DataPoint, typename Fit, typename WeightFunc> //, typename Fit, typename WeightFunction>
void testFunction(bool _bAddNoise = false)
{
// Define related structure
typedef typename DataPoint::Scalar Scalar;
typedef typename DataPoint::VectorType VectorType;
typedef typename DataPoint::MatrixType MatrixType;
typedef typename DataPoint::QuaternionType QuaternionType;
//generate sampled paraboloid
int nbPoints = Eigen::internal::random<int>(1000, 10000);
// VectorType coef = 100 * VectorType(Eigen::internal::random<Scalar>(-1,1), Eigen::internal::random<Scalar>(-1,1), 0);
// create saddle surface
VectorType coef = VectorType::Zero();
coef[0] = Eigen::internal::random<Scalar>(1,100);
coef[1] = -coef[0];
Scalar analysisScale = Eigen::internal::random<Scalar>(0.5, 1.0)/coef[0];
Scalar epsilon = 0.1; //testEpsilon<Scalar>();
if( _bAddNoise ) // relax a bit the testing threshold
epsilon = Scalar(0.002*MAX_NOISE);
vector<DataPoint> vectorPoints(nbPoints);
QuaternionType q;
for(unsigned int i = 0; i < vectorPoints.size(); ++i)
{
vectorPoints[i] = getPointOnParaboloid<DataPoint>(VectorType::Zero(),
coef,
q,
analysisScale*Scalar(1.2),
_bAddNoise);
}
// Test at central point if dNormal is correct
DataPoint point;
point.pos() = VectorType::Zero();
point.normal() = VectorType::Zero();
point.normal().z() = Scalar(0.1);
{
Fit fit;
fit.setWeightFunc(WeightFunc(analysisScale));
for(int i=0; i< 10 ; ++i)
{
fit.init(point.pos());
fit.compute(vectorPoints.cbegin(), vectorPoints.cend());
if(fit.isStable())
point.pos() = fit.project(point.pos());
}
if(fit.isStable())
{
MatrixType dN = fit.dNormal().template rightCols<DataPoint::Dim>();
Eigen::SelfAdjointEigenSolver<MatrixType> solver(dN);
VectorType eigValues = solver.eigenvalues();
MatrixType eigVectors = solver.eigenvectors();
if( ! (std::abs(eigValues[1]) < epsilon &&
eigValues[0] < Scalar(0.) &&
eigValues[2] > Scalar(0.) &&
std::abs(eigValues[0]+eigValues[2])/std::abs(2*eigValues[0]) < epsilon ))
{
cout << eigValues << endl;
exportOBJ("what.obj", vectorPoints);
}
VERIFY( std::abs(eigValues[1]) < epsilon );
VERIFY( eigValues[0] < Scalar(0.) );
VERIFY( eigValues[2] > Scalar(0.) );
VERIFY( std::abs(eigValues[0]+eigValues[2])/std::abs(2*eigValues[0]) < epsilon );
VERIFY( std::abs(Scalar(1.)-std::abs(VectorType(0,1,0).dot(eigVectors.col(0)))) < epsilon );
VERIFY( std::abs(Scalar(1.)-std::abs(VectorType(0,0,1).dot(eigVectors.col(1)))) < epsilon );
VERIFY( std::abs(Scalar(1.)-std::abs(VectorType(1,0,0).dot(eigVectors.col(2)))) < epsilon );
}
else
{
cout << "Warning: unstable fit" << endl;
}
}
}
template<typename Scalar, int Dim>
void callSubTests()
{
typedef PointPositionNormal<Scalar, Dim> Point;
typedef DistWeightFunc<Point, SmoothWeightKernel<Scalar> > WeightSmoothFunc;
typedef DistWeightFunc<Point, ConstantWeightKernel<Scalar> > WeightConstantFunc;
typedef Basket<Point,WeightSmoothFunc,OrientedSphereFit, OrientedSphereSpaceDer> FitSmoothOrientedSpaceDer;
typedef Basket<Point,WeightConstantFunc,OrientedSphereFit, OrientedSphereSpaceDer> FitConstantOrientedSpaceDer;
typedef Basket<Point,WeightSmoothFunc,OrientedSphereFit, OrientedSphereScaleSpaceDer> FitSmoothOrientedScaleSpaceDer;
typedef Basket<Point,WeightConstantFunc,OrientedSphereFit, OrientedSphereScaleSpaceDer> FitConstantOrientedScaleSpaceDer;
typedef Basket<Point,WeightSmoothFunc,OrientedSphereFit, OrientedSphereSpaceDer, MlsSphereFitDer> FitSmoothOrientedSpaceMlsDer;
typedef Basket<Point,WeightConstantFunc,OrientedSphereFit, OrientedSphereSpaceDer, MlsSphereFitDer> FitConstantOrientedSpaceMlsDer;
typedef Basket<Point,WeightSmoothFunc,OrientedSphereFit, OrientedSphereScaleSpaceDer, MlsSphereFitDer> FitSmoothOrientedScaleSpaceMlsDer;
typedef Basket<Point,WeightConstantFunc,OrientedSphereFit, OrientedSphereScaleSpaceDer, MlsSphereFitDer> FitConstantOrientedScaleSpaceMlsDer;
cout << "Testing with perfect plane (oriented / unoriented)..." << endl;
for(int i = 0; i < g_repeat; ++i)
{
CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedSpaceDer, WeightSmoothFunc>() ));
CALL_SUBTEST(( testFunction<Point, FitConstantOrientedSpaceDer, WeightConstantFunc>() ));
CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedScaleSpaceDer, WeightSmoothFunc>() ));
CALL_SUBTEST(( testFunction<Point, FitConstantOrientedScaleSpaceDer, WeightConstantFunc>() ));
// CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedSpaceMlsDer, WeightSmoothFunc>() ));
// CALL_SUBTEST(( testFunction<Point, FitConstantOrientedSpaceMlsDer, WeightConstantFunc>() ));
// CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedScaleSpaceMlsDer, WeightSmoothFunc>() ));
// CALL_SUBTEST(( testFunction<Point, FitConstantOrientedScaleSpaceMlsDer, WeightConstantFunc>() ));
}
cout << "Ok!" << endl;
}
int main(int argc, char** argv)
{
if(!init_testing(argc, argv))
{
return EXIT_FAILURE;
}
cout << "Test dnormal" << endl;
// callSubTests<float, 3>();
callSubTests<long double, 3>();
// callSubTests<double, 3>();
}
This diff is collapsed.
This diff is collapsed.
/*
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 test/Grenaille/mls_paraboloid_der.cpp
\brief Test validity of dNormal
*/
#include "../common/testing.h"
#include "../common/testUtils.h"
#include <vector>
using namespace std;
using namespace Grenaille;
template<typename DataPoint, typename Fit, typename WeightFunc> //, typename Fit, typename WeightFunction>
void testFunction(bool _bAddNoise = false)
{
// Define related structure
typedef typename DataPoint::Scalar Scalar;
typedef typename DataPoint::VectorType VectorType;
typedef typename DataPoint::QuaternionType QuaternionType;
//generate sampled paraboloid
int nbPoints = Eigen::internal::random<int>(100, 1000);
VectorType coef = 100 * VectorType(Eigen::internal::random<Scalar>(-1,1), Eigen::internal::random<Scalar>(-1,1), 0);
coef = VectorType(2, 0, 0);
// Scalar analysisScale = Scalar(.001);
Scalar analysisScale = Scalar(2.);
Scalar epsilon = testEpsilon<Scalar>();
if( _bAddNoise ) // relax a bit the testing threshold
epsilon = Scalar(0.002*MAX_NOISE);
vector<DataPoint> vectorPoints(nbPoints);
QuaternionType q;
for(unsigned int i = 0; i < vectorPoints.size(); ++i)
{
vectorPoints[i] = getPointOnParaboloid<DataPoint>(VectorType::Zero(),
coef,
q,
analysisScale*Scalar(1.2));
}
// Test at central point if the principal curvature directions are correct
DataPoint point;
point.pos() = VectorType::Zero();
point.normal() = VectorType::Zero();
point.normal().z() = Scalar(1.);
{
Fit fit;
fit.setWeightFunc(WeightFunc(analysisScale));
for(int i=0; i<5; ++i)
{
fit.init(point.pos());
fit.compute(vectorPoints.cbegin(), vectorPoints.cend());
point.pos() = fit.project(point.pos());
}
if(fit.isStable())
{
VectorType k1v = fit.k1Direction();
VectorType k2v = fit.k2Direction();
VectorType normal = fit.normal();
Scalar k1 = fit.k1();
Scalar k2 = fit.k2();
VectorType k1Theoritical = VectorType(1,0,0);
VectorType k2Theoritical = VectorType(0,1,0);
VectorType normalTheoritical = VectorType(0,0,1);
Scalar absDot1 = std::abs(k1Theoritical.dot(k1v));
Scalar absDot2 = std::abs(k2Theoritical.dot(k2v));
Scalar absDotN = std::abs(normalTheoritical.dot(normal));
std::cout << "k1="<< k1 << " ## " << k1v.transpose() << " ## err1=" << std::abs(Scalar(1.)-absDot1) << std::endl;
std::cout << "k2="<< k2 << " ## " << k2v.transpose() << " ## err2=" << std::abs(Scalar(1.)-absDot2) << std::endl;
std::cout << "N =" << normal.transpose() << " ## errN=" << std::abs(Scalar(1.)-absDotN) << std::endl << std::endl;
auto tmp = fit.dNormal();
std::cout << "dNormal:\n" << tmp << std::endl << std::endl;
//TODO(thib) to be continued...
}
}
}
template<typename Scalar, int Dim>
void callSubTests()
{
typedef PointPositionNormal<Scalar, Dim> Point;
typedef DistWeightFunc<Point, SmoothWeightKernel<Scalar> > WeightSmoothFunc;
typedef DistWeightFunc<Point, ConstantWeightKernel<Scalar> > WeightConstantFunc;
//TODO(thib) delete this
typedef Basket<Point,WeightSmoothFunc, OrientedSphereFit, OrientedSphereSpaceDer/*, MlsSphereFitDer*/, CurvatureEstimator> FitTest;
typedef Basket<Point,WeightSmoothFunc,OrientedSphereFit, OrientedSphereSpaceDer, MlsSphereFitDer, CurvatureEstimator> FitSmoothOrientedSpaceMlsDer;
typedef Basket<Point,WeightConstantFunc,OrientedSphereFit, OrientedSphereSpaceDer, MlsSphereFitDer, CurvatureEstimator> FitConstantOrientedSpaceMlsDer;
typedef Basket<Point,WeightSmoothFunc,OrientedSphereFit, OrientedSphereScaleSpaceDer, MlsSphereFitDer, CurvatureEstimator> FitSmoothOrientedScaleSpaceMlsDer;
typedef Basket<Point,WeightConstantFunc,OrientedSphereFit, OrientedSphereScaleSpaceDer, MlsSphereFitDer, CurvatureEstimator> FitConstantOrientedScaleSpaceMlsDer;
cout << "Testing with perfect primitive (oriented / unoriented)..." << endl;
for(int i = 0; i < g_repeat; ++i)
{
//TODO(thib) delete this
// CALL_SUBTEST(( testFunction<Point, FitTest, WeightSmoothFunc>() ));
CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedSpaceMlsDer, WeightSmoothFunc>() ));
// CALL_SUBTEST(( testFunction<Point, FitConstantOrientedSpaceMlsDer, WeightConstantFunc>() ));
// CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedScaleSpaceMlsDer, WeightSmoothFunc>() ));
// CALL_SUBTEST(( testFunction<Point, FitConstantOrientedScaleSpaceMlsDer, WeightConstantFunc>() ));
}
cout << "Ok!" << endl;
cout << "Testing with noise on position and normals (oriented / unoriented)..." << endl;
for(int i = 0; i < g_repeat; ++i)
{
// CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedSpaceMlsDer, WeightSmoothFunc>(true) ));
// CALL_SUBTEST(( testFunction<Point, FitConstantOrientedSpaceMlsDer, WeightConstantFunc>(true) ));
// CALL_SUBTEST(( testFunction<Point, FitSmoothOrientedScaleSpaceMlsDer, WeightSmoothFunc>(true) ));
// CALL_SUBTEST(( testFunction<Point, FitConstantOrientedScaleSpaceMlsDer, WeightConstantFunc>(true) ));
}
cout << "Ok!" << endl;
}
int main(int argc, char** argv)
{
if(!init_testing(argc, argv))
{
return EXIT_FAILURE;
}
cout << "Test normal estimation" << endl;
//callSubTests<double, 2>();
// callSubTests<float, 3>();
callSubTests<long double, 3>();
// callSubTests<double, 3>();
}
/*
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/.
*/
//TODO(thib) This is just a temporary test to compile mlsSphereFitDer
/*!
\file test/Grenaille/gls_sphere_der.cpp
\brief Test validity of GLS derivatives for a sphere
*/
#include "../common/testing.h"
#include "../common/testUtils.h"
#include <vector>
using namespace std;
using namespace Grenaille;
template<typename DataPoint, typename Fit, typename WeightFunc> //, typename Fit, typename WeightFunction>
void testFunction(bool _bAddPositionNoise = false, bool _bAddNormalNoise = false)
{
// Define related structure
typedef typename DataPoint::Scalar Scalar;
typedef typename DataPoint::VectorType VectorType;
typedef typename Fit::ScalarArray ScalarArray;
//generate sampled sphere
int nbPoints = Eigen::internal::random<int>(100, 1000);
Scalar radius = Eigen::internal::random<Scalar>(1,10);
VectorType center = VectorType::Random() * Eigen::internal::random<Scalar>(1, 10000);
Scalar analysisScale = Scalar(10.) * std::sqrt(Scalar(4. * M_PI) * radius * radius / nbPoints);
Scalar epsilon = testEpsilon<Scalar>();
//Scalar radisuEpsilon = epsilon * radius;
vector<DataPoint> vectorPoints(nbPoints);
for(unsigned int i = 0; i < vectorPoints.size(); ++i)
{
vectorPoints[i] = getPointOnSphere<DataPoint>(radius, center, _bAddPositionNoise, _bAddNormalNoise);
}
// Test for each point if the Derivatives of kappa are equal to 0
#pragma omp parallel for
for(int i = 0; i < int(vectorPoints.size()); ++i)
{
Fit fit;
fit.setWeightFunc(WeightFunc(analysisScale));
fit.init(vectorPoints[i].pos());
fit.compute(vectorPoints.cbegin(), vectorPoints.cend());
if(fit.isStable())
{
ScalarArray dkappa = fit.dkappa();
for(int i = 0; i < dkappa.size(); ++i)
{
VERIFY( Eigen::internal::isMuchSmallerThan(dkappa[i], Scalar(1.), epsilon) );
}
}
}
}
template<typename Scalar, int Dim>
void callSubTests()
{
typedef PointPositionNormal<Scalar, Dim> Point;
typedef DistWeightFunc<Point, SmoothWeightKernel<Scalar> > WeightSmoothFunc;
typedef Basket<Point, WeightSmoothFunc, OrientedSphereFit, GLSParam, OrientedSphereScaleSpaceDer, MlsSphereFitDer, GLSDer> FitSmoothOriented;
cout << "Testing with perfect sphere (oriented / unoriented)..." << endl;
for(int i = 0; i < g_repeat; ++i)
{
CALL_SUBTEST(( testFunction<Point, FitSmoothOriented, WeightSmoothFunc>() ));
}
cout << "Ok!" << endl;
/*cout << "Testing with noise on position and normals (oriented / unoriented)..." << endl;
for(int i = 0; i < g_repeat; ++i)
{
CALL_SUBTEST(( testFunction<Point, FitSmoothOriented, WeightSmoothFunc>(true, true) ));
}
cout << "Ok!" << endl;*/
}
int main(int argc, char** argv)
{
if(!init_testing(argc, argv))
{
return EXIT_FAILURE;
}
cout << "Test sphere derivatives with GLSParam and OrientedSphereFit..." << endl;
// callSubTests<double, 2>(); //TODO(thib) uncomment this
callSubTests<float, 3>();
callSubTests<long double, 3>();
callSubTests<double, 3>();
}
......@@ -31,14 +31,17 @@ void checkConsistency(const vector<DataPoint> vectorPoints, typename DataPoint::
const int nbPoints = vectorPoints.size();
// Note: FitScaleDer with only scale differentiation is not tested here since
// it does not compile with the MlsSphereFitDer extension...
// test that dPotential and dNormal wrt scale, space and both are the same regardless of the derivation type
#pragma omp parallel for
for(int i=0; i<nbPoints; ++i)
{
FitScaleDer fitScaleDer;
fitScaleDer.setWeightFunc(WeightFunc(analysisScale));
fitScaleDer.init(vectorPoints[i].pos());
fitScaleDer.compute(vectorPoints.cbegin(), vectorPoints.cend());
// FitScaleDer fitScaleDer;
// fitScaleDer.setWeightFunc(WeightFunc(analysisScale));
// fitScaleDer.init(vectorPoints[i].pos());
// fitScaleDer.compute(vectorPoints.cbegin(), vectorPoints.cend());
FitSpaceDer fitSpaceDer;
fitSpaceDer.setWeightFunc(WeightFunc(analysisScale));
......@@ -50,19 +53,19 @@ void checkConsistency(const vector<DataPoint> vectorPoints, typename DataPoint::
fitScaleSpaceDer.init(vectorPoints[i].pos());
fitScaleSpaceDer.compute(vectorPoints.cbegin(), vectorPoints.cend());
VERIFY( fitScaleDer.isStable()==fitSpaceDer.isStable() && fitScaleDer.isStable()==fitScaleSpaceDer.isStable() );
VERIFY( /*fitScaleDer.isStable()==fitSpaceDer.isStable() &&*/ fitSpaceDer.isStable()==fitScaleSpaceDer.isStable() );
if(fitScaleDer.isStable() && fitSpaceDer.isStable() && fitScaleSpaceDer.isStable())
if(/*fitScaleDer.isStable() &&*/ fitSpaceDer.isStable() && fitScaleSpaceDer.isStable())
{
typename FitScaleDer::ScalarArray dPotentialScaleDer = fitScaleDer.dPotential();
// typename FitScaleDer::ScalarArray dPotentialScaleDer = fitScaleDer.dPotential();
typename FitSpaceDer::ScalarArray dPotentialSpaceDer = fitSpaceDer.dPotential();
typename FitScaleSpaceDer::ScalarArray dPotentialScaleSpaceDer = fitScaleSpaceDer.dPotential();
typename FitSpaceDer::VectorArray dNormalSpaceDer = fitSpaceDer.dNormal();
typename FitScaleSpaceDer::VectorArray dNormalScaleSpaceDer = fitScaleSpaceDer.dNormal();
Scalar dt_potential1 = dPotentialScaleDer[0];
Scalar dt_potential2 = dPotentialScaleSpaceDer[0];
// Scalar dt_potential1 = dPotentialScaleDer[0];
// Scalar dt_potential2 = dPotentialScaleSpaceDer[0];
VectorType dx_potential1 = dPotentialSpaceDer.template tail<DataPoint::Dim>();
VectorType dx_potential2 = dPotentialScaleSpaceDer.template tail<DataPoint::Dim>();
......@@ -70,7 +73,7 @@ void checkConsistency(const vector<DataPoint> vectorPoints, typename DataPoint::
MatrixType dx_normal1 = dNormalSpaceDer.template rightCols<DataPoint::Dim>();
MatrixType dx_normal2 = dNormalScaleSpaceDer.template rightCols<DataPoint::Dim>();
VERIFY( std::abs(dt_potential1-dt_potential2) < epsilon );
// VERIFY( std::abs(dt_potential1-dt_potential2) < epsilon );
VERIFY( (dx_potential1-dx_potential2).norm() < epsilon );
for(int i=0; i<DataPoint::Dim; ++i)
......
#include "../common/testing.h"
#include "../common/testUtils.h"
#include <unsupported/Eigen/AutoDiff>
using namespace std;
using namespace Grenaille;
template<typename ScalarT, class VectorT, class MatrixT>
struct DistanceFunction
{
static ScalarT dist(const VectorT& q, ScalarT t) {return q.norm()/t;}
static VectorT dx_dist(const VectorT& q, ScalarT t) {return -q/(t*q.norm());}
static ScalarT dt_dist(const VectorT& q, ScalarT t) {return -q.norm()/(t*t);}
static MatrixT d2x_dist(const VectorT& q, ScalarT t) {return ScalarT(1.)/(t*q.norm())*( MatrixT::Identity() - q*q.transpose()/(q.norm()*q.norm()) );}
static ScalarT d2t_dist(const VectorT& q, ScalarT t) {return ScalarT(2.)*q.norm()/(t*t*t);}
static VectorT d2tx_dist(const VectorT& q, ScalarT t) {return q/(t*t*q.norm());}
};
template<typename DataPoint, typename WeightKernel>
void testFunctionAutoDiff()
{
// Define related structure
typedef typename WeightKernel::Scalar ScalarDiff;
typedef typename ScalarDiff::Scalar Scalar;
typedef typename DataPoint::VectorType VectorTypeDiff;
typedef typename DataPoint::MatrixType MatrixTypeDiff;
Scalar epsilon = testEpsilon<Scalar>();
Scalar scale = Eigen::internal::random<Scalar>(0.10, 10.0);
ScalarDiff t = ScalarDiff(scale, DataPoint::Dim+1, 0);
// random center
VectorTypeDiff center = Eigen::internal::random<Scalar>(1, 10000)*VectorTypeDiff::Random();
VectorTypeDiff x;
for(int i=0; i<DataPoint::Dim; ++i) {
x[i] = ScalarDiff(center[i].value(), DataPoint::Dim+1, i+1);
}
// rejection sampling inside a sphere of radius equal to scale
VectorTypeDiff p = scale*VectorTypeDiff::Random();
while ( p.norm() < 0.001*scale || 0.999*scale < p.norm())
p = scale*VectorTypeDiff::Random();
p += center;
// centered query
VectorTypeDiff q = p-x