summaryrefslogtreecommitdiff
diff options
context:
space:
mode:
-rw-r--r--src/Spatial_searching/include/gudhi/Kd_tree_search.h98
-rw-r--r--src/Subsampling/include/gudhi/pick_n_random_points.h14
-rw-r--r--src/Subsampling/include/gudhi/sparsify_point_set.h33
-rw-r--r--src/python/gudhi/subsampling.pyx2
-rwxr-xr-xsrc/python/test/test_subsampling.py16
5 files changed, 126 insertions, 37 deletions
diff --git a/src/Spatial_searching/include/gudhi/Kd_tree_search.h b/src/Spatial_searching/include/gudhi/Kd_tree_search.h
index 87969dd9..a50a8537 100644
--- a/src/Spatial_searching/include/gudhi/Kd_tree_search.h
+++ b/src/Spatial_searching/include/gudhi/Kd_tree_search.h
@@ -12,11 +12,12 @@
#ifndef KD_TREE_SEARCH_H_
#define KD_TREE_SEARCH_H_
+#include <gudhi/Debug_utils.h>
+
#include <CGAL/Orthogonal_k_neighbor_search.h>
#include <CGAL/Orthogonal_incremental_neighbor_search.h>
#include <CGAL/Search_traits.h>
#include <CGAL/Search_traits_adapter.h>
-#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/property_map.h>
#include <CGAL/version.h> // for CGAL_VERSION_NR
@@ -40,7 +41,6 @@
namespace Gudhi {
namespace spatial_searching {
-
/**
* \class Kd_tree_search Kd_tree_search.h gudhi/Kd_tree_search.h
* \brief Spatial tree data structure to perform (approximate) nearest and furthest neighbor search.
@@ -83,7 +83,8 @@ class Kd_tree_search {
typedef CGAL::Search_traits<
FT, Point,
typename Traits::Cartesian_const_iterator_d,
- typename Traits::Construct_cartesian_const_iterator_d> Traits_base;
+ typename Traits::Construct_cartesian_const_iterator_d,
+ typename Traits::Dimension> Traits_base;
typedef CGAL::Search_traits_adapter<
std::ptrdiff_t,
@@ -110,7 +111,76 @@ class Kd_tree_search {
/// of a point P and `second` is the squared distance between P and the query point.
typedef Incremental_neighbor_search INS_range;
- typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere;
+ // Because CGAL::Fuzzy_sphere takes the radius and not its square
+ struct Sphere_for_kdtree_search
+ {
+ typedef typename Traits::Point_d Point_d;
+ typedef typename Traits::FT FT;
+ typedef typename Traits::Dimension D;
+ typedef D Dimension;
+
+ private:
+ STraits traits;
+ Point_d c;
+ FT sqradmin, sqradmax;
+ bool use_max;
+
+ public:
+ // `prefer_max` means that we prefer outputting more points at squared distance between r2min and r2max,
+ // while `!prefer_max` means we prefer fewer.
+ Sphere_for_kdtree_search(Point_d const& c_, FT const& r2min, FT const& r2max, bool prefer_max=true, STraits const& traits_ = {})
+ : traits(traits_), c(c_), sqradmin(r2min), sqradmax(r2max), use_max(prefer_max)
+ { GUDHI_CHECK(r2min >= 0 && r2max >= r2min, "0 <= r2min <= r2max"); }
+
+ bool contains(std::ptrdiff_t i) const {
+ const Point_d& p = get(traits.point_property_map(), i);
+ auto ccci = traits.construct_cartesian_const_iterator_d_object();
+ return contains_point_given_as_coordinates(ccci(p), ccci(p, 0));
+ }
+
+ template <typename Coord_iterator>
+ bool contains_point_given_as_coordinates(Coord_iterator pi, Coord_iterator CGAL_UNUSED) const {
+ FT distance = 0;
+ auto ccci = traits.construct_cartesian_const_iterator_d_object();
+ auto ci = ccci(c);
+ auto ce = ccci(c, 0);
+ FT const& limit = use_max ? sqradmax : sqradmin;
+ while (ci != ce) {
+ distance += CGAL::square(*pi++ - *ci++);
+ // I think Clément advised to check the distance at every step instead of
+ // just at the end, especially when the dimension becomes large. Distance
+ // isn't part of the concept anyway.
+ if (distance > limit) return false;
+ }
+ return true;
+ }
+
+ bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT, D> const& rect) const {
+ auto ccci = traits.construct_cartesian_const_iterator_d_object();
+ FT distance = 0;
+ auto ci = ccci(c);
+ auto ce = ccci(c, 0);
+ for (int i = 0; ci != ce; ++i, ++ci) {
+ distance += CGAL::square(CGAL::max<FT>(CGAL::max<FT>(*ci - rect.max_coord(i), rect.min_coord(i) - *ci), 0 ));
+ if (distance > sqradmin) return false;
+ }
+ return true;
+ }
+
+
+ bool outer_range_contains(CGAL::Kd_tree_rectangle<FT, D> const& rect) const {
+ auto ccci = traits.construct_cartesian_const_iterator_d_object();
+ FT distance = 0;
+ auto ci = ccci(c);
+ auto ce = ccci(c, 0);
+ for (int i = 0; ci != ce; ++i, ++ci) {
+ distance += CGAL::square(CGAL::max<FT>(*ci - rect.min_coord(i), rect.max_coord(i) - *ci));
+ if (distance > sqradmax) return false;
+ }
+ return true;
+ }
+ };
+
/// \brief Constructor
/// @param[in] points Const reference to the point range. This range
/// is not copied, so it should not be destroyed or modified afterwards.
@@ -266,10 +336,26 @@ class Kd_tree_search {
/// @param[in] eps Approximation factor.
template <typename OutputIterator>
void all_near_neighbors(Point const& p,
- FT radius,
+ FT const& radius,
OutputIterator it,
FT eps = FT(0)) const {
- m_tree.search(it, Fuzzy_sphere(p, radius, eps, m_tree.traits()));
+ all_near_neighbors2(p, CGAL::square(radius - eps), CGAL::square(radius + eps), it);
+ }
+
+ /// \brief Search for all the neighbors in a ball. This is similar to `all_near_neighbors` but takes directly
+ /// the square of the minimum distance below which points must be considered neighbors and square of the
+ /// maximum distance above which they cannot be.
+ /// @param[in] p The query point.
+ /// @param[in] sq_radius_min The square of the minimum search radius
+ /// @param[in] sq_radius_max The square of the maximum search radius
+ /// @param[out] it The points that lie inside the sphere of center `p` and squared radius `sq_radius`.
+ /// Note: `it` is used this way: `*it++ = each_point`.
+ template <typename OutputIterator>
+ void all_near_neighbors2(Point const& p,
+ FT const& sq_radius_min,
+ FT const& sq_radius_max,
+ OutputIterator it) const {
+ m_tree.search(it, Sphere_for_kdtree_search(p, sq_radius_min, sq_radius_max, true, m_tree.traits()));
}
int tree_depth() const {
diff --git a/src/Subsampling/include/gudhi/pick_n_random_points.h b/src/Subsampling/include/gudhi/pick_n_random_points.h
index a67b2b84..e4246c29 100644
--- a/src/Subsampling/include/gudhi/pick_n_random_points.h
+++ b/src/Subsampling/include/gudhi/pick_n_random_points.h
@@ -11,7 +11,9 @@
#ifndef PICK_N_RANDOM_POINTS_H_
#define PICK_N_RANDOM_POINTS_H_
-#include <gudhi/Clock.h>
+#ifdef GUDHI_SUBSAMPLING_PROFILING
+# include <gudhi/Clock.h>
+#endif
#include <boost/range/size.hpp>
@@ -44,6 +46,12 @@ void pick_n_random_points(Point_container const &points,
Gudhi::Clock t;
#endif
+ std::random_device rd;
+ std::mt19937 g(rd());
+
+#if __cplusplus >= 201703L
+ std::sample(std::begin(points), std::end(points), output_it, final_size, g);
+#else
std::size_t nbP = boost::size(points);
if (final_size > nbP)
final_size = nbP;
@@ -51,14 +59,12 @@ void pick_n_random_points(Point_container const &points,
std::vector<int> landmarks(nbP);
std::iota(landmarks.begin(), landmarks.end(), 0);
- std::random_device rd;
- std::mt19937 g(rd());
-
std::shuffle(landmarks.begin(), landmarks.end(), g);
landmarks.resize(final_size);
for (int l : landmarks)
*output_it++ = points[l];
+#endif
#ifdef GUDHI_SUBSAMPLING_PROFILING
t.end();
diff --git a/src/Subsampling/include/gudhi/sparsify_point_set.h b/src/Subsampling/include/gudhi/sparsify_point_set.h
index b30cec80..4571b8f3 100644
--- a/src/Subsampling/include/gudhi/sparsify_point_set.h
+++ b/src/Subsampling/include/gudhi/sparsify_point_set.h
@@ -11,6 +11,13 @@
#ifndef SPARSIFY_POINT_SET_H_
#define SPARSIFY_POINT_SET_H_
+#include <boost/version.hpp>
+#if BOOST_VERSION < 106600
+# include <boost/function_output_iterator.hpp>
+#else
+# include <boost/iterator/function_output_iterator.hpp>
+#endif
+
#include <gudhi/Kd_tree_search.h>
#ifdef GUDHI_SUBSAMPLING_PROFILING
#include <gudhi/Clock.h>
@@ -27,7 +34,7 @@ namespace subsampling {
* \ingroup subsampling
* \brief Outputs a subset of the input points so that the
* squared distance between any two points
- * is greater than or equal to `min_squared_dist`.
+ * is greater than `min_squared_dist`.
*
* \tparam Kernel must be a model of the <a target="_blank"
* href="http://doc.cgal.org/latest/Spatial_searching/classSearchTraits.html">SearchTraits</a>
@@ -63,29 +70,15 @@ sparsify_point_set(
// Parse the input points, and add them if they are not too close to
// the other points
std::size_t pt_idx = 0;
- for (typename Point_range::const_iterator it_pt = input_pts.begin();
- it_pt != input_pts.end();
- ++it_pt, ++pt_idx) {
- if (dropped_points[pt_idx])
+ for (auto const& pt : input_pts) {
+ if (dropped_points[pt_idx++])
continue;
- *output_it++ = *it_pt;
-
- auto ins_range = points_ds.incremental_nearest_neighbors(*it_pt);
+ *output_it++ = pt;
// If another point Q is closer that min_squared_dist, mark Q to be dropped
- for (auto const& neighbor : ins_range) {
- std::size_t neighbor_point_idx = neighbor.first;
- // If the neighbor is too close, we drop the neighbor
- if (neighbor.second < min_squared_dist) {
- // N.B.: If neighbor_point_idx < pt_idx,
- // dropped_points[neighbor_point_idx] is already true but adding a
- // test doesn't make things faster, so why bother?
- dropped_points[neighbor_point_idx] = true;
- } else {
- break;
- }
- }
+ auto drop = [&dropped_points] (std::ptrdiff_t neighbor_point_idx) { dropped_points[neighbor_point_idx] = true; };
+ points_ds.all_near_neighbors2(pt, min_squared_dist, min_squared_dist, boost::make_function_output_iterator(std::ref(drop)));
}
#ifdef GUDHI_SUBSAMPLING_PROFILING
diff --git a/src/python/gudhi/subsampling.pyx b/src/python/gudhi/subsampling.pyx
index b11d07e5..46f32335 100644
--- a/src/python/gudhi/subsampling.pyx
+++ b/src/python/gudhi/subsampling.pyx
@@ -105,7 +105,7 @@ def pick_n_random_points(points=None, off_file='', nb_points=0):
def sparsify_point_set(points=None, off_file='', min_squared_dist=0.0):
"""Outputs a subset of the input points so that the squared distance
- between any two points is greater than or equal to min_squared_dist.
+ between any two points is greater than min_squared_dist.
:param points: The input point set.
:type points: Iterable[Iterable[float]]
diff --git a/src/python/test/test_subsampling.py b/src/python/test/test_subsampling.py
index 31f64e32..4019852e 100755
--- a/src/python/test/test_subsampling.py
+++ b/src/python/test/test_subsampling.py
@@ -141,12 +141,16 @@ def test_simple_sparsify_points():
# assert gudhi.sparsify_point_set(points = [], min_squared_dist = 0.0) == []
# assert gudhi.sparsify_point_set(points = [], min_squared_dist = 10.0) == []
assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=0.0) == point_set
- assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.0) == point_set
- assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.0) == [
+ assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=0.999) == point_set
+ assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.001) == [
[0, 1],
[1, 0],
]
- assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.01) == [[0, 1]]
+ assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.999) == [
+ [0, 1],
+ [1, 0],
+ ]
+ assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.001) == [[0, 1]]
assert (
len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=0.0))
@@ -157,11 +161,11 @@ def test_simple_sparsify_points():
== 5
)
assert (
- len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.0))
+ len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.1))
== 4
)
assert (
- len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=90.0))
+ len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=89.9))
== 3
)
assert (
@@ -169,7 +173,7 @@ def test_simple_sparsify_points():
== 2
)
assert (
- len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=325.0))
+ len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=324.9))
== 2
)
assert (