Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Several minor improvements #148

Open
wants to merge 6 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
17 changes: 14 additions & 3 deletions Ponca/src/Fitting/algebraicSphere.h
Original file line number Diff line number Diff line change
Expand Up @@ -58,8 +58,7 @@ class AlgebraicSphere : public T
//! \brief Is the implicit scalar field normalized using Pratt
bool m_isNormalized;

// results
public:
// results
Scalar m_uc {0}, /*!< \brief Constant parameter of the Algebraic hyper-sphere */
m_uq {0}; /*!< \brief Quadratic parameter of the Algebraic hyper-sphere */
VectorType m_ul {VectorType::Zero()}; /*!< \brief Linear parameter of the Algebraic hyper-sphere */
Expand Down Expand Up @@ -255,6 +254,18 @@ class AlgebraicSphere : public T
\warning The gradient is not normalized by default */
PONCA_MULTIARCH inline const VectorType& primitiveGradient () const { return m_ul; }

///\brief Mean curvature (inverse of the sphere radius). Return 0 for planes
PONCA_MULTIARCH inline Scalar meanCurvature () const { return Scalar(2)*m_uq; }

///\brief Read access to constant term \f$ u_c \f$
PONCA_MULTIARCH inline Scalar uc() const { return m_uc; }

///\brief Read access to linear term \f$ \mathbf{u}_l \f$
PONCA_MULTIARCH inline VectorType ul() const { return m_ul; }

///\brief Read access to quadratic term \f$ u_q \f$
PONCA_MULTIARCH inline Scalar uq() const { return m_uq; }

/*!
\brief Used to know if the fitting result to a plane
\return true if finalize() have been called and the fitting result to a plane
Expand All @@ -263,7 +274,7 @@ class AlgebraicSphere : public T
{
PONCA_MULTIARCH_STD_MATH(abs);
Scalar epsilon = Eigen::NumTraits<Scalar>::dummy_precision();
bool bPlanar = Eigen::internal::isMuchSmallerThan(abs(m_uq), Scalar(1.), epsilon);
bool bPlanar = Eigen::internal::isMuchSmallerThan(abs(meanCurvature()), Scalar(1.), epsilon);
bool bReady = Base::isReady();

return bReady && bPlanar;
Expand Down
42 changes: 21 additions & 21 deletions Ponca/src/SpatialPartitioning/KdTree/kdTree.h
Original file line number Diff line number Diff line change
Expand Up @@ -145,7 +145,7 @@ class KdTreeBase
/// \param points Input points
/// \param c Cast/Convert input point type to DataType
template<typename PointUserContainer, typename Converter>
inline void build(PointUserContainer&& points, Converter c);
inline void build(PointUserContainer&& points, Converter c, LeafSizeType min_cell_size = DEFAULT_LEAF_SIZE );

/// Convert a custom point container to the KdTree \ref PointContainer using \ref DataPoint default constructor
struct DefaultConverter
Expand All @@ -166,9 +166,9 @@ class KdTreeBase
/// \tparam PointUserContainer Input point container, transformed to PointContainer
/// \param points Input points
template<typename PointUserContainer>
inline void build(PointUserContainer&& points)
inline void build(PointUserContainer&& points, LeafSizeType min_cell_size = DEFAULT_LEAF_SIZE)
{
build(std::forward<PointUserContainer>(points), DefaultConverter());
build(std::forward<PointUserContainer>(points), DefaultConverter(), min_cell_size);
}

/// Clear tree data
Expand Down Expand Up @@ -224,13 +224,6 @@ class KdTreeBase
return m_min_cell_size;
}

/// Write leaf min size
inline void set_min_cell_size(LeafSizeType min_cell_size)
{
PONCA_DEBUG_ASSERT(min_cell_size > 0);
m_min_cell_size = min_cell_size;
}

// Index mapping -----------------------------------------------------------
public:
/// Return the point index associated with the specified sample index
Expand Down Expand Up @@ -289,7 +282,8 @@ public :

// Utilities ---------------------------------------------------------------
public:
inline bool valid() const;
/// \param ignore_duplicates By default, valid KdTree can have duplicated samples.
inline bool valid(bool ignore_duplicates = true) const;
inline void print(std::ostream& os, bool verbose = false) const;

// Data --------------------------------------------------------------------
Expand All @@ -298,7 +292,8 @@ public :
NodeContainer m_nodes;
IndexContainer m_indices;

LeafSizeType m_min_cell_size {64}; ///< Minimal number of points per leaf
static const LeafSizeType DEFAULT_LEAF_SIZE = 64;
LeafSizeType m_min_cell_size {DEFAULT_LEAF_SIZE}; ///< Minimal number of points per leaf
NodeIndexType m_leaf_count {0}; ///< Number of leaves in the Kdtree (computed during construction)

// Internal ----------------------------------------------------------------
Expand All @@ -315,7 +310,8 @@ public :
template<typename PointUserContainer, typename IndexUserContainer, typename Converter>
inline void buildWithSampling(PointUserContainer&& points,
IndexUserContainer sampling,
Converter c);
Converter c,
LeafSizeType min_cell_size = DEFAULT_LEAF_SIZE);

/// Generate a tree sampled from a custom contained type converted using a \ref KdTreeBase::DefaultConverter
/// \tparam PointUserContainer Input points, transformed to PointContainer
Expand All @@ -324,9 +320,10 @@ public :
/// \param sampling Samples used in the tree
template<typename PointUserContainer, typename IndexUserContainer>
inline void buildWithSampling(PointUserContainer&& points,
IndexUserContainer sampling)
IndexUserContainer sampling,
LeafSizeType min_cell_size = DEFAULT_LEAF_SIZE)
{
buildWithSampling(std::forward<PointUserContainer>(points), std::move(sampling), DefaultConverter());
buildWithSampling(std::forward<PointUserContainer>(points), std::move(sampling), DefaultConverter(), min_cell_size);
}

private:
Expand Down Expand Up @@ -361,10 +358,11 @@ class KdTreeDenseBase : public KdTreeBase<Traits>

/// Constructor generating a tree from a custom contained type converted using a \ref KdTreeBase::DefaultConverter
template<typename PointUserContainer>
inline explicit KdTreeDenseBase(PointUserContainer&& points)
inline explicit KdTreeDenseBase(PointUserContainer&& points,
typename Base::LeafSizeType min_cell_size = Base::DEFAULT_LEAF_SIZE)
: Base()
{
this->build(std::forward<PointUserContainer>(points));
this->build(std::forward<PointUserContainer>(points), min_cell_size);
}
};

Expand Down Expand Up @@ -397,10 +395,11 @@ class KdTreeSparseBase : public KdTreeBase<Traits>

/// Constructor generating a tree from a custom contained type converted using a \ref KdTreeBase::DefaultConverter
template<typename PointUserContainer>
inline explicit KdTreeSparseBase(PointUserContainer&& points)
inline explicit KdTreeSparseBase(PointUserContainer&& points,
typename Base::LeafSizeType min_cell_size = Base::DEFAULT_LEAF_SIZE)
: Base()
{
this->build(std::forward<PointUserContainer>(points));
this->build(std::forward<PointUserContainer>(points), min_cell_size);
}

/// Constructor generating a tree sampled from a custom contained type converted using a \ref KdTreeBase::DefaultConverter
Expand All @@ -409,10 +408,11 @@ class KdTreeSparseBase : public KdTreeBase<Traits>
/// \param point Input points
/// \param sampling Samples used in the tree
template<typename PointUserContainer, typename IndexUserContainer>
inline KdTreeSparseBase(PointUserContainer&& points, IndexUserContainer sampling)
inline KdTreeSparseBase(PointUserContainer&& points, IndexUserContainer sampling,
typename Base::LeafSizeType min_cell_size = Base::DEFAULT_LEAF_SIZE)
: Base()
{
this->buildWithSampling(std::forward<PointUserContainer>(points), std::move(sampling));
this->buildWithSampling(std::forward<PointUserContainer>(points), std::move(sampling), min_cell_size);
}

using Base::buildWithSampling;
Expand Down
18 changes: 13 additions & 5 deletions Ponca/src/SpatialPartitioning/KdTree/kdTree.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -8,11 +8,11 @@

template<typename Traits>
template<typename PointUserContainer, typename Converter>
inline void KdTreeBase<Traits>::build(PointUserContainer&& points, Converter c)
inline void KdTreeBase<Traits>::build(PointUserContainer&& points, Converter c, LeafSizeType min_cell_size)
{
IndexContainer ids(points.size());
std::iota(ids.begin(), ids.end(), 0);
this->buildWithSampling(std::forward<PointUserContainer>(points), std::move(ids), std::move(c));
this->buildWithSampling(std::forward<PointUserContainer>(points), std::move(ids), std::move(c), min_cell_size);
}

template<typename Traits>
Expand All @@ -25,11 +25,13 @@ void KdTreeBase<Traits>::clear()
}

template<typename Traits>
bool KdTreeBase<Traits>::valid() const
bool KdTreeBase<Traits>::valid(bool ignore_duplicates) const
{
// Check if point container is empty
if (m_points.empty())
return m_nodes.empty() && m_indices.empty();

// Check that at least one node has been created, and index container is not empty
if(m_nodes.empty() || m_indices.empty())
{
return false;
Expand All @@ -38,7 +40,7 @@ bool KdTreeBase<Traits>::valid() const
std::vector<bool> b(point_count(), false);
for(IndexType idx : m_indices)
{
if(idx < 0 || point_count() <= idx || b[idx])
if(idx < 0 || point_count() <= idx || ( ignore_duplicates ? false : b[idx]))
{
return false;
}
Expand Down Expand Up @@ -120,11 +122,15 @@ template<typename Traits>
template<typename PointUserContainer, typename IndexUserContainer, typename Converter>
inline void KdTreeBase<Traits>::buildWithSampling(PointUserContainer&& points,
IndexUserContainer sampling,
Converter c)
Converter c,
LeafSizeType min_cell_size)
{
PONCA_DEBUG_ASSERT(points.size() <= MAX_POINT_COUNT);
this->clear();

// Set cell size
this->m_min_cell_size = min_cell_size;

// Move, copy or convert input samples
c(std::forward<PointUserContainer>(points), m_points);

Expand All @@ -150,6 +156,8 @@ void KdTreeBase<Traits>::build_rec(NodeIndexType node_id, IndexType start, Index
node.set_is_leaf(
end-start <= m_min_cell_size ||
level >= Traits::MAX_DEPTH ||
// Stop descending if the content of the node is not one unique point
aabb.diagonal().squaredNorm() <= Eigen::NumTraits<Scalar>::epsilon() ||
Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@azaleostu I did this trick to stop the recursion when a node contains only several copies of the same point. Unfortunately it does not do the trick, and I cannot reproduce the segfault on my computer.
Could you please give it a try ?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I am getting consistent segfaults when there are duplicates. I haven't been able the find the source yet but it looks like duplicate samples can break the node structure.

For example in some runs I have been getting the same node index being processed multiple times, even at different depths:

Node 0, range: [0, 24] (depth: 1), samples: [7, 5, 5, 7, 6, 5, 6, 7, 8, 10, 10, 11, 0, 9, 4, 3, 2, 1, 0, 9, 9, 2, 0, 0]
Node 1, range: [0, 12] (depth: 2), samples: [11, 8, 6, 6, 7, 5, 5, 7, 5, 10, 10, 7]
Node 4, range: [4, 12] (depth: 3), samples: [5, 5, 5, 7, 7, 10, 10, 7]
Node 2, range: [12, 24] (depth: 2), samples: [0, 0, 0, 3, 2, 1, 0, 2, 9, 9, 4, 9]
Node 0, range: [12, 20] (depth: 3), samples: [1, 3, 0, 0, 2, 0, 0, 2]

There also seems to be some corrupted data after the build, like empty leaves, inconsistent child indices, negative child indices, etc.

KdTree:
  MaxNodes: 9223372036854775808
  MaxPoints: 8589934592
  MaxDepth: 32
  PointCount: 12
  SampleCount: 24
  NodeCount: 11
  Samples: [
    11, 8, 6, 6, 5, 5, 5, 7, 7, 10,
    10, 7, 1, 3, 0, 0, 2, 0, 0, 2,
    9, 9, 4, 9]
  Nodes:
    - Type: Inner
      SplitDim: 1
      SplitValue: -0.121189
      FirstChild: 9
    - Type: Leaf
      Start: 20
      Size: 4
    - Type: Inner
      SplitDim: 1
      SplitValue: 0.00595114
      FirstChild: 1701603584
    - Type: Leaf
      Start: 0
      Size: 4
    - Type: Inner
      SplitDim: 1
      SplitValue: 0.664357
      FirstChild: 2002541616
    - Type: Leaf
      Start: 4
      Size: 3
    - Type: Leaf
      Start: 7
      Size: 5
    - Type: Leaf
      Start: 0
      Size: 0
    - Type: Leaf
      Start: 0
      Size: 0
    - Type: Leaf
      Start: 12
      Size: 2
    - Type: Leaf
      Start: 14
      Size: 6

I'm still not sure where the issue might be but I think there could be some edge cases in the build algorithm when it comes to duplicates, I will try again when I can.

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, i ran to similar inconsistent behaviors.
Something I found: the reallocation of the nodeContainer breaks (when the number of nodes goes larger than the preallocation) the structure: nodes copy fails, and nodes are inconsistently initialized as leaf or internal, intervals are lost, ...

What I don't get is: how is this related to duplicates ?!

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm also unsure since it seems like the algorithm should work with duplicates, I see that the node array alloc is based on point_count() instead of sample_count(), maybe changing this could at least fix the realloc issues?

Copy link
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The fix could work indeed.
But it does not solve the reallocation problem, we need to fix this..

// Since we add 2 nodes per inner node we need to stop if we can't add
// them both
(NodeIndexType)m_nodes.size() > MAX_NODE_COUNT - 2);
Expand Down
52 changes: 32 additions & 20 deletions tests/common/kdtree_utils.h
Original file line number Diff line number Diff line change
Expand Up @@ -51,9 +51,11 @@ class MyPoint {
};

template<typename Scalar, typename VectorContainer>
bool check_range_neighbors(const VectorContainer& points, const std::vector<int>& sampling, int index, Scalar r, const std::vector<int>& neighbors)
bool check_range_neighbors(const VectorContainer& points, const std::vector<int>& sampling, int index, Scalar r,
const std::vector<int>& neighbors,
bool allow_duplicates = false)
{
if (has_duplicate(neighbors))
if (!allow_duplicates && has_duplicate(neighbors))
{
return false;
}
Expand Down Expand Up @@ -103,9 +105,11 @@ bool check_range_neighbors(const VectorContainer& points, const std::vector<int>
}

template<typename Scalar, typename VectorType, typename VectorContainer>
bool check_range_neighbors(const VectorContainer& points, const std::vector<int>& sampling, const VectorType& point, Scalar r, const std::vector<int>& neighbors)
bool check_range_neighbors(const VectorContainer& points, const std::vector<int>& sampling,
const VectorType& point, Scalar r, const std::vector<int>& neighbors,
bool allow_duplicates = false)
{
if (has_duplicate(neighbors))
if (!allow_duplicates && has_duplicate(neighbors))
{
return false;
}
Expand Down Expand Up @@ -147,14 +151,15 @@ bool check_range_neighbors(const VectorContainer& points, const std::vector<int>
}

template<typename Scalar, typename VectorContainer>
bool check_k_nearest_neighbors(const VectorContainer& points, int index, int k, const std::vector<int>& neighbors)
bool check_k_nearest_neighbors(const VectorContainer& points, int index, int k,
const std::vector<int>& neighbors, bool allow_duplicates = false)
{
if (int(points.size()) > k && int(neighbors.size()) != k)
{
return false;
}

if (has_duplicate(neighbors))
if (!allow_duplicates && has_duplicate(neighbors))
{
return false;
}
Expand Down Expand Up @@ -190,14 +195,15 @@ bool check_k_nearest_neighbors(const VectorContainer& points, int index, int k,
}

template<typename Scalar, typename VectorContainer>
bool check_k_nearest_neighbors(const VectorContainer& points, const std::vector<int>& sampling, int index, int k, const std::vector<int>& neighbors)
bool check_k_nearest_neighbors(const VectorContainer& points, const std::vector<int>& sampling, int index, int k,
const std::vector<int>& neighbors, bool allow_duplicates = false)
{
if (int(points.size()) > k && int(neighbors.size()) != k)
{
return false;
}

if (has_duplicate(neighbors))
if (!allow_duplicates && has_duplicate(neighbors))
{
return false;
}
Expand Down Expand Up @@ -242,14 +248,16 @@ bool check_k_nearest_neighbors(const VectorContainer& points, const std::vector<
}

template<typename Scalar, typename VectorType, typename VectorContainer>
bool check_k_nearest_neighbors(const VectorContainer& points, const std::vector<int>& sampling, const VectorType& point, int k, const std::vector<int>& neighbors)
bool check_k_nearest_neighbors(const VectorContainer& points, const std::vector<int>& sampling,
const VectorType& point, int k, const std::vector<int>& neighbors,
bool allow_duplicates = false)
{
if (int(sampling.size()) >= k && int(neighbors.size()) != k)
{
return false;
}

if (has_duplicate(neighbors))
if (!allow_duplicates && has_duplicate(neighbors))
{
return false;
}
Expand Down Expand Up @@ -287,14 +295,15 @@ bool check_k_nearest_neighbors(const VectorContainer& points, const std::vector<


template<typename Scalar, typename VectorType, typename VectorContainer>
bool check_k_nearest_neighbors(const VectorContainer& points, const VectorType& point, int k, const std::vector<int>& neighbors)
bool check_k_nearest_neighbors(const VectorContainer& points, const VectorType& point, int k,
const std::vector<int>& neighbors, bool allow_duplicates = false)
{
if (int(points.size()) >= k && int(neighbors.size()) != k)
{
return false;
}

if (has_duplicate(neighbors))
if (!allow_duplicates && has_duplicate(neighbors))
{
return false;
}
Expand Down Expand Up @@ -323,24 +332,27 @@ bool check_k_nearest_neighbors(const VectorContainer& points, const VectorType&


template<typename Scalar, typename VectorContainer>
bool check_nearest_neighbor(const VectorContainer& points, int index, int nearest)
bool check_nearest_neighbor(const VectorContainer& points, int index, int nearest, bool allow_duplicates = false)
{
return check_k_nearest_neighbors<Scalar, VectorContainer>(points, index, 1, { nearest });
return check_k_nearest_neighbors<Scalar, VectorContainer>(points, index, 1, { nearest }, allow_duplicates);
}
template<typename Scalar, typename VectorType, typename VectorContainer>
bool check_nearest_neighbor(const VectorContainer& points, const VectorType& point, int nearest)
bool check_nearest_neighbor(const VectorContainer& points, const VectorType& point, int nearest,
bool allow_duplicates = false)
{
return check_k_nearest_neighbors<Scalar, VectorType, VectorContainer>(points, point, 1, { nearest });
return check_k_nearest_neighbors<Scalar, VectorType, VectorContainer>(points, point, 1, { nearest }, allow_duplicates);
}

template<typename Scalar, typename VectorType, typename VectorContainer>
bool check_nearest_neighbor(const VectorContainer& points, const std::vector<int>& sampling, const VectorType& point, int nearest)
bool check_nearest_neighbor(const VectorContainer& points, const std::vector<int>& sampling,
const VectorType& point, int nearest, bool allow_duplicates = false)
{
return check_k_nearest_neighbors<Scalar, VectorContainer>(points, sampling, point, 1, { nearest });
return check_k_nearest_neighbors<Scalar, VectorContainer>(points, sampling, point, 1, { nearest }, allow_duplicates);
}

template<typename Scalar, typename VectorType, typename VectorContainer>
bool check_nearest_neighbor(const VectorContainer& points, const std::vector<int>& sampling, int index, int nearest)
bool check_nearest_neighbor(const VectorContainer& points, const std::vector<int>& sampling, int index, int nearest,
bool allow_duplicates = false)
{
return check_k_nearest_neighbors<Scalar, VectorType, VectorContainer>(points, sampling, index, 1, { nearest });
return check_k_nearest_neighbors<Scalar, VectorType, VectorContainer>(points, sampling, index, 1, { nearest }, allow_duplicates);
}
1 change: 1 addition & 0 deletions tests/src/CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -23,3 +23,4 @@ add_multi_test(weight_kernel.cpp)
add_multi_test(queries_range.cpp)
add_multi_test(queries_nearest.cpp)
add_multi_test(queries_knearest.cpp)
add_multi_test(kdtree.cpp)
1 change: 1 addition & 0 deletions tests/src/algebraicsphere_primitive.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -71,6 +71,7 @@ void testFunction()
// move basis center to a portion of the radius of the sphere
f2.changeBasis(fitInitPos + Scalar(0.1)*fit.radius()*VectorType::Random());

VERIFY( fit.algebraicSphere().meanCurvature() - f2.algebraicSphere().meanCurvature() < epsilon );
VERIFY( fit.radius() - f2.radius() < epsilon );
VERIFY( (fit.center() - f2.center() ).norm() < epsilon );
VERIFY( (fit.project(candidate) - f2.project(candidate)).norm() - Scalar(1.) < epsilon );
Expand Down
2 changes: 1 addition & 1 deletion tests/src/fit_line.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -35,7 +35,7 @@ void testFunction(bool _bAddPositionNoise = false)
typedef typename DataPoint::VectorType VectorType;

//generate sampled line
int nbPoints = Eigen::internal::random<int>(100, 1000);
int nbPoints = Eigen::internal::random<int>(500, 1000);

VectorType _direction = VectorType::Random().normalized();

Expand Down
9 changes: 6 additions & 3 deletions tests/src/fit_radius_curvature_center.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -54,11 +54,14 @@ subTestSpatial<true, 3>::eval(const Fit& _fit,
Scalar kmin = _fit.kmin();
Scalar kmax = _fit.kmax();
Scalar kmean = (kmin + kmax) / Scalar (2.);
Scalar radius = Scalar(1.) / kmean;
Scalar radius1 = Scalar(1.) / kmean;
Scalar radius2 = Scalar(1.) / _fit.algebraicSphere().meanCurvature();

//Test value
VERIFY( (radius > _radiusMin - _radiusEpsilon) &&
(radius < _radiusMax + _radiusEpsilon) );
VERIFY( (radius1 > _radiusMin - _radiusEpsilon) &&
(radius1 < _radiusMax + _radiusEpsilon) );
VERIFY( (radius2 > _radiusMin - _radiusEpsilon) &&
(radius2 < _radiusMax + _radiusEpsilon) );
}


Expand Down
Loading
Loading