Skip to content

Commit b15562e

Browse files
nmelladonmellado
committed
Add wrapper to Ponca v1.2
1 parent 39e34a8 commit b15562e

3 files changed

Lines changed: 261 additions & 0 deletions

File tree

‎src/wrapper-versions/PoncaV1x2.sh‎

Lines changed: 9 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,9 @@
1+
set -o errexit
2+
pushd external/ponca
3+
git checkout v1.2
4+
popd
5+
cp wrapper-versions/PoncaV1x2/* .
6+
rm -f *o
7+
R -e "Rcpp::compileAttributes('..')"
8+
R CMD INSTALL ..
9+
R --vanilla < ../tests/testthat/test-CRAN.R

‎src/wrapper-versions/PoncaV1x2/curvatureEstimation.cpp‎

Lines changed: 204 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,204 @@
1+
#include "curvatureEstimation.h"
2+
3+
#include "external/ponca/Ponca/Ponca"
4+
#include <iostream>
5+
6+
#define DIMENSION 3
7+
8+
using namespace Ponca;
9+
10+
template<typename InMType>
11+
class MyPointMap
12+
{
13+
public:
14+
enum {Dim = DIMENSION};
15+
using Scalar = double;
16+
typedef Eigen::Matrix<Scalar, Dim, 1> VectorType;
17+
typedef Eigen::Matrix<Scalar, Dim, Dim> MatrixType;
18+
typedef Eigen::VectorBlock<InMType> InnerVectorType;
19+
20+
PONCA_MULTIARCH inline MyPointMap(InMType &mat, int _pId)
21+
: m_pos (mat.col(_pId).head(3)),
22+
m_normal (mat.col(_pId).tail(3))
23+
{}
24+
25+
PONCA_MULTIARCH inline const InnerVectorType& pos() const { return m_pos; }
26+
PONCA_MULTIARCH inline const InnerVectorType& normal() const { return m_normal; }
27+
28+
public:
29+
InnerVectorType m_pos, m_normal;
30+
// VectorType m_pos, m_normal;
31+
};
32+
33+
// struct PassThroughConverter{
34+
// inline void operator()( const std::vector< MyPointMap > &&i, std::vector< MyPointMap > & o ) {
35+
// o = std::move(i);
36+
// }
37+
// };
38+
39+
class MyPointSimple
40+
{
41+
public:
42+
enum {Dim = DIMENSION};
43+
using Scalar = double;
44+
typedef Eigen::Matrix<Scalar, Dim, 1> VectorType;
45+
typedef Eigen::Matrix<Scalar, Dim, Dim> MatrixType;
46+
47+
PONCA_MULTIARCH inline MyPointSimple(const VectorType& p, const VectorType& n)
48+
: m_pos (p), m_normal(n)
49+
{}
50+
PONCA_MULTIARCH inline const VectorType& pos() const { return m_pos; }
51+
PONCA_MULTIARCH inline const VectorType& normal() const { return m_normal; }
52+
53+
private:
54+
VectorType m_pos, m_normal;
55+
};
56+
57+
58+
59+
/// Generate acceleration structure
60+
Ponca::KdTree<MyPointSimple> tree;
61+
62+
#define MIN_NOISE 0.99
63+
#define MAX_NOISE 1.01
64+
/*! \brief Generate points on a plane */
65+
template <typename DataPoint>
66+
[[nodiscard]] DataPoint getPointOnPlane(const typename DataPoint::VectorType& _vPosition,
67+
const typename DataPoint::Scalar& _width,
68+
const typename DataPoint::Scalar& _height,
69+
const typename DataPoint::VectorType& _localxAxis,
70+
const typename DataPoint::VectorType& _localyAxis,
71+
const bool _bAddPositionNoise = true)
72+
{
73+
using Scalar = typename DataPoint::Scalar;
74+
using VectorType = typename DataPoint::VectorType;
75+
76+
const Scalar u = Eigen::internal::random<Scalar>(-_width / Scalar(2), _width / Scalar(2));
77+
const Scalar v = Eigen::internal::random<Scalar>(-_height / Scalar(2), _height / Scalar(2));
78+
79+
VectorType vRandomPosition = _vPosition + u * _localxAxis + v * _localyAxis;
80+
81+
if (_bAddPositionNoise)
82+
{
83+
vRandomPosition = vRandomPosition +
84+
VectorType::Random().normalized() * Eigen::internal::random<Scalar>(0., 1. - MIN_NOISE);
85+
}
86+
87+
return DataPoint(vRandomPosition, _localxAxis.cross(_localyAxis));
88+
}
89+
90+
void generatePointClouds(Eigen::MatrixXd& points,
91+
Eigen::MatrixXd& queries,
92+
double dataScale)
93+
{
94+
MyPointSimple::VectorType position = MyPointSimple::VectorType::Random();
95+
96+
for (int i = 0; i != points.rows(); ++i)
97+
{
98+
auto p = getPointOnPlane<MyPointSimple>(position, dataScale,dataScale, {1,0,0}, {0,1,0});
99+
points.row(i) << p.pos().x(), p.pos().y(),p.pos().z(),p.normal().x(),p.normal().y(),p.normal().z();
100+
}
101+
102+
for (int i = 0; i != queries.rows(); ++i)
103+
{
104+
auto p = getPointOnPlane<MyPointSimple>(position, dataScale,dataScale, {1,0,0}, {0,1,0});
105+
queries.row(i) << p.pos().x(), p.pos().y(),p.pos().z();
106+
}
107+
// reset KdTree
108+
tree.build(std::vector<MyPointSimple>());
109+
}
110+
111+
bool buildKdTree(const Eigen::MatrixXd& points)
112+
{
113+
int nPoints = points.rows();
114+
115+
/// Bind dataset to Ponca representation
116+
std::vector<MyPointSimple> data;
117+
data.reserve(nPoints);
118+
for (int i = 0; i != nPoints; ++i)
119+
{
120+
data.emplace_back(points.row(i).head(3),points.row(i).tail(3));
121+
}
122+
tree.build(data);
123+
124+
return nPoints != 0;
125+
}
126+
127+
struct ComputeReturnType
128+
{
129+
// number of fits
130+
int nbFit{0};
131+
// mean number of neighbors
132+
int kNeiMean{0};
133+
};
134+
135+
136+
template <typename Fit, bool range, typename Param>
137+
ComputeReturnType computeFit(const Eigen::MatrixXd& queries, Param p)
138+
{
139+
ComputeReturnType ret;
140+
141+
if (tree.point_count() == 0)
142+
{
143+
std::cerr<< "KdTree has not been initialized" << std::endl;
144+
return ret;
145+
}
146+
147+
using W = typename Fit::WeightFunction;
148+
using Point = typename Fit::DataPoint;
149+
using Vector = typename Point::VectorType;
150+
using Scalar = typename Point::Scalar;
151+
152+
int nQueries = queries.rows();
153+
154+
// compute queries
155+
for (int i = 0; i != nQueries; ++i)
156+
{
157+
Vector q(queries.row(i).head(3));
158+
Fit f;
159+
f.setWeightFunc(W(Scalar(p)));
160+
f.init(q);
161+
if (range)
162+
f.computeWithIds(tree.range_neighbors(q, p), tree.point_data());
163+
else
164+
f.computeWithIds(tree.k_nearest_neighbors(q, p), tree.point_data());
165+
if (f.isStable())
166+
{
167+
ret.nbFit++;
168+
ret.kNeiMean += f.getNumNeighbors();
169+
}
170+
}
171+
ret.kNeiMean /= ret.nbFit;
172+
173+
return ret;
174+
}
175+
176+
using NF = DistWeightFunc<MyPointSimple, SmoothWeightKernel<double> > ;
177+
using ASOBasket = Basket<MyPointSimple, NF, OrientedSphereFit>;
178+
using ASOFit = BasketDiff<ASOBasket, FitSpaceDer, OrientedSphereDer, MlsSphereFitDer>;
179+
using PlaneFit = Ponca::Basket<MyPointSimple, NF, CovariancePlaneFit>;
180+
181+
int asoCurvatureEstimation(const Eigen::MatrixXd& queries, double scale, int& meanNeiSize)
182+
{
183+
auto ret = computeFit<ASOFit, true>(queries, scale);
184+
meanNeiSize = ret.kNeiMean;
185+
return ret.nbFit;
186+
}
187+
188+
int planeFit(const Eigen::MatrixXd& queries, double scale, int& meanNeiSize)
189+
{
190+
auto ret = computeFit<PlaneFit, true>(queries, scale);
191+
meanNeiSize = ret.kNeiMean;
192+
return ret.nbFit;
193+
}
194+
int asoCurvatureEstimation(const Eigen::MatrixXd& queries, int k)
195+
{
196+
auto ret = computeFit<ASOFit, false>(queries, k);
197+
return ret.nbFit;
198+
}
199+
200+
int planeFit(const Eigen::MatrixXd& queries, int k)
201+
{
202+
auto ret = computeFit<PlaneFit, false>(queries, k);
203+
return ret.nbFit;
204+
}

‎src/wrapper-versions/PoncaV1x2/curvatureEstimation.h‎

Lines changed: 48 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,48 @@
1+
#pragma once
2+
3+
#include <Eigen/Dense>
4+
5+
6+
/// \brief Default function use to generate the point cloud and the queries for the tests
7+
void generatePointClouds(Eigen::MatrixXd& points,
8+
Eigen::MatrixXd& queries,
9+
double dataScale);
10+
11+
/// \brief Compute the KdTree from an input set of points.
12+
///
13+
/// \note The KdTree is made accessible to the estimation functions through a global variable.
14+
///
15+
/// \note This function can be used either to measure the performances of the KdTree construction,
16+
/// or to prepare the KdTree for the estimation functions below
17+
///
18+
/// \return false if points is empty
19+
bool buildKdTree(const Eigen::MatrixXd& points);
20+
21+
22+
/*********************************************************************************************/
23+
/* Estimation functions */
24+
/* (require a preliminary call to generateData and then buildKdTree) */
25+
/* Several estimation functions can be called sequentially without recomputing the KdTree */
26+
/*********************************************************************************************/
27+
/// \brief Run principal curvature estimation using Algebraic Shape operator on range queries
28+
/// \warning requires a Kdtree, \see buildKdTree
29+
///
30+
/// \param meanNeiSize mean number of points in the neighborhood
31+
/// \return the number of points correctly analyzed, or -1 if the KdTree has not been initialized first
32+
int asoCurvatureEstimation(const Eigen::MatrixXd& queries, double scale, int &meanNeiSize);
33+
/// \brief Run covariance plane fitting on range queries
34+
/// \warning requires a Kdtree, \see buildKdTree
35+
///
36+
/// \param meanNeiSize mean number of points in the neighborhood
37+
/// \return the number of points correctly analyzed, or -1 if the KdTree has not been initialized first
38+
int planeFit(const Eigen::MatrixXd& queries, double scale, int &meanNeiSize);
39+
/// \brief Run principal curvature estimation using Algebraic Shape operator on k-neighbors
40+
/// \warning requires a Kdtree, \see buildKdTree
41+
///
42+
/// \return the number of points correctly analyzed, or -1 if the KdTree has not been initialized first
43+
int asoCurvatureEstimation(const Eigen::MatrixXd& queries, int k);
44+
/// \brief Run covariance plane fitting on k-neighbors
45+
/// \warning requires a Kdtree, \see buildKdTree
46+
///
47+
/// \return the number of points correctly analyzed, or -1 if the KdTree has not been initialized first
48+
int planeFit(const Eigen::MatrixXd& queries, int k);

0 commit comments

Comments
 (0)