songsenand
/
cgal
mirror of https://github.com/CGAL/cgal
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

Merge pull request #2548 from cjamin/Faster_dd_spatial_searching-cjamin 

Spatial_searching: faster Kd_tree
This commit is contained in:
Sebastien Loriot 2019-12-20 09:27:09 +01:00 committed by GitHub
parent 58f8a3eec2 c328c66299
commit a4551ab619
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
30 changed files with 1347 additions and 422 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Classification/include/CGAL/Classification/Point_set_neighborhood.h
View File

@ -81,7 +81,7 @@ class Point_set_neighborhood
typedef Search_traits_adapter <boost::uint32_t, My_point_property_map, SearchTraits_3> Search_traits;
typedef Sliding_midpoint<Search_traits> Splitter;
typedef Distance_adapter<boost::uint32_t, My_point_property_map, Euclidean_distance<SearchTraits_3> > Distance;
typedef Kd_tree<Search_traits, Splitter, Tag_true> Tree;
typedef Kd_tree<Search_traits, Splitter, Tag_true, Tag_true> Tree;
typedef Fuzzy_sphere<Search_traits> Sphere;
typedef Orthogonal_k_neighbor_search<Search_traits, Distance, Splitter, Tree> Knn;

23
Installation/CHANGES.md
View File

@ -27,6 +27,29 @@ Release date: June 2020
does not allow any intersection, except for the configuration of two constraints having a single
common endpoints, for convience.
### dD Spatial Searching
- Improved the performance of the kd-tree in some cases:
- Not storing the points coordinates inside the tree usually
generates a lot of cache misses, leading to non-optimal
performance. This is the case for example
when indices are stored inside the tree, or to a lesser extent when the points
coordinates are stored in a dynamically allocated array (e.g., `Epick_d`
with dynamic dimension) &mdash; we says "to a lesser extent" because the points
are re-created by the kd-tree in a cache-friendly order after its construction,
so the coordinates are more likely to be stored in a near-optimal order
on the heap.
In these cases, the new `EnablePointsCache` template parameter of the
`CGAL::Kd_tree` class can be set to `CGAL::Tag_true`. The points coordinates
will then be cached in an optimal way. This will increase memory
consumption but provides better search performance. See the updated
`GeneralDistance` and `FuzzyQueryItem`
concepts for additional requirements when using such a cache.
- In most cases (e.g., Euclidean distance), the distance computation
algorithm knows before its end that the distance will be greater
than or equal to some given value. This is used in the (orthogonal)
k-NN search to interrupt some distance computations before its end,
saving precious milliseconds, in particular in medium-to-high dimension.
[Release 5.0](https://github.com/CGAL/cgal/releases/tag/releases%2FCGAL-5.0)
-----------

2
Point_set_processing_3/include/CGAL/Point_set_processing_3/internal/Search_traits_vertex_handle_3.h
View File

@ -117,8 +117,10 @@ struct Construct_cartesian_const_iterator_vertex_handle_3
template <class Vertex_handle>
struct Euclidean_distance_vertex_handle_3
{
typedef double FT;
typedef CGAL::internal::Point_vertex_handle_3<Vertex_handle> Point_vertex_handle_3;
typedef Point_vertex_handle_3 Query_item;
typedef Point_vertex_handle_3 Point_d;
double transformed_distance(const Point_vertex_handle_3& p1, const Point_vertex_handle_3& p2) const {
double distx = p1.x()-p2.x();

19
STL_Extension/include/CGAL/Compact_container.h
View File

@ -30,6 +30,7 @@
#include <CGAL/iterator.h>
#include <CGAL/CC_safe_handle.h>
#include <CGAL/Time_stamper.h>
#include <CGAL/Has_member.h>
#include <boost/mpl/if.hpp>
@ -82,24 +83,6 @@
namespace CGAL {
#define CGAL_GENERATE_MEMBER_DETECTOR(X) \
template<typename T> class has_##X { \
struct Fallback { int X; }; \
struct Derived : T, Fallback { }; \
\
template<typename U, U> struct Check; \
\
typedef char ArrayOfOne[1]; \
typedef char ArrayOfTwo[2]; \
\
template<typename U> static ArrayOfOne & func( \
Check<int Fallback::*, &U::X> *); \
template<typename U> static ArrayOfTwo & func(...); \
public: \
typedef has_##X type; \
enum { value = sizeof(func<Derived>(0)) == 2 }; \
} // semicolon is after the macro call
#define CGAL_INIT_COMPACT_CONTAINER_BLOCK_SIZE 14
#define CGAL_INCREMENT_COMPACT_CONTAINER_BLOCK_SIZE 16

2
Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/K_neighbor_query.h
View File

@ -95,7 +95,7 @@ namespace Point_set {
CGAL::Sliding_midpoint<Search_traits>;
using Search_tree =
CGAL::Kd_tree<Search_traits, Splitter, CGAL::Tag_true>;
CGAL::Kd_tree<Search_traits, Splitter, CGAL::Tag_true, CGAL::Tag_true>;
using Neighbor_search =
CGAL::Orthogonal_k_neighbor_search<

2
Shape_detection/include/CGAL/Shape_detection/Region_growing/Region_growing_on_point_set/Sphere_neighbor_query.h
View File

@ -101,7 +101,7 @@ namespace Point_set {
= CGAL::Fuzzy_sphere<Search_traits>;
using Tree
= CGAL::Kd_tree<Search_traits, Splitter, CGAL::Tag_true>;
= CGAL::Kd_tree<Search_traits, Splitter, CGAL::Tag_true, CGAL::Tag_true>;
/// \endcond
/// @}

18
Spatial_searching/doc/Spatial_searching/CGAL/Euclidean_distance.h
View File

@ -64,6 +64,24 @@ Returns the squared Euclidean distance between `q` and `p`.
*/
FT transformed_distance(Query_item q, Point_d p) const;
/*!
Returns the transformed distance between `q` and the point whose Cartesian
coordinates are contained in the range [`begin`, `end`).
*/
template <typename Coord_iterator>
FT transformed_distance_from_coordinates(
Query_item q, Coord_iterator begin, Coord_iterator end) const;
/*!
Returns the transformed distance between `q` and the point whose Cartesian
coordinates are contained in the range [`begin`, `end`), or any value
\f$ \geq \f$ `stop_if_geq_to_this` if the transformed distance is
\f$ \geq \f$ `stop_if_geq_to_this`.
*/
template <typename Coord_iterator>
FT interruptible_transformed_distance(
Query_item q, Coord_iterator begin, Coord_iterator end, FT stop_if_geq_to_this) const;
/*!
Returns the squared Euclidean distance between `q` and
the point on the boundary of `r` closest to `q`.

7
Spatial_searching/doc/Spatial_searching/CGAL/Fuzzy_iso_box.h
View File

@ -73,6 +73,13 @@ Test whether the fuzzy iso box contains `p`.
*/
bool contains(Point_d p) const;
/*!
Test whether the fuzzy iso box contains the point whose Cartesian
coordinates are contained in the range [`begin`, `end`).
*/
template <typename Coord_iterator>
bool contains_point_given_as_coordinates(Coord_iterator begin, Coord_iterator end) const;
/*!
Test whether the inner box intersects the rectangle
associated with a node of a tree.

7
Spatial_searching/doc/Spatial_searching/CGAL/Fuzzy_sphere.h
View File

@ -75,6 +75,13 @@ less than \f$ r\f$.
*/
bool contains(const Point_d& p) const;
/*!
Test whether the fuzzy sphere contains the point whose Cartesian coordinates
are contained in the range [`begin`, `end`).
*/
template <typename Coord_iterator>
bool contains_point_given_as_coordinates(Coord_iterator begin, Coord_iterator end) const;
/*!
Test whether the inner sphere intersects the rectangle
associated with a node of a tree.

4
Spatial_searching/doc/Spatial_searching/CGAL/Incremental_neighbor_search.h
View File

@ -20,8 +20,8 @@ and `Euclidean_distance<Traits>` otherwise.
The default type is `Sliding_midpoint<Traits>`.
\tparam SpatialTree must be a model of the concept `SpatialTree`.
The default type is `Kd_tree<Traits,Splitter,Tag_false>`. The
template argument `Tag_false` makes that the tree is built with unextended nodes.
The default type is `Kd_tree<Traits, Splitter, Tag_false, Tag_false>`. The
third template argument `Tag_false` makes that the tree is built with unextended nodes.
\sa `CGAL::Orthogonal_incremental_neighbor_search<Traits, OrthogonalDistance, Splitter, SpatialTree>`

7
Spatial_searching/doc/Spatial_searching/CGAL/K_neighbor_search.h
View File

@ -11,14 +11,17 @@ extended or unextended nodes.
\tparam Traits must be an implementation of the concept `SearchTraits`,
for example `Simple_cartesian<double>`.
\tparam Splitter must be a model of the
\tparam GeneralDistance must be a model of the
concept `GeneralDistance`. If `Traits` is
`Search_traits_adapter<Key,PointPropertyMap,BaseTraits>`
the default type is `Distance_adapter<Key,PointPropertyMap,Euclidean_distance<BaseTraits> >`,
and `Euclidean_distance<Traits>` otherwise.
\tparam Splitter must be a model of the concept `Splitter`.
The default type is `Sliding_midpoint<Traits>`.
\tparam SpatialTree must be an implementation of the concept `SpatialTree`.
The default type is `Kd_tree<Traits, Splitter, Tag_false>`. The
The default type is `Kd_tree<Traits, Splitter, Tag_false, Tag_false>`. The third
template argument `Tag_false` makes that the tree is built with unextended nodes.
\sa `CGAL::Orthogonal_k_neighbor_search<Traits, OrthogonalDistance, Splitter, SpatialTree>`

17
Spatial_searching/doc/Spatial_searching/CGAL/Kd_tree.h
View File

@ -14,13 +14,28 @@ It defaults to `Sliding_midpoint<Traits>`.
\tparam UseExtendedNode must be `Tag_true`, if the
tree shall be built with extended nodes, and `Tag_false` otherwise.
\tparam EnablePointsCache can be `Tag_true` or `Tag_false`.
Not storing the points coordinates inside the tree usually generates a
lot of cache misses, leading to non-optimal performance. This is the case
for example when indices are stored inside the tree,
or to a lesser extent when the points coordinates are stored
in a dynamically allocated array (e.g., `Epick_d` with dynamic
dimension) &mdash; we says "to a lesser extent" because the points
are re-created by the kd-tree in a cache-friendly order after its construction,
so the coordinates are more likely to be stored in a near-optimal order on the
heap. When EnablePointsCache` is set to `Tag_true`, the points
coordinates will be cached in an optimal way. This will
increase memory consumption but provide better search performance.
See also the `GeneralDistance` and `FuzzyQueryItem` concepts for
additional requirements when using such a cache.
\sa `CGAL::Kd_tree_node<Traits>`
\sa `CGAL::Search_traits_2<Kernel>`
\sa `CGAL::Search_traits_3<Kernel>`
\sa `CGAL::Search_traits<FT_,Point,CartesianIterator,ConstructCartesianIterator>`
*/
template< typename Traits, typename Splitter, typename UseExtendedNode >
template< typename Traits, typename Splitter, typename UseExtendedNode, typename EnablePointsCache >
class Kd_tree {
public:

2
Spatial_searching/doc/Spatial_searching/CGAL/Orthogonal_incremental_neighbor_search.h
View File

@ -18,7 +18,7 @@ and `Euclidean_distance<Traits>` otherwise.
The default type is `Sliding_midpoint<Traits>`.
\tparam SpatialTree must be a model of the concept `SpatialTree`.
The default type is `Kd_tree<Traits, Splitter, Tag_true>`. The
The default type is `Kd_tree<Traits, Splitter, Tag_true, Tag_false>`. The third
template argument must be `Tag_true` because orthogonal search needs extended
kd tree nodes.

4
Spatial_searching/doc/Spatial_searching/CGAL/Orthogonal_k_neighbor_search.h
View File

@ -16,13 +16,11 @@ concept `OrthogonalDistance`. If `Traits` is
the default type is `Distance_adapter<Key,PointPropertyMap,Euclidean_distance<BaseTraits> >`,
and `Euclidean_distance<Traits>` otherwise.
The default type is `Euclidean_distance<Traits>`.
\tparam Splitter must be a model of the concept `Splitter`.
The default type is `Sliding_midpoint<Traits>`.
\tparam SpatialTree must be a model of the concept `SpatialTree`.
The default type is `Kd_tree<Traits, Splitter, Tag_true>`. The
The default type is `Kd_tree<Traits, Splitter, Tag_true, Tag_false>`. The third
template argument must be `Tag_true` because orthogonal search needs extended
kd tree nodes.

9
Spatial_searching/doc/Spatial_searching/Concepts/FuzzyQueryItem.h
View File

@ -40,6 +40,15 @@ Test whether the query item contains `p`.
*/
bool contains(Point_d p) const;
/*!
Optional: must be defined when used with a `Kd_tree` where `EnablePointsCache`
is set to `Tag_true`.
Test whether the query item contains the point whose Cartesian coordinates
are contained in the range [`begin`, `end`).
*/
template <typename Coord_iterator>
bool contains_point_given_as_coordinates(Coord_iterator begin, Coord_iterator end) const;
/*!
Test whether the inner approximation of the spatial object intersects a rectangle
associated with a node of a tree.

28
Spatial_searching/doc/Spatial_searching/Concepts/GeneralDistance.h
View File

@ -49,6 +49,34 @@ Returns the transformed distance between `q` and `r`.
*/
FT transformed_distance(Query_item q, Point_d r);
/*!
Optional: must be defined when used with a `Kd_tree` where `EnablePointsCache`
is set to `Tag_true`.
Returns the transformed distance between `q` and the point whose Cartesian
coordinates are contained in the range [`begin`, `end`).
*/
template <typename Coord_iterator>
FT transformed_distance_from_coordinates(
Query_item q, Coord_iterator begin, Coord_iterator end) const;
/*!
Optional: in most cases (e.g., Euclidean distance), the distance computation
algorithm knows before its end that the distance will be greater than or equal
to some given value. In this function, the computation can be stopped when the
distance is going to be greater than or equal to `stop_if_geq_to_this`. In this case,
the only requirement of the return value it to be \f$ \geq \f$ `stop_if_geq_to_this`.
Note that points cache does not have to be activated to enable this optimization.
Returns the transformed distance between `q` and the point whose Cartesian
coordinates are contained in the range [`begin`, `end`), or any value
\f$ \geq \f$ `stop_if_geq_to_this` if the transformed distance is
\f$ \geq \f$ `stop_if_geq_to_this`.
*/
template <typename Coord_iterator>
FT interruptible_transformed_distance(
Query_item q, Coord_iterator begin, Coord_iterator end, FT stop_if_geq_to_this) const;
/*!
Returns the transformed distance between `q` and
the point on the boundary of `r` closest to `q`.

15
Spatial_searching/doc/Spatial_searching/Spatial_searching.txt
View File

@ -503,6 +503,21 @@ instead of the distance itself. For instance for the Euclidean
distance, to avoid the expensive computation of square roots, squared
distances are used instead of the Euclidean distance itself.
Not storing the points coordinates inside the tree usually generates a
lot of cache misses, leading to non-optimal performance. This is the case
for example when indices are stored inside the tree,
or to a lesser extent when the points coordinates are stored
in a dynamically allocated array (e.g., `Epick_d` with dynamic
dimension) &mdash; we says "to a lesser extent" because the points
are re-created by the kd-tree in a cache-friendly order after its construction,
so the coordinates are more likely to be stored in a near-optimal order on the
heap.
In this case, the `EnablePointsCache` template parameter of the `Kd_tree` class can be
set to `Tag_true`. The points coordinates will then be cached in an optimal
way. This will increase memory consumption but provide better search
performance. See also the `GeneralDistance` and `FuzzyQueryItem` concepts for
additional requirements when using such a cache.
\section Spatial_searchingImplementationHistory Implementation History
The initial implementation of this package was done by Hans Tangelder and

1
Spatial_searching/examples/Spatial_searching/Distance.h
View File

@ -1,5 +1,6 @@
struct Distance {
typedef Point Query_item;
typedef Point Point_d;
typedef double FT;
typedef CGAL::Dimension_tag<3> D;

206
Spatial_searching/include/CGAL/Euclidean_distance.h
View File

@ -9,6 +9,7 @@
//
//
// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>)
// Clement Jamin (clement.jamin.pro@gmail.com)
#ifndef CGAL_EUCLIDEAN_DISTANCE_H
@ -21,6 +22,7 @@
#include <CGAL/number_utils.h>
#include <CGAL/internal/Get_dimension_tag.h>
#include <vector>
#include <iterator>
namespace CGAL {
@ -52,92 +54,146 @@ namespace CGAL {
// default constructor
Euclidean_distance(const SearchTraits& traits_=SearchTraits()):traits(traits_) {}
inline FT transformed_distance(const Query_item& q, const Point_d& p) const {
return transformed_distance(q,p, D());
inline FT transformed_distance(const Query_item& q, const Point_d& p) const
{
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d p_begin = construct_it(p), p_end = construct_it(p, 0);
return transformed_distance_from_coordinates(q, p_begin, p_end);
}
//Dynamic version for runtime dimension
inline FT transformed_distance(const Query_item& q, const Point_d& p, Dynamic_dimension_tag) const {
template <typename Coord_iterator>
inline FT transformed_distance_from_coordinates(const Query_item& q,
Coord_iterator it_coord_begin, Coord_iterator it_coord_end) const
{
return transformed_distance_from_coordinates(q, it_coord_begin, it_coord_end, D());
}
// Static dim = 2 loop unrolled
template <typename Coord_iterator>
inline FT transformed_distance_from_coordinates(const Query_item& q,
Coord_iterator it_coord_begin, Coord_iterator /*unused*/,
Dimension_tag<2>) const
{
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q);
FT distance = square(*qit - *it_coord_begin);
qit++; it_coord_begin++;
distance += square(*qit - *it_coord_begin);
return distance;
}
// Static dim = 3 loop unrolled
template <typename Coord_iterator>
inline FT transformed_distance_from_coordinates(const Query_item& q,
Coord_iterator it_coord_begin, Coord_iterator /*unused*/,
Dimension_tag<3>) const
{
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q);
FT distance = square(*qit - *it_coord_begin);
qit++; it_coord_begin++;
distance += square(*qit - *it_coord_begin);
qit++; it_coord_begin++;
distance += square(*qit - *it_coord_begin);
return distance;
}
// Other cases: static dim > 3 or dynamic dim
template <typename Coord_iterator, typename Dim>
inline FT transformed_distance_from_coordinates(const Query_item& q,
Coord_iterator it_coord_begin, Coord_iterator /*unused*/,
Dim) const
{
FT distance = FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
qe = construct_it(q,1), pit = construct_it(p);
for(; qit != qe; qit++, pit++){
distance += ((*qit)-(*pit))*((*qit)-(*pit));
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), qe = construct_it(q, 1);
for (; qit != qe; ++qit, ++it_coord_begin)
{
FT diff = (*qit) - (*it_coord_begin);
distance += diff*diff;
}
return distance;
}
//Generic version for DIM > 3
template < int DIM >
inline FT transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<DIM>) const {
// During the computation, if the partially-computed distance `pcd` gets greater or equal
// to `stop_if_geq_to_this`, the computation is stopped and `pcd` is returned
template <typename Coord_iterator>
inline FT interruptible_transformed_distance(const Query_item& q,
Coord_iterator it_coord_begin, Coord_iterator /*unused*/,
FT stop_if_geq_to_this) const
{
FT distance = FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
qe = construct_it(q,1), pit = construct_it(p);
for(; qit != qe; qit++, pit++){
distance += ((*qit)-(*pit))*((*qit)-(*pit));
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q), qe = construct_it(q, 1);
if (std::distance(qit, qe) >= 6)
{
// Every 4 coordinates, the current partially-computed distance
// is compared to stop_if_geq_to_this
// Note: the concept SearchTraits specifies that Cartesian_const_iterator_d
// must be a random-access iterator
typename SearchTraits::Cartesian_const_iterator_d qe_minus_5 = qe;
std::advance(qe, -5);
for (;;)
{
FT diff = (*qit) - (*it_coord_begin);
distance += diff*diff;
++qit; ++it_coord_begin;
diff = (*qit) - (*it_coord_begin);
distance += diff*diff;
++qit; ++it_coord_begin;
diff = (*qit) - (*it_coord_begin);
distance += diff*diff;
++qit; ++it_coord_begin;
diff = (*qit) - (*it_coord_begin);
distance += diff*diff;
++qit, ++it_coord_begin;
if (distance >= stop_if_geq_to_this)
return distance;
if (std::distance(qe_minus_5, qit) >= 0)
break;
}
}
for (; qit != qe; ++qit, ++it_coord_begin)
{
FT diff = (*qit) - (*it_coord_begin);
distance += diff*diff;
}
return distance;
}
//DIM = 2 loop unrolled
inline FT transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<2> ) const {
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),pit = construct_it(p);
FT distance = square(*qit - *pit);
qit++;pit++;
distance += square(*qit - *pit);
return distance;
}
//DIM = 3 loop unrolled
inline FT transformed_distance(const Query_item& q, const Point_d& p, Dimension_tag<3> ) const {
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),pit = construct_it(p);
FT distance = square(*qit - *pit);
qit++;pit++;
distance += square(*qit - *pit);
qit++;pit++;
distance += square(*qit - *pit);
return distance;
}
inline FT min_distance_to_rectangle(const Query_item& q,
const Kd_tree_rectangle<FT,D>& r) const {
const Kd_tree_rectangle<FT, D>& r) const {
FT distance = FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
qe = construct_it(q,1);
for(unsigned int i = 0;qit != qe; i++, qit++){
if((*qit) < r.min_coord(i))
qe = construct_it(q, 1);
for (unsigned int i = 0; qit != qe; i++, qit++) {
if ((*qit) < r.min_coord(i))
distance +=
(r.min_coord(i)-(*qit))*(r.min_coord(i)-(*qit));
(r.min_coord(i) - (*qit))*(r.min_coord(i) - (*qit));
else if ((*qit) > r.max_coord(i))
distance +=
((*qit)-r.max_coord(i))*((*qit)-r.max_coord(i));
((*qit) - r.max_coord(i))*((*qit) - r.max_coord(i));
}
return distance;
}
inline FT min_distance_to_rectangle(const Query_item& q,
const Kd_tree_rectangle<FT,D>& r,std::vector<FT>& dists) const {
const Kd_tree_rectangle<FT, D>& r, std::vector<FT>& dists) const {
FT distance = FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
qe = construct_it(q,1);
for(unsigned int i = 0;qit != qe; i++, qit++){
if((*qit) < r.min_coord(i)){
dists[i] = (r.min_coord(i)-(*qit));
qe = construct_it(q, 1);
for (unsigned int i = 0; qit != qe; i++, qit++) {
if ((*qit) < r.min_coord(i)) {
dists[i] = (r.min_coord(i) - (*qit));
distance += dists[i] * dists[i];
}
else if ((*qit) > r.max_coord(i)){
dists[i] = ((*qit)-r.max_coord(i));
else if ((*qit) > r.max_coord(i)) {
dists[i] = ((*qit) - r.max_coord(i));
distance += dists[i] * dists[i];
}
@ -146,33 +202,33 @@ namespace CGAL {
}
inline FT max_distance_to_rectangle(const Query_item& q,
const Kd_tree_rectangle<FT,D>& r) const {
FT distance=FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
const Kd_tree_rectangle<FT, D>& r) const {
FT distance = FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
qe = construct_it(q,1);
for(unsigned int i = 0;qit != qe; i++, qit++){
if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0))
distance += (r.max_coord(i)-(*qit))*(r.max_coord(i)-(*qit));
qe = construct_it(q, 1);
for (unsigned int i = 0; qit != qe; i++, qit++) {
if ((*qit) <= (r.min_coord(i) + r.max_coord(i)) / FT(2.0))
distance += (r.max_coord(i) - (*qit))*(r.max_coord(i) - (*qit));
else
distance += ((*qit)-r.min_coord(i))*((*qit)-r.min_coord(i));
distance += ((*qit) - r.min_coord(i))*((*qit) - r.min_coord(i));
};
return distance;
}
inline FT max_distance_to_rectangle(const Query_item& q,
const Kd_tree_rectangle<FT,D>& r,std::vector<FT>& dists ) const {
FT distance=FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
const Kd_tree_rectangle<FT, D>& r, std::vector<FT>& dists) const {
FT distance = FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d qit = construct_it(q),
qe = construct_it(q,1);
for(unsigned int i = 0;qit != qe; i++, qit++){
if ((*qit) <= (r.min_coord(i)+r.max_coord(i))/FT(2.0)){
dists[i] = (r.max_coord(i)-(*qit));
qe = construct_it(q, 1);
for (unsigned int i = 0; qit != qe; i++, qit++) {
if ((*qit) <= (r.min_coord(i) + r.max_coord(i)) / FT(2.0)) {
dists[i] = (r.max_coord(i) - (*qit));
distance += dists[i] * dists[i];
}
else{
dists[i] = ((*qit)-r.min_coord(i));
else {
dists[i] = ((*qit) - r.min_coord(i));
distance += dists[i] * dists[i];
}
};

12
Spatial_searching/include/CGAL/Fuzzy_iso_box.h
View File

@ -120,6 +120,17 @@ namespace CGAL {
return true;
}
template <typename Coord_iterator>
bool contains_point_given_as_coordinates(Coord_iterator it_coord_begin, Coord_iterator /*unused*/) const {
Construct_cartesian_const_iterator_d construct_it=traits.construct_cartesian_const_iterator_d_object();
Cartesian_const_iterator_d minit= min_begin, maxit = max_begin;
for (unsigned int i = 0; i < dim; ++i, ++it_coord_begin, ++minit, ++maxit) {
if ( ((*it_coord_begin) < (*minit)) || ((*it_coord_begin) > (*maxit)) )
return false;
}
return true;
}
bool inner_range_intersects(const Kd_tree_rectangle<FT,Dimension>& rectangle) const {
// test whether the box eroded by 'eps' intersects 'rectangle'
Cartesian_const_iterator_d minit= min_begin, maxit = max_begin;
@ -131,7 +142,6 @@ namespace CGAL {
return true;
}
bool outer_range_contains(const Kd_tree_rectangle<FT,Dimension>& rectangle) const {
// test whether the box dilated by 'eps' contains 'rectangle'
Cartesian_const_iterator_d minit= min_begin, maxit = max_begin;

15
Spatial_searching/include/CGAL/Fuzzy_sphere.h
View File

@ -72,6 +72,21 @@ public:
return (distance <= sq_radius);
}
template <typename Coord_iterator>
bool contains_point_given_as_coordinates(Coord_iterator it_coord_begin, Coord_iterator /*unused*/) const {
// test whether the distance between c and p is less than the radius
FT distance = FT(0);
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it =
traits.construct_cartesian_const_iterator_d_object();
typename SearchTraits::Cartesian_const_iterator_d cit = construct_it(c),
end = construct_it(c, 0);
for (; cit != end && (distance <= sq_radius); ++cit, ++it_coord_begin)
distance += ((*cit) - (*it_coord_begin))*((*cit) - (*it_coord_begin));
return (distance < sq_radius);
}
bool inner_range_intersects(const Kd_tree_rectangle<FT,Dimension>& rectangle) const {
// test whether the sphere with radius (r-eps) intersects 'rectangle', i.e.
// if the minimal distance of c to 'rectangle' is less than r-eps

125
Spatial_searching/include/CGAL/Incremental_neighbor_search.h
View File

@ -9,6 +9,7 @@
//
//
// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>)
// Clement Jamin (clement.jamin.pro@gmail.com)
#ifndef CGAL_INCREMENTAL_NEIGHBOR_SEARCH_H
#define CGAL_INCREMENTAL_NEIGHBOR_SEARCH_H
@ -16,21 +17,23 @@
#include <CGAL/license/Spatial_searching.h>
#include <CGAL/disable_warnings.h>
#include <CGAL/Kd_tree_node.h>
#include <CGAL/Kd_tree_rectangle.h>
#include <CGAL/Euclidean_distance.h>
#include <CGAL/internal/Search_helpers.h>
#include <cstring>
#include <list>
#include <queue>
#include <memory>
#include <CGAL/Kd_tree_node.h>
#include <CGAL/Kd_tree_rectangle.h>
#include <CGAL/Euclidean_distance.h>
#include <iterator> // for std::distance
namespace CGAL {
template <class SearchTraits,
class Distance_=typename internal::Spatial_searching_default_distance<SearchTraits>::type,
class Splitter_ = Sliding_midpoint<SearchTraits>,
class Tree_=Kd_tree<SearchTraits, Splitter_, Tag_false> >
class Tree_=Kd_tree<SearchTraits, Splitter_, Tag_false, Tag_false> >
class Incremental_neighbor_search {
public:
@ -264,6 +267,9 @@ namespace CGAL {
FT rd;
internal::Distance_helper<Distance, SearchTraits> m_distance_helper;
int m_dim;
Tree const& m_tree;
class Priority_higher {
@ -324,6 +330,8 @@ namespace CGAL {
Iterator_implementation(const Tree& tree, const Query_item& q,const Distance& tr,
FT Eps, bool search_nearest)
: query_point(q), search_nearest_neighbour(search_nearest),
m_distance_helper(tr, tree.traits()),
m_tree(tree),
PriorityQueue(Priority_higher(search_nearest)),
Item_PriorityQueue(Distance_smaller(search_nearest)),
distance(tr), reference_count(1), number_of_internal_nodes_visited(0),
@ -332,6 +340,10 @@ namespace CGAL {
{
if (tree.empty()) return;
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it =
m_tree.traits().construct_cartesian_const_iterator_d_object();
m_dim = static_cast<int>(std::distance(construct_it(q), construct_it(q, 0)));
multiplication_factor= distance.transformed_distance(FT(1)+Eps);
Node_box *bounding_box = new Node_box((tree.bounding_box()));
@ -422,6 +434,81 @@ namespace CGAL {
delete The_item_top;
}
// With cache
bool search_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_true, bool search_furthest)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
typename std::vector<FT>::const_iterator cache_point_begin = m_tree.cache_begin() + m_dim*(it_node_point - m_tree.begin());
for (; it_node_point != node->end(); ++it_node_point)
{
number_of_items_visited++;
FT distance_to_query_point =
m_distance_helper.transformed_distance_from_coordinates(
query_point, *it_node_point, cache_point_begin, cache_point_begin + m_dim);
Point_with_transformed_distance *NN_Candidate =
new Point_with_transformed_distance(*it_node_point, distance_to_query_point);
Item_PriorityQueue.push(NN_Candidate);
cache_point_begin += m_dim;
}
// old top of PriorityQueue has been processed,
// hence update rd
bool next_neighbour_found;
if (!(PriorityQueue.empty()))
{
rd = PriorityQueue.top()->second;
next_neighbour_found = (search_furthest ?
multiplication_factor*rd < Item_PriorityQueue.top()->second
: multiplication_factor*rd > Item_PriorityQueue.top()->second);
}
else // priority queue empty => last neighbour found
{
next_neighbour_found = true;
}
number_of_neighbours_computed++;
return next_neighbour_found;
}
// Without cache
bool search_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_false, bool search_furthest)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
for (; it_node_point != it_node_point_end; ++it_node_point)
{
number_of_items_visited++;
FT distance_to_query_point =
distance.transformed_distance(query_point, *it_node_point);
Point_with_transformed_distance *NN_Candidate =
new Point_with_transformed_distance(*it_node_point, distance_to_query_point);
Item_PriorityQueue.push(NN_Candidate);
}
// old top of PriorityQueue has been processed,
// hence update rd
bool next_neighbour_found;
if (!(PriorityQueue.empty()))
{
rd = PriorityQueue.top()->second;
next_neighbour_found = (search_furthest ?
multiplication_factor*rd < Item_PriorityQueue.top()->second
: multiplication_factor*rd > Item_PriorityQueue.top()->second);
}
else // priority queue empty => last neighbour found
{
next_neighbour_found = true;
}
number_of_neighbours_computed++;
return next_neighbour_found;
}
void
Compute_the_next_nearest_neighbour()
{
@ -501,36 +588,14 @@ namespace CGAL {
}
}
delete B;
// N is a leaf
typename Tree::Leaf_node_const_handle node =
static_cast<typename Tree::Leaf_node_const_handle>(N);
number_of_leaf_nodes_visited++;
if (node->size() > 0) {
for (typename Tree::iterator it = node->begin(); it != node->end(); it++) {
number_of_items_visited++;
FT distance_to_query_point=
distance.transformed_distance(query_point,*it);
Point_with_transformed_distance *NN_Candidate=
new Point_with_transformed_distance(*it,distance_to_query_point);
Item_PriorityQueue.push(NN_Candidate);
}
// old top of PriorityQueue has been processed,
// hence update rd
if (!PriorityQueue.empty()) {
rd = PriorityQueue.top()->second;
if (search_nearest_neighbour)
next_neighbour_found =
(multiplication_factor*rd >
Item_PriorityQueue.top()->second);
else
next_neighbour_found =
(multiplication_factor*rd <
Item_PriorityQueue.top()->second);
}
else // priority queue empty => last neighbour found
{
next_neighbour_found=true;
};
number_of_neighbours_computed++;
typename internal::Has_points_cache<Tree, internal::has_Enable_points_cache<Tree>::type::value>::type dummy;
next_neighbour_found = search_in_leaf(node, dummy, !search_nearest_neighbour);
}
} // next_neighbour_found or priority queue is empty
// in the latter case also the item priority queue is empty

138
Spatial_searching/include/CGAL/K_neighbor_search.h
View File

@ -8,7 +8,8 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>)
// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>),
// Clement Jamin (clement.jamin.pro@gmail.com)
#ifndef CGAL_K_NEIGHBOR_SEARCH_H
#define CGAL_K_NEIGHBOR_SEARCH_H
@ -19,8 +20,9 @@
#include <CGAL/internal/K_neighbor_search.h>
#include <CGAL/internal/Get_dimension_tag.h>
#include <CGAL/internal/Search_helpers.h>
#include <iterator> // for std::distance
namespace CGAL {
@ -30,9 +32,11 @@ class K_neighbor_search;
template <class SearchTraits,
class Distance= typename internal::Spatial_searching_default_distance<SearchTraits>::type,
class Splitter= Sliding_midpoint<SearchTraits> ,
class Tree= Kd_tree<SearchTraits, Splitter, Tag_true> >
class K_neighbor_search: public internal::K_neighbor_search<SearchTraits,Distance,Splitter,Tree> {
class Tree= Kd_tree<SearchTraits, Splitter, Tag_true, Tag_false> >
class K_neighbor_search: public internal::K_neighbor_search<SearchTraits,Distance,Splitter,Tree>
{
typedef internal::K_neighbor_search<SearchTraits,Distance,Splitter,Tree> Base;
typedef typename Tree::Point_d Point;
public:
typedef typename Base::FT FT;
@ -40,9 +44,12 @@ public:
K_neighbor_search(const Tree& tree, const typename Base::Query_item& q,
unsigned int k=1, FT Eps=FT(0.0), bool Search_nearest=true, const Distance& d=Distance(),bool sorted=true)
: Base(q,k,Eps,Search_nearest,d)
: Base(q,k,Eps,Search_nearest,d),
m_distance_helper(this->distance_instance, tree.traits()),
m_tree(tree)
{
if (tree.empty()) return;
compute_neighbors_general(tree.root(),tree.bounding_box());
if (sorted) this->queue.sort();
};
@ -51,6 +58,9 @@ private:
typedef typename Base::Node_const_handle Node_const_handle;
using Base::branch;
internal::Distance_helper<Distance, SearchTraits> m_distance_helper;
Tree const& m_tree;
void
compute_neighbors_general(typename Base::Node_const_handle N, const Kd_tree_rectangle<FT,D>& r)
{
@ -121,15 +131,123 @@ private:
static_cast<typename Tree::Leaf_node_const_handle>(N);
this->number_of_leaf_nodes_visited++;
if (node->size() > 0)
for (typename Tree::iterator it = node->begin(); it != node->end(); it++) {
this->number_of_items_visited++;
FT distance_to_query_object =
this->distance_instance.transformed_distance(this->query_object,*it);
this->queue.insert(std::make_pair(&(*it),distance_to_query_object));
{
typename internal::Has_points_cache<Tree, internal::has_Enable_points_cache<Tree>::type::value>::type dummy;
if (this->search_nearest)
search_nearest_in_leaf(node, dummy);
else
search_furthest_in_leaf(node, dummy);
}
}
}
// With cache
void search_nearest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_true)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
int dim = m_tree.dim();
typename std::vector<FT>::const_iterator cache_point_begin = m_tree.cache_begin() + dim*(it_node_point - m_tree.begin());
// As long as the queue is not full, the node is just inserted
for (; !this->queue.full() && it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
m_distance_helper.transformed_distance_from_coordinates(
this->query_object, *it_node_point, cache_point_begin, cache_point_begin + dim);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
cache_point_begin += dim;
}
// Now that the queue is full, we can gain time by keeping track
// of the current worst distance to interrupt the distance computation earlier
FT worst_dist = this->queue.top().second;
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
m_distance_helper.interruptible_transformed_distance(
this->query_object, *it_node_point, cache_point_begin, cache_point_begin + dim, worst_dist);
if (distance_to_query_object < worst_dist)
{
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
worst_dist = this->queue.top().second;
}
cache_point_begin += dim;
}
}
// Without cache
void search_nearest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_false)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
// As long as the queue is not full, the node is just inserted
for (; !this->queue.full() && it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
this->distance_instance.transformed_distance(this->query_object, *it_node_point);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
}
// Now that the queue is full, we can gain time by keeping track
// of the current worst distance to interrupt the distance computation earlier
FT worst_dist = this->queue.top().second;
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
m_distance_helper.interruptible_transformed_distance(
this->query_object, *it_node_point, worst_dist);
if (distance_to_query_object < worst_dist)
{
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
worst_dist = this->queue.top().second;
}
}
}
// With cache
void search_furthest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_true)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
int dim = m_tree.dim();
typename std::vector<FT>::const_iterator cache_point_begin = m_tree.cache_begin() + dim*(it_node_point - m_tree.begin());
// In furthest search mode, the interruptible distance cannot be used to optimize
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
m_distance_helper.transformed_distance_from_coordinates(
this->query_object, *it_node_point, cache_point_begin, cache_point_begin + dim);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
cache_point_begin += dim;
}
}
// Without cache
void search_furthest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_false)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
// In furthest search mode, the interruptible distance cannot be used to optimize
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
this->distance_instance.transformed_distance(this->query_object, *it_node_point);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
}
}
}; // class
} // namespace CGAL

61
Spatial_searching/include/CGAL/Kd_tree.h
View File

@ -8,7 +8,8 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>),
// : Waqar Khan <wkhan@mpi-inf.mpg.de>
// : Waqar Khan <wkhan@mpi-inf.mpg.de>,
// Clement Jamin (clement.jamin.pro@gmail.com)
#ifndef CGAL_KD_TREE_H
#define CGAL_KD_TREE_H
@ -37,7 +38,11 @@
namespace CGAL {
//template <class SearchTraits, class Splitter_=Median_of_rectangle<SearchTraits>, class UseExtendedNode = Tag_true >
template <class SearchTraits, class Splitter_=Sliding_midpoint<SearchTraits>, class UseExtendedNode = Tag_true >
template <
class SearchTraits,
class Splitter_=Sliding_midpoint<SearchTraits>,
class UseExtendedNode = Tag_true,
class EnablePointsCache = Tag_false>
class Kd_tree {
public:
@ -47,11 +52,11 @@ public:
typedef typename Splitter::Container Point_container;
typedef typename SearchTraits::FT FT;
typedef Kd_tree_node<SearchTraits, Splitter, UseExtendedNode > Node;
typedef Kd_tree_leaf_node<SearchTraits, Splitter, UseExtendedNode > Leaf_node;
typedef Kd_tree_internal_node<SearchTraits, Splitter, UseExtendedNode > Internal_node;
typedef Kd_tree<SearchTraits, Splitter> Tree;
typedef Kd_tree<SearchTraits, Splitter,UseExtendedNode> Self;
typedef Kd_tree_node<SearchTraits, Splitter, UseExtendedNode, EnablePointsCache> Node;
typedef Kd_tree_leaf_node<SearchTraits, Splitter, UseExtendedNode, EnablePointsCache> Leaf_node;
typedef Kd_tree_internal_node<SearchTraits, Splitter, UseExtendedNode, EnablePointsCache> Internal_node;
typedef Kd_tree<SearchTraits, Splitter, UseExtendedNode, EnablePointsCache> Tree;
typedef Kd_tree<SearchTraits, Splitter, UseExtendedNode, EnablePointsCache> Self;
typedef Node* Node_handle;
typedef const Node* Node_const_handle;
@ -69,6 +74,8 @@ public:
typedef typename internal::Get_dimension_tag<SearchTraits>::Dimension D;
typedef EnablePointsCache Enable_points_cache;
private:
SearchTraits traits_;
Splitter split;
@ -88,12 +95,18 @@ private:
Kd_tree_rectangle<FT,D>* bbox;
std::vector<Point_d> pts;
// Store a contiguous copy of the point coordinates
// for faster queries (reduce the number of cache misses)
std::vector<FT> points_cache;
// Instead of storing the points in arrays in the Kd_tree_node
// we put all the data in a vector in the Kd_tree.
// and we only store an iterator range in the Kd_tree_node.
//
std::vector<const Point_d*> data;
// Dimension of the points
int dim_;
#ifdef CGAL_HAS_THREADS
mutable CGAL_MUTEX building_mutex;//mutex used to protect const calls inducing build()
@ -103,7 +116,7 @@ private:
// protected copy constructor
Kd_tree(const Tree& tree)
: traits_(tree.traits_),built_(tree.built_)
: traits_(tree.traits_),built_(tree.built_),dim_(-1)
{};
@ -258,13 +271,13 @@ public:
CGAL_assertion(!removed_);
const Point_d& p = *pts.begin();
typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits_.construct_cartesian_const_iterator_d_object();
int dim = static_cast<int>(std::distance(ccci(p), ccci(p,0)));
dim_ = static_cast<int>(std::distance(ccci(p), ccci(p,0)));
data.reserve(pts.size());
for(unsigned int i = 0; i < pts.size(); i++){
data.push_back(&pts[i]);
}
Point_container c(dim, data.begin(), data.end(),traits_);
Point_container c(dim_, data.begin(), data.end(),traits_);
bbox = new Kd_tree_rectangle<FT,D>(c.bounding_box());
if (c.size() <= split.bucket_size()){
tree_root = create_leaf_node(c);
@ -275,9 +288,18 @@ public:
//Reorder vector for spatial locality
std::vector<Point_d> ptstmp;
ptstmp.resize(pts.size());
for (std::size_t i = 0; i < pts.size(); ++i){
for (std::size_t i = 0; i < pts.size(); ++i)
ptstmp[i] = *data[i];
// Cache?
if (Enable_points_cache::value)
{
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it = traits_.construct_cartesian_const_iterator_d_object();
points_cache.reserve(dim_ * pts.size());
for (std::size_t i = 0; i < pts.size(); ++i)
points_cache.insert(points_cache.end(), construct_it(ptstmp[i]), construct_it(ptstmp[i], 0));
}
for(std::size_t i = 0; i < leaf_nodes.size(); ++i){
std::ptrdiff_t tmp = leaf_nodes[i].begin() - pts.begin();
leaf_nodes[i].data = ptstmp.begin() + tmp;
@ -289,6 +311,12 @@ public:
built_ = true;
}
// Only correct when build() has been called
int dim() const
{
return dim_;
}
private:
//any call to this function is for the moment not threadsafe
void const_build() const {
@ -463,7 +491,7 @@ public:
const_build();
}
Kd_tree_rectangle<FT,D> b(*bbox);
return tree_root->search(it,q,b);
return tree_root->search(it,q,b,begin(),cache_begin(),dim_);
}
return it;
}
@ -479,7 +507,7 @@ public:
const_build();
}
Kd_tree_rectangle<FT,D> b(*bbox);
return tree_root->search_any_point(q,b);
return tree_root->search_any_point(q,b,begin(),cache_begin(),dim_);
}
return boost::none;
}
@ -538,6 +566,13 @@ public:
return *bbox;
}
typename std::vector<FT>::const_iterator
cache_begin() const
{
return points_cache.begin();
}
const_iterator
begin() const
{

225
Spatial_searching/include/CGAL/Kd_tree_node.h
View File

@ -19,31 +19,37 @@
#include <CGAL/Splitters.h>
#include <CGAL/Compact_container.h>
#include <CGAL/Has_member.h>
#include <CGAL/internal/Search_helpers.h>
#include <boost/cstdint.hpp>
namespace CGAL {
template <class SearchTraits, class Splitter, class UseExtendedNode>
CGAL_GENERATE_MEMBER_DETECTOR(contains_point_given_as_coordinates);
template <class SearchTraits, class Splitter, class UseExtendedNode, class EnablePointsCache>
class Kd_tree;
template < class TreeTraits, class Splitter, class UseExtendedNode >
template < class TreeTraits, class Splitter, class UseExtendedNode, class EnablePointsCache >
class Kd_tree_node {
friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
friend class Kd_tree<TreeTraits, Splitter, UseExtendedNode, EnablePointsCache>;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_handle Internal_node_handle;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Internal_node_const_handle Internal_node_const_handle;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_handle Leaf_node_handle;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Leaf_node_const_handle Leaf_node_const_handle;
typedef Kd_tree<TreeTraits, Splitter, UseExtendedNode, EnablePointsCache> Kdt;
typedef typename Kdt::Node_handle Node_handle;
typedef typename Kdt::Node_const_handle Node_const_handle;
typedef typename Kdt::Internal_node_handle Internal_node_handle;
typedef typename Kdt::Internal_node_const_handle Internal_node_const_handle;
typedef typename Kdt::Leaf_node_handle Leaf_node_handle;
typedef typename Kdt::Leaf_node_const_handle Leaf_node_const_handle;
typedef typename TreeTraits::Point_d Point_d;
typedef typename TreeTraits::FT FT;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Point_d_iterator Point_d_iterator;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::D D;
typedef typename Kdt::Separator Separator;
typedef typename Kdt::Point_d_iterator Point_d_iterator;
typedef typename Kdt::iterator iterator;
typedef typename Kdt::D D;
bool leaf;
@ -183,15 +189,19 @@ namespace CGAL {
template <class OutputIterator, class FuzzyQueryItem>
OutputIterator
search(OutputIterator it, const FuzzyQueryItem& q,
Kd_tree_rectangle<FT,D>& b) const
Kd_tree_rectangle<FT,D>& b,
typename Kdt::const_iterator tree_points_begin,
typename std::vector<FT>::const_iterator cache_begin,
int dim) const
{
if (is_leaf()) {
Leaf_node_const_handle node =
static_cast<Leaf_node_const_handle>(this);
if (node->size()>0)
for (iterator i=node->begin(); i != node->end(); i++)
if (q.contains(*i))
{*it++=*i;}
if (node->size() > 0)
{
typename internal::Has_points_cache<Kdt, internal::has_Enable_points_cache<Kdt>::type::value>::type dummy;
it = search_in_leaf(node, q, tree_points_begin, cache_begin, dim, it, dummy);
}
}
else {
Internal_node_const_handle node =
@ -204,12 +214,12 @@ namespace CGAL {
it=node->lower()->tree_items(it);
else
if (q.inner_range_intersects(b))
it=node->lower()->search(it,q,b);
it=node->lower()->search(it,q,b,tree_points_begin,cache_begin,dim);
if (q.outer_range_contains(b_upper))
it=node->upper()->tree_items(it);
else
if (q.inner_range_intersects(b_upper))
it=node->upper()->search(it,q,b_upper);
it=node->upper()->search(it,q,b_upper,tree_points_begin,cache_begin,dim);
};
return it;
}
@ -218,16 +228,20 @@ namespace CGAL {
template <class FuzzyQueryItem>
boost::optional<Point_d>
search_any_point(const FuzzyQueryItem& q,
Kd_tree_rectangle<FT,D>& b) const
Kd_tree_rectangle<FT,D>& b,
typename Kdt::const_iterator tree_points_begin,
typename std::vector<FT>::const_iterator cache_begin,
int dim) const
{
boost::optional<Point_d> result = boost::none;
if (is_leaf()) {
Leaf_node_const_handle node =
static_cast<Leaf_node_const_handle>(this);
if (node->size()>0)
for (iterator i=node->begin(); i != node->end(); i++)
if (q.contains(*i))
{ result = *i; break; }
{
typename internal::Has_points_cache<Kdt, internal::has_Enable_points_cache<Kdt>::type::value>::type dummy;
result = search_any_point_in_leaf(node, q, tree_points_begin, cache_begin, dim, dummy);
}
}
else {
Internal_node_const_handle node =
@ -237,26 +251,138 @@ namespace CGAL {
node->split_bbox(b, b_upper);
if (q.inner_range_intersects(b)) {
result = node->lower()->search_any_point(q,b);
result = node->lower()->search_any_point(q,b,tree_points_begin,cache_begin,dim);
if(result)
return result;
}
if (q.inner_range_intersects(b_upper))
result = node->upper()->search_any_point(q,b_upper);
result = node->upper()->search_any_point(q,b_upper,tree_points_begin,cache_begin,dim);
}
return result;
}
private:
// If contains_point_given_as_coordinates does not exist in `FuzzyQueryItem`
template <typename FuzzyQueryItem>
bool contains(
const FuzzyQueryItem& q,
Point_d const& p,
typename std::vector<FT>::const_iterator /*it_coord_begin*/,
typename std::vector<FT>::const_iterator /*it_coord_end*/,
Tag_false /*has_contains_point_given_as_coordinates*/) const
{
return q.contains(p);
}
// ... or if it exists
template <typename FuzzyQueryItem>
bool contains(
const FuzzyQueryItem& q,
Point_d const& /*p*/,
typename std::vector<FT>::const_iterator it_coord_begin,
typename std::vector<FT>::const_iterator it_coord_end,
Tag_true /*has_contains_point_given_as_coordinates*/) const
{
return q.contains_point_given_as_coordinates(it_coord_begin, it_coord_end);
}
// With cache
template<class FuzzyQueryItem, class OutputIterator>
OutputIterator search_in_leaf(
Leaf_node_const_handle node,
const FuzzyQueryItem &q,
typename Kdt::const_iterator tree_points_begin,
typename std::vector<FT>::const_iterator cache_begin,
int dim,
OutputIterator oit,
Tag_true /*has_points_cache*/) const
{
typename Kdt::iterator it_node_point = node->begin(), it_node_point_end = node->end();
typename std::vector<FT>::const_iterator cache_point_it = cache_begin + dim*(it_node_point - tree_points_begin);
for (; it_node_point != it_node_point_end; ++it_node_point, cache_point_it += dim)
{
Boolean_tag<has_contains_point_given_as_coordinates<FuzzyQueryItem>::value> dummy;
if (contains(q, *it_node_point, cache_point_it, cache_point_it + dim, dummy))
*oit++ = *it_node_point;
}
return oit;
}
// Without cache
template<class FuzzyQueryItem, class OutputIterator>
OutputIterator search_in_leaf(
Leaf_node_const_handle node,
const FuzzyQueryItem &q,
typename Kdt::const_iterator /*tree_points_begin*/,
typename std::vector<FT>::const_iterator /*cache_begin*/,
int /*dim*/,
OutputIterator oit,
Tag_false /*has_points_cache*/) const
{
for (iterator i = node->begin(); i != node->end(); ++i)
{
if (q.contains(*i))
*oit++ = *i;
}
return oit;
}
// With cache
template<class FuzzyQueryItem>
boost::optional<Point_d> search_any_point_in_leaf(
Leaf_node_const_handle node,
const FuzzyQueryItem &q,
typename Kdt::const_iterator tree_points_begin,
typename std::vector<FT>::const_iterator cache_begin,
int dim,
Tag_true /*has_points_cache*/) const
{
boost::optional<Point_d> result = boost::none;
typename Kdt::iterator it_node_point = node->begin(), it_node_point_end = node->end();
typename std::vector<FT>::const_iterator cache_point_it = cache_begin + dim*(it_node_point - tree_points_begin);
for (; it_node_point != it_node_point_end; ++it_node_point, cache_point_it += dim)
{
Boolean_tag<has_contains_point_given_as_coordinates<FuzzyQueryItem>::value> dummy;
if (contains(q, *it_node_point, cache_point_it, cache_point_it + dim, dummy))
{
result = *it_node_point;
break;
}
}
return result;
}
// Without cache
template<class FuzzyQueryItem>
boost::optional<Point_d> search_any_point_in_leaf(
Leaf_node_const_handle node,
const FuzzyQueryItem &q,
typename Kdt::const_iterator /*tree_points_begin*/,
typename std::vector<FT>::const_iterator /*cache_begin*/,
int /*dim*/,
Tag_false /*has_points_cache*/) const
{
boost::optional<Point_d> result = boost::none;
for (iterator i = node->begin(); i != node->end(); ++i)
{
if (q.contains(*i))
{
result = *i;
break;
}
}
return result;
}
};
template < class TreeTraits, class Splitter, class UseExtendedNode >
class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{
template < class TreeTraits, class Splitter, class UseExtendedNode, class EnablePointsCache >
class Kd_tree_leaf_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache >{
friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
friend class Kd_tree<TreeTraits, Splitter, UseExtendedNode, EnablePointsCache>;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::iterator iterator;
typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base;
typedef typename Kd_tree<TreeTraits, Splitter, UseExtendedNode, EnablePointsCache>::iterator iterator;
typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache> Base;
typedef typename TreeTraits::Point_d Point_d;
private:
@ -315,18 +441,20 @@ namespace CGAL {
template < class TreeTraits, class Splitter, class UseExtendedNode>
class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode >{
template < class TreeTraits, class Splitter, class UseExtendedNode, class EnablePointsCache>
class Kd_tree_internal_node : public Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache >{
friend class Kd_tree<TreeTraits,Splitter,UseExtendedNode>;
friend class Kd_tree<TreeTraits, Splitter, UseExtendedNode, EnablePointsCache>;
typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode> Base;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_handle Node_handle;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Node_const_handle Node_const_handle;
typedef Kd_tree<TreeTraits, Splitter, UseExtendedNode, EnablePointsCache> Kdt;
typedef Kd_tree_node< TreeTraits, Splitter, UseExtendedNode, EnablePointsCache> Base;
typedef typename Kdt::Node_handle Node_handle;
typedef typename Kdt::Node_const_handle Node_const_handle;
typedef typename TreeTraits::FT FT;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::Separator Separator;
typedef typename Kd_tree<TreeTraits,Splitter,UseExtendedNode>::D D;
typedef typename Kdt::Separator Separator;
typedef typename Kdt::D D;
private:
@ -463,18 +591,21 @@ namespace CGAL {
}
};//internal node
template < class TreeTraits, class Splitter>
class Kd_tree_internal_node<TreeTraits,Splitter,Tag_false> : public Kd_tree_node< TreeTraits, Splitter, Tag_false >{
template < class TreeTraits, class Splitter, class EnablePointsCache>
class Kd_tree_internal_node<TreeTraits,Splitter,Tag_false,EnablePointsCache>
: public Kd_tree_node< TreeTraits, Splitter, Tag_false, EnablePointsCache >
{
friend class Kd_tree<TreeTraits, Splitter, Tag_false, EnablePointsCache>;
friend class Kd_tree<TreeTraits,Splitter,Tag_false>;
typedef Kd_tree<TreeTraits, Splitter, Tag_false, EnablePointsCache> Kdt;
typedef Kd_tree_node< TreeTraits, Splitter, Tag_false> Base;
typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_handle Node_handle;
typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Node_const_handle Node_const_handle;
typedef Kd_tree_node< TreeTraits, Splitter, Tag_false, EnablePointsCache> Base;
typedef typename Kdt::Node_handle Node_handle;
typedef typename Kdt::Node_const_handle Node_const_handle;
typedef typename TreeTraits::FT FT;
typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::Separator Separator;
typedef typename Kd_tree<TreeTraits,Splitter,Tag_false>::D D;
typedef typename Kdt::Separator Separator;
typedef typename Kdt::D D;
private:

174
Spatial_searching/include/CGAL/Orthogonal_incremental_neighbor_search.h
View File

@ -9,6 +9,7 @@
//
//
// Author(s) : Hans Tangelder (<hanst@cs.uu.nl>)
// Clement Jamin (clement.jamin.pro@gmail.com)
#ifndef CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH
#define CGAL_ORTHOGONAL_INCREMENTAL_NEIGHBOR_SEARCH
@ -16,21 +17,23 @@
#include <CGAL/license/Spatial_searching.h>
#include <CGAL/disable_warnings.h>
#include <CGAL/Kd_tree.h>
#include <CGAL/Euclidean_distance.h>
#include <CGAL/tuple.h>
#include <CGAL/internal/Search_helpers.h>
#include <cstring>
#include <list>
#include <queue>
#include <memory>
#include <CGAL/Kd_tree.h>
#include <CGAL/Euclidean_distance.h>
#include <CGAL/tuple.h>
#include <iterator> // for std::distance
namespace CGAL {
template <class SearchTraits,
class Distance_= typename internal::Spatial_searching_default_distance<SearchTraits>::type,
class Splitter_ = Sliding_midpoint<SearchTraits>,
class Tree_= Kd_tree<SearchTraits, Splitter_, Tag_true> >
class Tree_= Kd_tree<SearchTraits, Splitter_, Tag_true, Tag_false> >
class Orthogonal_incremental_neighbor_search {
public:
@ -72,7 +75,8 @@ namespace CGAL {
Distance_vector dists;
Distance Orthogonal_distance_instance;
Distance orthogonal_distance_instance;
internal::Distance_helper<Distance, SearchTraits> m_distance_helper;
FT multiplication_factor;
@ -84,6 +88,8 @@ namespace CGAL {
FT rd;
int m_dim;
Tree const& m_tree;
class Priority_higher {
public:
@ -139,27 +145,30 @@ namespace CGAL {
FT Eps=FT(0.0), bool search_nearest=true)
: traits(tree.traits()),number_of_neighbours_computed(0), number_of_internal_nodes_visited(0),
number_of_leaf_nodes_visited(0), number_of_items_visited(0),
Orthogonal_distance_instance(tr), multiplication_factor(Orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)),
orthogonal_distance_instance(tr),
m_distance_helper(orthogonal_distance_instance, traits),
multiplication_factor(orthogonal_distance_instance.transformed_distance(FT(1.0)+Eps)),
query_point(q), search_nearest_neighbour(search_nearest),
m_tree(tree),
PriorityQueue(Priority_higher(search_nearest)), Item_PriorityQueue(Distance_smaller(search_nearest)),
reference_count(1)
{
if (tree.empty()) return;
if (m_tree.empty()) return;
typename SearchTraits::Construct_cartesian_const_iterator_d ccci=traits.construct_cartesian_const_iterator_d_object();
int dim = static_cast<int>(std::distance(ccci(q), ccci(q,0)));
m_dim = static_cast<int>(std::distance(ccci(q), ccci(q,0)));
dists.resize(dim);
for(int i=0 ; i<dim ; ++i){
dists.resize(m_dim);
for(int i=0 ; i<m_dim ; ++i){
dists[i] = 0;
}
if (search_nearest){
distance_to_root=
Orthogonal_distance_instance.min_distance_to_rectangle(q, tree.bounding_box(),dists);
Node_with_distance *The_Root = new Node_with_distance(tree.root(),
orthogonal_distance_instance.min_distance_to_rectangle(q, m_tree.bounding_box(),dists);
Node_with_distance *The_Root = new Node_with_distance(m_tree.root(),
distance_to_root, dists);
PriorityQueue.push(The_Root);
@ -169,9 +178,9 @@ namespace CGAL {
}
else{
distance_to_root=
Orthogonal_distance_instance.max_distance_to_rectangle(q,
tree.bounding_box(), dists);
Node_with_distance *The_Root = new Node_with_distance(tree.root(),
orthogonal_distance_instance.max_distance_to_rectangle(q,
m_tree.bounding_box(), dists);
Node_with_distance *The_Root = new Node_with_distance(m_tree.root(),
distance_to_root, dists);
PriorityQueue.push(The_Root);
@ -253,6 +262,83 @@ namespace CGAL {
delete The_item_top;
}
// With cache
bool search_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_true, bool search_furthest)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
typename std::vector<FT>::const_iterator cache_point_begin = m_tree.cache_begin() + m_dim*(it_node_point - m_tree.begin());
for (; it_node_point != it_node_point_end; ++it_node_point)
{
number_of_items_visited++;
FT distance_to_query_point =
m_distance_helper.transformed_distance_from_coordinates(
query_point, *it_node_point, cache_point_begin, cache_point_begin + m_dim);
Point_with_transformed_distance *NN_Candidate =
new Point_with_transformed_distance(*it_node_point, distance_to_query_point);
Item_PriorityQueue.push(NN_Candidate);
cache_point_begin += m_dim;
}
// old top of PriorityQueue has been processed,
// hence update rd
bool next_neighbour_found;
if (!(PriorityQueue.empty()))
{
rd = CGAL::cpp11::get<1>(*PriorityQueue.top());
next_neighbour_found = (search_furthest ?
multiplication_factor*rd < Item_PriorityQueue.top()->second
: multiplication_factor*rd > Item_PriorityQueue.top()->second);
}
else // priority queue empty => last neighbour found
{
next_neighbour_found = true;
}
number_of_neighbours_computed++;
return next_neighbour_found;
}
// Without cache
bool search_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_false, bool search_furthest)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
for (; it_node_point != it_node_point_end; ++it_node_point)
{
number_of_items_visited++;
FT distance_to_query_point=
orthogonal_distance_instance.transformed_distance(query_point, *it_node_point);
Point_with_transformed_distance *NN_Candidate =
new Point_with_transformed_distance(*it_node_point, distance_to_query_point);
Item_PriorityQueue.push(NN_Candidate);
}
// old top of PriorityQueue has been processed,
// hence update rd
bool next_neighbour_found;
if (!(PriorityQueue.empty()))
{
rd = CGAL::cpp11::get<1>(*PriorityQueue.top());
next_neighbour_found = (search_furthest ?
multiplication_factor*rd < Item_PriorityQueue.top()->second
: multiplication_factor*rd > Item_PriorityQueue.top()->second);
}
else // priority queue empty => last neighbour found
{
next_neighbour_found=true;
}
number_of_neighbours_computed++;
return next_neighbour_found;
}
void
Compute_the_next_nearest_neighbour()
{
@ -283,7 +369,7 @@ namespace CGAL {
FT diff2 = val - node->lower_high_value();
if (diff1 + diff2 < FT(0.0)) {
new_rd=
Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
CGAL_assertion(new_rd >= copy_rd);
dists[new_cut_dim] = diff1;
@ -295,7 +381,7 @@ namespace CGAL {
}
else { // compute new distance
new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
new_rd=orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
CGAL_assertion(new_rd >= copy_rd);
dists[new_cut_dim] = diff2;
Node_with_distance *Lower_Child =
@ -310,29 +396,8 @@ namespace CGAL {
static_cast<typename Tree::Leaf_node_const_handle>(N);
number_of_leaf_nodes_visited++;
if (node->size() > 0) {
for (typename Tree::iterator it=node->begin(); it != node->end(); it++) {
number_of_items_visited++;
FT distance_to_query_point=
Orthogonal_distance_instance.transformed_distance(query_point,*it);
Point_with_transformed_distance *NN_Candidate=
new Point_with_transformed_distance(*it,distance_to_query_point);
Item_PriorityQueue.push(NN_Candidate);
}
// old top of PriorityQueue has been processed,
// hence update rd
if (!(PriorityQueue.empty())) {
rd = std::get<1>(*PriorityQueue.top());
next_neighbour_found =
(multiplication_factor*rd >
Item_PriorityQueue.top()->second);
}
else // priority queue empty => last neighbour found
{
next_neighbour_found=true;
}
number_of_neighbours_computed++;
typename internal::Has_points_cache<Tree, internal::has_Enable_points_cache<Tree>::type::value>::type dummy;
next_neighbour_found = search_in_leaf(node, dummy, false);
}
} // next_neighbour_found or priority queue is empty
// in the latter case also the item priority quee is empty
@ -370,7 +435,7 @@ namespace CGAL {
if (diff1 + diff2 < FT(0.0)) {
diff1 = val - node->upper_high_value();
new_rd=
Orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
orthogonal_distance_instance.new_distance(copy_rd,dst,diff1,new_cut_dim);
Node_with_distance *Lower_Child =
new Node_with_distance(node->lower(), copy_rd, dists);
PriorityQueue.push(Lower_Child);
@ -381,7 +446,7 @@ namespace CGAL {
}
else { // compute new distance
diff2 = val - node->lower_low_value();
new_rd=Orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
new_rd=orthogonal_distance_instance.new_distance(copy_rd,dst,diff2,new_cut_dim);
Node_with_distance *Upper_Child =
new Node_with_distance(node->upper(), copy_rd, dists);
PriorityQueue.push(Upper_Child);
@ -395,29 +460,8 @@ namespace CGAL {
static_cast<typename Tree::Leaf_node_const_handle>(N);
number_of_leaf_nodes_visited++;
if (node->size() > 0) {
for (typename Tree::iterator it=node->begin(); it != node->end(); it++) {
number_of_items_visited++;
FT distance_to_query_point=
Orthogonal_distance_instance.transformed_distance(query_point,*it);
Point_with_transformed_distance *NN_Candidate=
new Point_with_transformed_distance(*it,distance_to_query_point);
Item_PriorityQueue.push(NN_Candidate);
}
// old top of PriorityQueue has been processed,
// hence update rd
if (!(PriorityQueue.empty())) {
rd = std::get<1>(*PriorityQueue.top());
next_neighbour_found =
(multiplication_factor*rd <
Item_PriorityQueue.top()->second);
}
else // priority queue empty => last neighbour found
{
next_neighbour_found=true;
}
number_of_neighbours_computed++;
typename internal::Has_points_cache<Tree, internal::has_Enable_points_cache<Tree>::type::value>::type dummy;
next_neighbour_found = search_in_leaf(node, dummy, true);
}
} // next_neighbour_found or priority queue is empty
// in the latter case also the item priority quee is empty

209
Spatial_searching/include/CGAL/Orthogonal_k_neighbor_search.h
View File

@ -8,7 +8,9 @@
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Gael Guennebaud (gael.guennebaud@inria.fr), Hans Tangelder (<hanst@cs.uu.nl>)
// Author(s) : Gael Guennebaud (gael.guennebaud@inria.fr),
// Hans Tangelder (<hanst@cs.uu.nl>),
// Clement Jamin (clement.jamin.pro@gmail.com)
#ifndef CGAL_ORTHOGONAL_K_NEIGHBOR_SEARCH_H
#define CGAL_ORTHOGONAL_K_NEIGHBOR_SEARCH_H
@ -18,36 +20,49 @@
#include <CGAL/disable_warnings.h>
#include <CGAL/internal/K_neighbor_search.h>
#include <CGAL/internal/Search_helpers.h>
#include <iterator> // for std::distance
namespace CGAL {
template <class SearchTraits,
class Distance= typename internal::Spatial_searching_default_distance<SearchTraits>::type,
class Splitter= Sliding_midpoint<SearchTraits> ,
class Tree= Kd_tree<SearchTraits, Splitter, Tag_true> >
class Orthogonal_k_neighbor_search: public internal::K_neighbor_search<SearchTraits,Distance,Splitter,Tree> {
class Tree= Kd_tree<SearchTraits, Splitter, Tag_true, Tag_false> >
class Orthogonal_k_neighbor_search: public internal::K_neighbor_search<SearchTraits,Distance,Splitter,Tree>
{
typedef internal::K_neighbor_search<SearchTraits,Distance,Splitter,Tree> Base;
typedef typename Tree::Point_d Point;
typename SearchTraits::Cartesian_const_iterator_d query_object_it;
std::vector<typename Base::FT> dists;
public:
typedef typename Base::FT FT;
private:
typename SearchTraits::Cartesian_const_iterator_d query_object_it;
internal::Distance_helper<Distance, SearchTraits> m_distance_helper;
std::vector<FT> dists;
int m_dim;
Tree const& m_tree;
public:
Orthogonal_k_neighbor_search(const Tree& tree, const typename Base::Query_item& q,
unsigned int k=1, FT Eps=FT(0.0), bool Search_nearest=true, const Distance& d=Distance(),bool sorted=true)
: Base(q,k,Eps,Search_nearest,d)
: Base(q,k,Eps,Search_nearest,d),
m_distance_helper(this->distance_instance, tree.traits()),
m_tree(tree)
{
if (tree.empty()) return;
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it=tree.traits().construct_cartesian_const_iterator_d_object();
query_object_it = construct_it(this->query_object);
int dim = static_cast<int>(std::distance(query_object_it, construct_it(this->query_object,0)));
m_dim = static_cast<int>(std::distance(query_object_it, construct_it(this->query_object,0)));
dists.resize(dim);
for(int i=0;i<dim;i++)
dists.resize(m_dim);
for(int i=0;i<m_dim;i++)
dists[i]=0;
FT distance_to_root;
@ -60,28 +75,106 @@ public:
compute_furthest_neighbors_orthogonally(tree.root(), distance_to_root);
}
if (sorted) this->queue.sort();
}
private:
// With cache
void search_nearest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_true)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
typename std::vector<FT>::const_iterator cache_point_begin = m_tree.cache_begin() + m_dim*(it_node_point - m_tree.begin());
// As long as the queue is not full, the node is just inserted
for (; !this->queue.full() && it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
m_distance_helper.transformed_distance_from_coordinates(
this->query_object, *it_node_point, cache_point_begin, cache_point_begin + m_dim);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
cache_point_begin += m_dim;
}
// Now that the queue is full, we can gain time by keeping track
// of the current worst distance to interrupt the distance computation earlier
FT worst_dist = this->queue.top().second;
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
m_distance_helper.interruptible_transformed_distance(
this->query_object, *it_node_point, cache_point_begin, cache_point_begin + m_dim, worst_dist);
if (distance_to_query_object < worst_dist)
{
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
worst_dist = this->queue.top().second;
}
cache_point_begin += m_dim;
}
}
// Without cache
void search_nearest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_false)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
// As long as the queue is not full, the node is just inserted
for (; !this->queue.full() && it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
this->distance_instance.transformed_distance(this->query_object, *it_node_point);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
}
// Now that the queue is full, we can gain time by keeping track
// of the current worst distance to interrupt the distance computation earlier
FT worst_dist = this->queue.top().second;
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
m_distance_helper.interruptible_transformed_distance(
this->query_object, *it_node_point, worst_dist);
if (distance_to_query_object < worst_dist)
{
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
worst_dist = this->queue.top().second;
}
}
}
void compute_nearest_neighbors_orthogonally(typename Base::Node_const_handle N, FT rd)
{
if (!(N->is_leaf()))
if (N->is_leaf())
{
// n is a leaf
typename Tree::Leaf_node_const_handle node =
static_cast<typename Tree::Leaf_node_const_handle>(N);
this->number_of_leaf_nodes_visited++;
if (node->size() > 0)
{
typename internal::Has_points_cache<Tree, internal::has_Enable_points_cache<Tree>::type::value>::type dummy;
search_nearest_in_leaf(node, dummy);
}
}
else
{
typename Tree::Internal_node_const_handle node =
static_cast<typename Tree::Internal_node_const_handle>(N);
this->number_of_internal_nodes_visited++;
int new_cut_dim=node->cutting_dimension();
int new_cut_dim = node->cutting_dimension();
typename Base::Node_const_handle bestChild, otherChild;
FT new_off;
FT val = *(query_object_it + new_cut_dim);
FT diff1 = val - node->upper_low_value();
FT diff2 = val - node->lower_high_value();
if ( (diff1 + diff2 < FT(0.0)) )
if ((diff1 + diff2 < FT(0.0)))
{
new_off = diff1;
bestChild = node->lower();
@ -89,46 +182,71 @@ private:
}
else // compute new distance
{
new_off= diff2;
new_off = diff2;
bestChild = node->upper();
otherChild = node->lower();
}
compute_nearest_neighbors_orthogonally(bestChild, rd);
FT dst=dists[new_cut_dim];
FT new_rd = this->distance_instance.new_distance(rd,dst,new_off,new_cut_dim);
dists[new_cut_dim]=new_off;
FT dst = dists[new_cut_dim];
FT new_rd = this->distance_instance.new_distance(rd, dst, new_off, new_cut_dim);
dists[new_cut_dim] = new_off;
if (this->branch_nearest(new_rd))
{
compute_nearest_neighbors_orthogonally(otherChild, new_rd);
}
dists[new_cut_dim]=dst;
dists[new_cut_dim] = dst;
}
else
}
// With cache
void search_furthest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_true)
{
// n is a leaf
typename Tree::Leaf_node_const_handle node =
static_cast<typename Tree::Leaf_node_const_handle>(N);
this->number_of_leaf_nodes_visited++;
bool full = this->queue.full();
FT worst_dist = this->queue.top().second;
if (node->size() > 0)
{
for (typename Tree::iterator it=node->begin(); it != node->end(); it++)
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
typename std::vector<FT>::const_iterator cache_point_begin = m_tree.cache_begin() + m_dim*(it_node_point - m_tree.begin());
// In furthest search mode, the interruptible distance cannot be used to optimize
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object=
this->distance_instance.transformed_distance(this->query_object,*it);
if(!full || distance_to_query_object < worst_dist)
this->queue.insert(std::make_pair(&(*it),distance_to_query_object));
FT distance_to_query_object =
m_distance_helper.transformed_distance_from_coordinates(
this->query_object, *it_node_point, cache_point_begin, cache_point_begin + m_dim);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
cache_point_begin += m_dim;
}
}
// Without cache
void search_furthest_in_leaf(typename Tree::Leaf_node_const_handle node, Tag_false)
{
typename Tree::iterator it_node_point = node->begin(), it_node_point_end = node->end();
// In furthest search mode, the interruptible distance cannot be used to optimize
for (; it_node_point != it_node_point_end; ++it_node_point)
{
this->number_of_items_visited++;
FT distance_to_query_object =
this->distance_instance.transformed_distance(this->query_object, *it_node_point);
this->queue.insert(std::make_pair(&(*it_node_point), distance_to_query_object));
}
}
void compute_furthest_neighbors_orthogonally(typename Base::Node_const_handle N, FT rd)
{
if (!(N->is_leaf()))
if (N->is_leaf())
{
// n is a leaf
typename Tree::Leaf_node_const_handle node =
static_cast<typename Tree::Leaf_node_const_handle>(N);
this->number_of_leaf_nodes_visited++;
if (node->size() > 0)
{
typename internal::Has_points_cache<Tree, internal::has_Enable_points_cache<Tree>::type::value>::type dummy;
search_furthest_in_leaf(node, dummy);
}
}
else
{
typename Tree::Internal_node_const_handle node =
static_cast<typename Tree::Internal_node_const_handle>(N);
@ -163,23 +281,6 @@ private:
compute_furthest_neighbors_orthogonally(otherChild, new_rd);
dists[new_cut_dim]=dst;
}
else
{
// n is a leaf
typename Tree::Leaf_node_const_handle node =
static_cast<typename Tree::Leaf_node_const_handle>(N);
this->number_of_leaf_nodes_visited++;
if (node->size() > 0)
{
for (typename Tree::iterator it=node->begin(); it != node->end(); it++)
{
this->number_of_items_visited++;
FT distance_to_query_object=
this->distance_instance.transformed_distance(this->query_object,*it);
this->queue.insert(std::make_pair(&(*it),distance_to_query_object));
}
}
}
}
}; // class

2
Spatial_searching/include/CGAL/Search_traits_adapter.h
View File

@ -269,7 +269,6 @@ public:
template <class Point_with_info,class PointPropertyMap,class Base_distance>
class Distance_adapter : public Base_distance {
PointPropertyMap ppmap;
typedef typename Base_distance::FT FT;
public:
@ -279,6 +278,7 @@ public:
using Base_distance::transformed_distance;
typedef typename Base_distance::FT FT;
typedef Point_with_info Point_d;
typedef typename Base_distance::Query_item Query_item;

240
Spatial_searching/include/CGAL/internal/Search_helpers.h Normal file
View File

@ -0,0 +1,240 @@
// Copyright (c) 2002,2011 Utrecht University (The Netherlands).
// All rights reserved.
//
// This file is part of CGAL (www.cgal.org).
//
// $URL$
// $Id$
// SPDX-License-Identifier: GPL-3.0-or-later OR LicenseRef-Commercial
//
//
// Author(s) : Clement Jamin (clement.jamin.pro@gmail.com)
#ifndef CGAL_INTERNAL_SEARCH_HELPERS_H
#define CGAL_INTERNAL_SEARCH_HELPERS_H
#include <CGAL/license/Spatial_searching.h>
#include <CGAL/Has_member.h>
#include <boost/mpl/has_xxx.hpp>
namespace CGAL {
namespace internal {
// Helper struct to know at compile-time if there is a cache of the points
// stored in the tree
template <typename Tree, bool has_enable_points_cache>
struct Has_points_cache;
template <typename Tree>
struct Has_points_cache<Tree, true>
{
typedef typename Tree::Enable_points_cache type;
static const bool value = type::value;
};
template <typename Tree>
struct Has_points_cache<Tree, false>
{
typedef Tag_false type;
static const bool value = false;
};
CGAL_GENERATE_MEMBER_DETECTOR(transformed_distance_from_coordinates);
CGAL_GENERATE_MEMBER_DETECTOR(interruptible_transformed_distance);
BOOST_MPL_HAS_XXX_TRAIT_NAMED_DEF(has_Enable_points_cache, Enable_points_cache, false)
// If transformed_distance_from_coordinates does not exist in `Distance`
template <typename Distance, typename SearchTraits, bool has_transformed_distance_from_coordinates>
class Transformed_distance_from_coordinates
{
typedef typename Distance::FT FT;
typedef typename Distance::Point_d Point;
typedef typename Distance::Query_item Query_item;
public:
Transformed_distance_from_coordinates(Distance const& distance)
: m_distance(distance)
{}
FT operator() (
const Query_item& q,
Point const& p,
typename std::vector<FT>::const_iterator /*it_coord_begin*/,
typename std::vector<FT>::const_iterator /*it_coord_end*/) const
{
return m_distance.transformed_distance(q, p);
}
private:
Distance const& m_distance;
};
// ... or if it exists
template <typename Distance, typename SearchTraits>
class Transformed_distance_from_coordinates<Distance, SearchTraits, true>
{
typedef typename Distance::FT FT;
typedef typename Distance::Point_d Point;
typedef typename Distance::Query_item Query_item;
public:
Transformed_distance_from_coordinates(Distance const& distance)
: m_distance(distance)
{}
FT operator() (
const Query_item& q,
Point const& /*p*/,
typename std::vector<FT>::const_iterator it_coord_begin,
typename std::vector<FT>::const_iterator it_coord_end) const
{
return m_distance.transformed_distance_from_coordinates(q, it_coord_begin, it_coord_end);
}
private:
Distance const& m_distance;
};
// If interruptible_transformed_distance does not exist in `Distance`
template <typename Distance, typename SearchTraits, bool has_interruptible_transformed_distance>
class Interruptible_transformed_distance
{
typedef typename Distance::FT FT;
typedef typename Distance::Point_d Point;
typedef typename Distance::Query_item Query_item;
public:
typedef Transformed_distance_from_coordinates<Distance, SearchTraits,
has_transformed_distance_from_coordinates<Distance>::value> Tdfc;
Interruptible_transformed_distance(
SearchTraits const&, Distance const& distance, Tdfc const& tdfc)
: m_distance(distance), m_ref_to_tdfc(tdfc)
{}
FT operator() (
const Query_item& q,
Point const& p,
FT) const
{
return m_distance.transformed_distance(q, p);
}
FT operator() (
const Query_item& q,
Point const& p,
typename std::vector<FT>::const_iterator it_coord_begin,
typename std::vector<FT>::const_iterator it_coord_end,
FT) const
{
return m_ref_to_tdfc(q, p, it_coord_begin, it_coord_end);
}
private:
Distance const& m_distance;
Tdfc const& m_ref_to_tdfc;
};
// ... or if it exists
template <typename Distance, typename SearchTraits>
class Interruptible_transformed_distance<Distance, SearchTraits, true>
{
typedef typename Distance::FT FT;
typedef typename Distance::Point_d Point;
typedef typename Distance::Query_item Query_item;
public:
typedef Transformed_distance_from_coordinates<Distance, SearchTraits,
has_transformed_distance_from_coordinates<Distance>::value> Tdfc;
Interruptible_transformed_distance(
SearchTraits const& traits, Distance const& distance, Tdfc const&)
: m_traits(traits), m_distance(distance)
{}
FT operator() (
const Query_item& q,
Point const& p,
FT stop_if_geq_to_this) const
{
typename SearchTraits::Construct_cartesian_const_iterator_d construct_it =
m_traits.construct_cartesian_const_iterator_d_object();
return m_distance.interruptible_transformed_distance(
q, construct_it(p), construct_it(p, 0), stop_if_geq_to_this);
}
FT operator() (
const Query_item& q,
Point const& /*p*/,
typename std::vector<FT>::const_iterator it_coord_begin,
typename std::vector<FT>::const_iterator it_coord_end,
FT stop_if_geq_to_this) const
{
return m_distance.interruptible_transformed_distance(
q, it_coord_begin, it_coord_end, stop_if_geq_to_this);
}
private:
SearchTraits const& m_traits;
Distance const& m_distance;
};
template <typename Distance, typename SearchTraits>
class Distance_helper
{
typedef typename Distance::FT FT;
typedef typename Distance::Point_d Point;
typedef typename Distance::Query_item Query_item;
public:
Distance_helper(Distance const& distance, SearchTraits const& traits)
: m_distance(distance), m_tdfc(m_distance), m_itd(traits, distance, m_tdfc)
{}
FT
transformed_distance_from_coordinates(
const Query_item& q,
Point const& p,
typename std::vector<FT>::const_iterator it_coord_begin,
typename std::vector<FT>::const_iterator it_coord_end)
{
return m_tdfc(q, p, it_coord_begin, it_coord_end);
}
FT
interruptible_transformed_distance(
const Query_item& q,
Point const& p,
FT stop_if_geq_to_this)
{
return m_itd(q, p, stop_if_geq_to_this);
}
FT
interruptible_transformed_distance(
const Query_item& q,
Point const& p,
typename std::vector<FT>::const_iterator it_coord_begin,
typename std::vector<FT>::const_iterator it_coord_end,
FT stop_if_geq_to_this)
{
return m_itd(q, p, it_coord_begin, it_coord_end, stop_if_geq_to_this);
}
private:
Distance const& m_distance;
Transformed_distance_from_coordinates<Distance, SearchTraits,
has_transformed_distance_from_coordinates<Distance>::value> m_tdfc;
Interruptible_transformed_distance<Distance, SearchTraits,
has_interruptible_transformed_distance<Distance>::value> m_itd;
}; // Distance_helper
}} // namespace CGAL::internal
#endif // CGAL_INTERNAL_SEARCH_HELPERSS_H

1
Spatial_searching/test/Spatial_searching/Distance.h
View File

@ -1,4 +1,5 @@
struct Distance {
typedef Point Point_d;
typedef Point Query_item;
typedef double FT;