boost/geometry
2
0
Fork 0
mirror of https://github.com/boostorg/geometry.git synced 2026-02-12 00:02:09 +00:00
Code Issues Projects Releases Wiki Activity

Merge branch 'develop' into feature/gc3

Browse Source
This commit is contained in:
Adam Wulkiewicz
2021-09-22 23:29:00 +02:00
parent 4b72c60ecf 2c811b427c
commit 8c9b11bb85
82 changed files with 2867 additions and 1373 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
include/boost/geometry/algorithms/detail/buffer/piece_border.hpp
View File

@@ -24,7 +24,6 @@
#include <boost/geometry/algorithms/comparable_distance.hpp>
#include <boost/geometry/algorithms/equals.hpp>
#include <boost/geometry/algorithms/expand.hpp>
#include <boost/geometry/algorithms/detail/buffer/buffer_box.hpp>
#include <boost/geometry/algorithms/detail/buffer/buffer_policies.hpp>
#include <boost/geometry/algorithms/detail/expand_by_epsilon.hpp>
#include <boost/geometry/strategies/cartesian/turn_in_ring_winding.hpp>

15
include/boost/geometry/algorithms/detail/equals/collect_vectors.hpp
View File

@@ -39,7 +39,8 @@
#include <boost/geometry/views/detail/closed_clockwise_view.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/spherical/ssf.hpp>
#include <boost/geometry/strategies/normalize.hpp>
@@ -61,11 +62,11 @@ struct collected_vector
: nyi::not_implemented_tag
{};
// compatible with side_by_triangle cartesian strategy
// compatible with side_robust cartesian strategy
template <typename T, typename Geometry, typename CT, typename CSTag>
struct collected_vector
<
T, Geometry, strategy::side::side_by_triangle<CT>, CSTag
T, Geometry, strategy::side::side_robust<CT>, CSTag
>
{
typedef T type;
@@ -156,6 +157,14 @@ private:
//T dx_0, dy_0;
};
template <typename T, typename Geometry, typename CT, typename CSTag>
struct collected_vector
<
T, Geometry, strategy::side::side_by_triangle<CT>, CSTag
>
: collected_vector<T, Geometry, strategy::side::side_robust<CT>, CSTag>
{};
// Compatible with spherical_side_formula which currently
// is the default spherical_equatorial and geographic strategy
// so CSTag is spherical_equatorial_tag or geographic_tag

111
include/boost/geometry/core/coordinate_promotion.hpp Normal file
View File

@@ -0,0 +1,111 @@
// Boost.Geometry
// Copyright (c) 2021 Barend Gehrels, Amsterdam, the Netherlands.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_GEOMETRY_CORE_COORDINATE_PROMOTION_HPP
#define BOOST_GEOMETRY_CORE_COORDINATE_PROMOTION_HPP
#include <boost/geometry/core/coordinate_type.hpp>
// TODO: move this to a future headerfile implementing traits for these types
#include <boost/rational.hpp>
#include <boost/multiprecision/cpp_bin_float.hpp>
#include <boost/multiprecision/cpp_int.hpp>
namespace boost { namespace geometry
{
namespace traits
{
// todo
} // namespace traits
/*!
\brief Meta-function converting, if necessary, to "a floating point" type
\details
- if input type is integer, type is double
- else type is input type
\ingroup utility
*/
// TODO: replace with, or call, promoted_to_floating_point
template <typename T, typename PromoteIntegerTo = double>
struct promote_floating_point
{
typedef std::conditional_t
<
std::is_integral<T>::value,
PromoteIntegerTo,
T
> type;
};
// TODO: replace with promoted_to_floating_point
template <typename Geometry>
struct fp_coordinate_type
{
typedef typename promote_floating_point
<
typename coordinate_type<Geometry>::type
>::type type;
};
namespace detail
{
// Promote any integral type to double. Floating point
// and other user defined types stay as they are, unless specialized.
// TODO: we shold add a coordinate_promotion traits for promotion to
// floating point or (larger) integer types.
template <typename Type>
struct promoted_to_floating_point
{
using type = std::conditional_t
<
std::is_integral<Type>::value, double, Type
>;
};
// Boost.Rational goes to double
template <typename T>
struct promoted_to_floating_point<boost::rational<T>>
{
using type = double;
};
// Any Boost.Multiprecision goes to double (for example int128_t),
// unless specialized
template <typename Backend>
struct promoted_to_floating_point<boost::multiprecision::number<Backend>>
{
using type = double;
};
// Boost.Multiprecision binary floating point numbers are used as FP.
template <unsigned Digits>
struct promoted_to_floating_point
<
boost::multiprecision::number
<
boost::multiprecision::cpp_bin_float<Digits>
>
>
{
using type = boost::multiprecision::number
<
boost::multiprecision::cpp_bin_float<Digits>
>;
};
}
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_CORE_COORDINATE_PROMOTION_HPP

10
include/boost/geometry/core/coordinate_type.hpp
View File

@@ -22,7 +22,6 @@
#include <boost/geometry/core/point_type.hpp>
#include <boost/geometry/core/static_assert.hpp>
#include <boost/geometry/core/tag.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/type_traits_std.hpp>
@@ -94,15 +93,6 @@ struct coordinate_type
>::type type;
};
template <typename Geometry>
struct fp_coordinate_type
{
typedef typename promote_floating_point
<
typename coordinate_type<Geometry>::type
>::type type;
};
/*!
\brief assert_coordinate_type_equal, a compile-time check for equality of two coordinate types
\ingroup utility

3
include/boost/geometry/core/radian_access.hpp
View File

@@ -27,8 +27,7 @@
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/util/math.hpp>

7
include/boost/geometry/extensions/algorithms/distance_info.hpp
View File

@@ -5,8 +5,8 @@
// Copyright (c) 2009-2013 Mateusz Loskot, London, UK.
// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland.
// This file was modified by Oracle on 2020.
// Modifications copyright (c) 2020, Oracle and/or its affiliates.
// This file was modified by Oracle on 2020-2021.
// Modifications copyright (c) 2020-2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Use, modification and distribution is subject to the Boost Software License,
@@ -30,10 +30,9 @@
#include <boost/geometry/geometries/concepts/check.hpp>
#include <boost/geometry/strategies/distance.hpp>
#include <boost/geometry/strategies/default_distance_result.hpp>
#include <boost/geometry/strategies/distance/cartesian.hpp>
#include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
#include <boost/geometry/extensions/strategies/cartesian/distance_info.hpp>
#include <boost/geometry/util/math.hpp>

16
include/boost/geometry/extensions/gis/io/wkb/detail/endian.hpp
View File

@@ -2,8 +2,8 @@
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
// This file was modified by Oracle on 2020.
// Modifications copyright (c) 2020, Oracle and/or its affiliates.
// This file was modified by Oracle on 2020-2021.
// Modifications copyright (c) 2020-2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Use, modification and distribution is subject to the Boost Software License,
@@ -31,7 +31,7 @@
#include <boost/config.hpp>
#include <boost/cstdint.hpp>
#include <boost/static_assert.hpp>
#include <boost/detail/endian.hpp>
#include <boost/predef/other/endian.h>
#if CHAR_BIT != 8
#error Platforms with CHAR_BIT != 8 are not supported
@@ -50,10 +50,16 @@ namespace detail { namespace endian
struct big_endian_tag {};
struct little_endian_tag {};
#ifdef BOOST_BIG_ENDIAN
#if defined(BOOST_ENDIAN_BIG_BYTE_AVAILABLE)
typedef big_endian_tag native_endian_tag;
#else
#elif defined(BOOST_ENDIAN_LITTLE_BYTE_AVAILABLE)
typedef little_endian_tag native_endian_tag;
#elif defined(BOOST_ENDIAN_BIG_WORD_BYTE_AVAILABLE)
#error Word-swapped big-endian not supported
#elif defined(BOOST_ENDIAN_LITTLE_WORD_BYTE_AVAILABLE)
#error Word-swapped little-endian not supported
#else
#error Unknown endian memory ordering
#endif
// Unrolled loops for loading and storing streams of bytes.

50
include/boost/geometry/extensions/nsphere/geometries/concepts/check.hpp
View File

@@ -1,50 +0,0 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland.
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_GEOMETRY_EXTENSIONS_NSPHERE_GEOMETRIES_CONCEPTS_CHECK_HPP
#define BOOST_GEOMETRY_EXTENSIONS_NSPHERE_GEOMETRIES_CONCEPTS_CHECK_HPP
#include <boost/geometry/geometries/concepts/check.hpp>
namespace boost { namespace geometry
{
#ifndef DOXYGEN_NO_DISPATCH
namespace dispatch
{
template <typename Geometry>
struct check<Geometry, nsphere_tag, true>
: detail::concept_check::check<concepts::ConstNsphere<Geometry> >
{};
template <typename Geometry>
struct check<Geometry, nsphere_tag, false>
: detail::concept_check::check<concepts::Nsphere<Geometry> >
{};
} // namespace dispatch
#endif
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_EXTENSIONS_NSPHERE_GEOMETRIES_CONCEPTS_CHECK_HPP

16
include/boost/geometry/extensions/nsphere/geometries/concepts/nsphere_concept.hpp
View File

@@ -4,6 +4,10 @@
// Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
// This file was modified by Oracle on 2021.
// Modifications copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -117,6 +121,18 @@ public :
};
template <typename Geometry>
struct concept_type<Geometry, nsphere_tag>
{
using type = Nsphere<Geometry>;
};
template <typename Geometry>
struct concept_type<Geometry const, nsphere_tag>
{
using type = ConstNsphere<Geometry>;
};
}}} // namespace boost::geometry::concepts

5
include/boost/geometry/extensions/nsphere/nsphere.hpp
View File

@@ -5,6 +5,10 @@
// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
// Copyright (c) 2011-2015 Adam Wulkiewicz, Lodz, Poland.
// This file was modified by Oracle on 2021.
// Modifications copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
@@ -22,7 +26,6 @@
#include <boost/geometry/extensions/nsphere/core/tags.hpp>
#include <boost/geometry/extensions/nsphere/core/topological_dimension.hpp>
#include <boost/geometry/extensions/nsphere/geometries/concepts/check.hpp>
#include <boost/geometry/extensions/nsphere/geometries/concepts/nsphere_concept.hpp>
#include <boost/geometry/extensions/nsphere/geometries/nsphere.hpp>

107
include/boost/geometry/index/detail/maxmin_heap.hpp Normal file
View File

@@ -0,0 +1,107 @@
// Boost.Geometry
// Copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Licensed under the Boost Software License version 1.0.
// http://www.boost.org/users/license.html
#ifndef BOOST_GEOMETRY_INDEX_DETAIL_MAXMIN_HEAP_HPP
#define BOOST_GEOMETRY_INDEX_DETAIL_MAXMIN_HEAP_HPP
#include <boost/geometry/index/detail/minmax_heap.hpp>
namespace boost { namespace geometry { namespace index { namespace detail
{
template <typename It, typename Compare>
inline void push_maxmin_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
minmax_heap_detail::push_heap<max_call, min_call>(first, last, comp);
}
template <typename It>
inline void push_maxmin_heap(It first, It last)
{
using namespace minmax_heap_detail;
minmax_heap_detail::push_heap<max_call, min_call>(first, last, std::less<>());
}
template <typename It, typename Compare>
inline void pop_top_maxmin_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
pop_heap<max_call, min_call>(first, first, last, comp);
}
template <typename It>
inline void pop_top_maxmin_heap(It first, It last)
{
using namespace minmax_heap_detail;
pop_heap<max_call, min_call>(first, first, last, std::less<>());
}
template <typename It, typename Compare>
inline void pop_bottom_maxmin_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
It bottom = minmax_heap_detail::bottom_heap<max_call>(first, last, comp);
pop_heap<max_call, min_call>(first, bottom, last, comp);
}
template <typename It>
inline void pop_bottom_maxmin_heap(It first, It last)
{
using namespace minmax_heap_detail;
auto&& comp = std::less<>();
It bottom = minmax_heap_detail::bottom_heap<max_call>(first, last, comp);
pop_heap<max_call, min_call>(first, bottom, last, comp);
}
template <typename It, typename Compare>
inline void make_maxmin_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::make_heap<max_call, min_call>(first, last, comp);
}
template <typename It>
inline void make_maxmin_heap(It first, It last)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::make_heap<max_call, min_call>(first, last, std::less<>());
}
template <typename It, typename Compare>
inline bool is_maxmin_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::is_heap<max_call>(first, last, comp);
}
template <typename It>
inline bool is_maxmin_heap(It first, It last)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::is_heap<max_call>(first, last, std::less<>());
}
template <typename It, typename Compare>
inline decltype(auto) bottom_maxmin_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
return *minmax_heap_detail::bottom_heap<max_call>(first, last, comp);
}
template <typename It>
inline decltype(auto) bottom_maxmin_heap(It first, It last)
{
using namespace minmax_heap_detail;
return *minmax_heap_detail::bottom_heap<max_call>(first, last, std::less<>());
}
}}}} // namespace boost::geometry::index::detail
#endif // BOOST_GEOMETRY_INDEX_DETAIL_MAXMIN_HEAP_HPP

519
include/boost/geometry/index/detail/minmax_heap.hpp Normal file
View File

@@ -0,0 +1,519 @@
// Boost.Geometry
// Copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Licensed under the Boost Software License version 1.0.
// http://www.boost.org/users/license.html
#ifndef BOOST_GEOMETRY_INDEX_DETAIL_MINMAX_HEAP_HPP
#define BOOST_GEOMETRY_INDEX_DETAIL_MINMAX_HEAP_HPP
#include <iterator>
#include <type_traits>
#include <utility>
#ifdef _MSC_VER // msvc and clang-win
#include <intrin.h>
#endif
namespace boost { namespace geometry { namespace index { namespace detail
{
// Resources:
// https://en.wikipedia.org/wiki/Min-max_heap
// http://akira.ruc.dk/~keld/teaching/algoritmedesign_f03/Artikler/02/Atkinson86.pdf
// https://probablydance.com/2020/08/31/on-modern-hardware-the-min-max-heap-beats-a-binary-heap/
// https://stackoverflow.com/questions/6531543/efficient-implementation-of-binary-heaps
// https://stackoverflow.com/questions/994593/how-to-do-an-integer-log2-in-c
namespace minmax_heap_detail
{
template <typename T>
using bitsize = std::integral_constant<std::size_t, sizeof(T) * CHAR_BIT>;
template <typename It>
using diff_t = typename std::iterator_traits<It>::difference_type;
template <typename It>
using val_t = typename std::iterator_traits<It>::value_type;
// TODO: In C++20 use std::bit_width()
template <typename T, std::enable_if_t<!std::is_integral<T>::value || (bitsize<T>::value != 32 && bitsize<T>::value != 64), int> = 0>
inline int level(T i)
{
++i;
int r = 0;
while (i >>= 1) { ++r; }
return r;
}
//template <typename T>
//inline int level(T i)
//{
// using std::log2;
// return int(log2(i + 1));
//}
#ifdef _MSC_VER // msvc and clang-win
template <typename T, std::enable_if_t<std::is_integral<T>::value && bitsize<T>::value == 32, int> = 0>
inline int level(T i)
{
unsigned long r = 0;
_BitScanReverse(&r, (unsigned long)(i + 1));
return int(r);
}
template <typename T, std::enable_if_t<std::is_integral<T>::value && bitsize<T>::value == 64, int> = 0>
inline int level(T i)
{
unsigned long r = 0;
#ifdef _WIN64
_BitScanReverse64(&r, (unsigned long long)(i + 1));
#else
if (_BitScanReverse(&r, (unsigned long)((i + 1) >> 32))) { r += 32; }
else { _BitScanReverse(&r, (unsigned long)(i + 1)); }
#endif
return int(r);
}
#elif defined(__clang__) || defined(__GNUC__)
// Only available in gcc-10 and clang-10
//#elif defined(__has_builtin) && __has_builtin(__builtin_clzl) && __has_builtin(__builtin_clzll)
template <typename T, std::enable_if_t<std::is_integral<T>::value && bitsize<T>::value == 32, int> = 0>
inline int level(T i)
{
return 31 - __builtin_clzl((unsigned long)(i + 1));
}
template <typename T, std::enable_if_t<std::is_integral<T>::value && bitsize<T>::value == 64, int> = 0>
inline int level(T i)
{
return 63 - __builtin_clzll((unsigned long long)(i + 1));
}
#else
template <typename T, std::enable_if_t<std::is_integral<T>::value && bitsize<T>::value == 32, int> = 0>
inline int level(T i)
{
++i;
int r = 0;
if (i >= 65536) { r += 16; i >>= 16; }
if (i >= 256) { r += 8; i >>= 8; }
if (i >= 16) { r += 4; i >>= 4; }
if (i >= 4) { r += 2; i >>= 2; }
if (i >= 2) { r += 1; i >>= 1; }
return r;
}
template <typename T, std::enable_if_t<std::is_integral<T>::value && bitsize<T>::value == 64, int> = 0>
inline int level(T i)
{
++i;
int r = 0;
if (i >= 4294967296ll) { r += 32; i >>= 32; }
if (i >= 65536ll) { r += 16; i >>= 16; }
if (i >= 256ll) { r += 8; i >>= 8; }
if (i >= 16ll) { r += 4; i >>= 4; }
if (i >= 4ll) { r += 2; i >>= 2; }
if (i >= 2ll) { r += 1; i >>= 1; }
return r;
}
#endif
// min/max functions only differ in the order of arguments in comp
struct min_call
{
template <typename Compare, typename T1, typename T2>
bool operator()(Compare&& comp, T1 const& v1, T2 const& v2) const
{
return comp(v1, v2);
}
};
struct max_call
{
template <typename Compare, typename T1, typename T2>
bool operator()(Compare&& comp, T1 const& v1, T2 const& v2) const
{
return comp(v2, v1);
}
};
template <typename Call, typename It, typename Compare>
inline void push_heap2(It first, diff_t<It> c, val_t<It> val, Compare comp)
{
while (c > 2)
{
diff_t<It> const g = (c - 3) >> 2; // grandparent index
if (! Call()(comp, val, *(first + g)))
{
break;
}
*(first + c) = std::move(*(first + g));
c = g;
}
*(first + c) = std::move(val);
}
template <typename MinCall, typename MaxCall, typename It, typename Compare>
inline void push_heap1(It first, diff_t<It> c, val_t<It> val, Compare comp)
{
diff_t<It> const p = (c - 1) >> 1; // parent index
if (MinCall()(comp, *(first + p), val))
{
*(first + c) = std::move(*(first + p));
return push_heap2<MaxCall>(first, p, std::move(val), comp);
}
else
{
return push_heap2<MinCall>(first, c, std::move(val), comp);
}
}
template <typename MinCall, typename MaxCall, typename It, typename Compare>
inline void push_heap(It first, It last, Compare comp)
{
diff_t<It> const size = last - first;
if (size < 2)
{
return;
}
diff_t<It> c = size - 1; // back index
val_t<It> val = std::move(*(first + c));
if (level(c) % 2 == 0) // is min level
{
push_heap1<MinCall, MaxCall>(first, c, std::move(val), comp);
}
else
{
push_heap1<MaxCall, MinCall>(first, c, std::move(val), comp);
}
}
template <typename Call, typename It, typename Compare>
inline diff_t<It> pick_grandchild4(It first, diff_t<It> f, Compare comp)
{
It it = first + f;
diff_t<It> m1 = Call()(comp, *(it), *(it + 1)) ? f : f + 1;
diff_t<It> m2 = Call()(comp, *(it + 2), *(it + 3)) ? f + 2 : f + 3;
return Call()(comp, *(first + m1), *(first + m2)) ? m1 : m2;
}
//template <typename Call, typename It, typename Compare>
//inline diff_t<It> pick_descendant(It first, diff_t<It> f, diff_t<It> l, Compare comp)
//{
// diff_t<It> m = f;
// for (++f; f != l; ++f)
// {
// if (Call()(comp, *(first + f), *(first + m)))
// {
// m = f;
// }
// }
// return m;
//}
template <typename Call, typename It, typename Compare>
inline void pop_heap1(It first, diff_t<It> p, diff_t<It> size, val_t<It> val, Compare comp)
{
if (size >= 7) // grandparent of 4 grandchildren is possible
{
diff_t<It> const last_g = (size - 3) >> 2; // grandparent of the element behind back
while (p < last_g) // p is grandparent of 4 grandchildren
{
diff_t<It> const ll = 4 * p + 3;
diff_t<It> const m = pick_grandchild4<Call>(first, ll, comp);
if (! Call()(comp, *(first + m), val))
{
break;
}
*(first + p) = std::move(*(first + m));
diff_t<It> par = (m - 1) >> 1;
if (Call()(comp, *(first + par), val))
{
using std::swap;
swap(*(first + par), val);
}
p = m;
}
}
if (size >= 2 && p <= ((size - 2) >> 1)) // at least one child
{
diff_t<It> const l = 2 * p + 1;
diff_t<It> m = l; // left child
if (size >= 3 && p <= ((size - 3) >> 1)) // at least two children
{
// m = left child
diff_t<It> m2 = l + 1; // right child
if (size >= 4 && p <= ((size - 4) >> 2)) // at least two children and one grandchild
{
diff_t<It> const ll = 2 * l + 1;
m = ll; // left most grandchild
// m2 = right child
if (size >= 5 && p <= ((size - 5) >> 2)) // at least two children and two grandchildren
{
m = Call()(comp, *(first + ll), *(first + ll + 1)) ? ll : (ll + 1); // greater of the left grandchildren
// m2 = right child
if (size >= 6 && p <= ((size - 6) >> 2)) // at least two children and three grandchildren
{
// m = greater of the left grandchildren
m2 = ll + 2; // third grandchild
}
}
}
m = Call()(comp, *(first + m), *(first + m2)) ? m : m2;
}
if (Call()(comp, *(first + m), val))
{
*(first + p) = std::move(*(first + m));
if (m >= 3 && p <= ((m - 3) >> 2)) // is p grandparent of m
{
diff_t<It> par = (m - 1) >> 1;
if (Call()(comp, *(first + par), val))
{
using std::swap;
swap(*(first + par), val);
}
}
p = m;
}
}
*(first + p) = std::move(val);
}
template <typename MinCall, typename MaxCall, typename It, typename Compare>
inline void pop_heap(It first, It el, It last, Compare comp)
{
diff_t<It> size = last - first;
if (size < 2)
{
return;
}
--last;
val_t<It> val = std::move(*last);
*last = std::move(*el);
// Ignore the last element
--size;
diff_t<It> p = el - first;
if (level(p) % 2 == 0) // is min level
{
pop_heap1<MinCall>(first, p, size, std::move(val), comp);
}
else
{
pop_heap1<MaxCall>(first, p, size, std::move(val), comp);
}
}
template <typename MinCall, typename MaxCall, typename It, typename Compare>
inline void make_heap(It first, It last, Compare comp)
{
diff_t<It> size = last - first;
diff_t<It> p = size / 2;
if (p <= 0)
{
return;
}
int level_p = level(p - 1);
diff_t<It> level_f = (diff_t<It>(1) << level_p) - 1;
while (p > 0)
{
--p;
val_t<It> val = std::move(*(first + p));
if (level_p % 2 == 0) // is min level
{
pop_heap1<MinCall>(first, p, size, std::move(val), comp);
}
else
{
pop_heap1<MaxCall>(first, p, size, std::move(val), comp);
}
if (p == level_f)
{
--level_p;
level_f >>= 1;
}
}
}
template <typename Call, typename It, typename Compare>
inline bool is_heap(It first, It last, Compare comp)
{
diff_t<It> const size = last - first;
diff_t<It> pow2 = 4;
bool is_min_level = false;
for (diff_t<It> i = 1; i < size; ++i)
{
if (i == pow2 - 1)
{
pow2 *= 2;
is_min_level = ! is_min_level;
}
diff_t<It> const p = (i - 1) >> 1;
if (is_min_level)
{
if (Call()(comp, *(first + p), *(first + i)))
{
return false;
}
}
else
{
if (Call()(comp, *(first + i), *(first + p)))
{
return false;
}
}
if (i >= 3)
{
diff_t<It> const g = (p - 1) >> 1;
if (is_min_level)
{
if (Call()(comp, *(first + i), *(first + g)))
{
return false;
}
}
else
{
if (Call()(comp, *(first + g), *(first + i)))
{
return false;
}
}
}
}
return true;
}
template <typename Call, typename It, typename Compare>
inline It bottom_heap(It first, It last, Compare comp)
{
diff_t<It> const size = last - first;
return size <= 1 ? first :
size == 2 ? (first + 1) :
Call()(comp, *(first + 1), *(first + 2)) ? (first + 2) : (first + 1);
}
} // namespace minmax_heap_detail
template <typename It, typename Compare>
inline void push_minmax_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
minmax_heap_detail::push_heap<min_call, max_call>(first, last, comp);
}
template <typename It>
inline void push_minmax_heap(It first, It last)
{
using namespace minmax_heap_detail;
minmax_heap_detail::push_heap<min_call, max_call>(first, last, std::less<>());
}
template <typename It, typename Compare>
inline void pop_top_minmax_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
pop_heap<min_call, max_call>(first, first, last, comp);
}
template <typename It>
inline void pop_top_minmax_heap(It first, It last)
{
using namespace minmax_heap_detail;
pop_heap<min_call, max_call>(first, first, last, std::less<>());
}
template <typename It, typename Compare>
inline void pop_bottom_minmax_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
It bottom = minmax_heap_detail::bottom_heap<min_call>(first, last, comp);
pop_heap<min_call, max_call>(first, bottom, last, comp);
}
template <typename It>
inline void pop_bottom_minmax_heap(It first, It last)
{
using namespace minmax_heap_detail;
auto&& comp = std::less<>();
It bottom = minmax_heap_detail::bottom_heap<min_call>(first, last, comp);
pop_heap<min_call, max_call>(first, bottom, last, comp);
}
template <typename It, typename Compare>
inline void make_minmax_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::make_heap<min_call, max_call>(first, last, comp);
}
template <typename It>
inline void make_minmax_heap(It first, It last)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::make_heap<min_call, max_call>(first, last, std::less<>());
}
template <typename It, typename Compare>
inline bool is_minmax_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::is_heap<min_call>(first, last, comp);
}
template <typename It>
inline bool is_minmax_heap(It first, It last)
{
using namespace minmax_heap_detail;
return minmax_heap_detail::is_heap<min_call>(first, last, std::less<>());
}
template <typename It, typename Compare>
inline decltype(auto) bottom_minmax_heap(It first, It last, Compare comp)
{
using namespace minmax_heap_detail;
return *minmax_heap_detail::bottom_heap<min_call>(first, last, comp);
}
template <typename It>
inline decltype(auto) bottom_minmax_heap(It first, It last)
{
using namespace minmax_heap_detail;
return *minmax_heap_detail::bottom_heap<min_call>(first, last, std::less<>());
}
}}}} // namespace boost::geometry::index::detail
#endif // BOOST_GEOMETRY_INDEX_DETAIL_MINMAX_HEAP_HPP

8
include/boost/geometry/index/detail/predicates.hpp
View File

@@ -622,11 +622,13 @@ struct predicates_check_impl<std::tuple<Ts...>, Tag, First, Last>
}
};
template <typename Tag, std::size_t First, std::size_t Last, typename Predicates, typename Value, typename Indexable, typename Strategy>
template <typename Tag, typename Predicates, typename Value, typename Indexable, typename Strategy>
inline bool predicates_check(Predicates const& p, Value const& v, Indexable const& i, Strategy const& s)
{
return detail::predicates_check_impl<Predicates, Tag, First, Last>
::apply(p, v, i, s);
return detail::predicates_check_impl
<
Predicates, Tag, 0, predicates_length<Predicates>::value
>::apply(p, v, i, s);
}
// ------------------------------------------------------------------ //

150
include/boost/geometry/index/detail/priority_dequeue.hpp Normal file
View File

@@ -0,0 +1,150 @@
// Boost.Geometry
// Copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Licensed under the Boost Software License version 1.0.
// http://www.boost.org/users/license.html
#ifndef BOOST_GEOMETRY_INDEX_DETAIL_PRIORITY_DEQUEUE_HPP
#define BOOST_GEOMETRY_INDEX_DETAIL_PRIORITY_DEQUEUE_HPP
#include <vector>
#include <boost/geometry/index/detail/maxmin_heap.hpp>
namespace boost { namespace geometry { namespace index { namespace detail
{
template
<
typename T,
typename Container = std::vector<T>,
typename Compare = std::less<typename Container::value_type>
>
class priority_dequeue
{
public:
using container_type = Container;
using value_compare = Compare;
using value_type = typename Container::value_type;
using size_type = typename Container::size_type;
using reference = typename Container::reference;
using const_reference = typename Container::const_reference;
priority_dequeue()
: c(), comp()
{}
explicit priority_dequeue(Compare const& compare)
: c(), comp(compare)
{}
priority_dequeue(Compare const& compare, Container const& cont)
: c(cont), comp(compare)
{
make_maxmin_heap(c.begin(), c.end(), comp);
}
priority_dequeue(Compare const& compare, Container&& cont)
: c(std::move(cont)), comp(compare)
{
make_maxmin_heap(c.begin(), c.end(), comp);
}
template <typename It>
priority_dequeue(It first, It last)
: c(first, last), comp()
{
make_maxmin_heap(c.begin(), c.end(), comp);
}
template <typename It>
priority_dequeue(It first, It last, Compare const& compare)
: c(first, last), comp(compare)
{
make_maxmin_heap(c.begin(), c.end(), comp);
}
template <typename It>
priority_dequeue(It first, It last, Compare const& compare, Container const& cont)
: c(cont), comp(compare)
{
c.insert(first, last);
make_maxmin_heap(c.begin(), c.end(), comp);
}
template <typename It>
priority_dequeue(It first, It last, Compare const& compare, Container&& cont)
: c(std::move(cont)), comp(compare)
{
c.insert(first, last);
make_maxmin_heap(c.begin(), c.end(), comp);
}
const_reference top() const
{
return *c.begin();
}
const_reference bottom() const
{
return bottom_maxmin_heap(c.begin(), c.end(), comp);
}
void push(const value_type& value)
{
c.push_back(value);
push_maxmin_heap(c.begin(), c.end(), comp);
}
void push(value_type&& value)
{
c.push_back(std::move(value));
push_maxmin_heap(c.begin(), c.end(), comp);
}
template <typename ...Args>
void emplace(Args&& ...args)
{
c.emplace_back(std::forward<Args>(args)...);
push_maxmin_heap(c.begin(), c.end(), comp);
}
void pop_top()
{
pop_top_maxmin_heap(c.begin(), c.end(), comp);
c.pop_back();
}
void pop_bottom()
{
pop_bottom_maxmin_heap(c.begin(), c.end(), comp);
c.pop_back();
}
bool empty() const
{
return c.empty();
}
size_t size() const
{
return c.size();
}
void swap(priority_dequeue& other)
{
using std::swap;
std::swap(c, other.c);
std::swap(comp, other.comp);
}
protected:
Container c;
Compare comp;
};
}}}} // namespace boost::geometry::index::detail
#endif // BOOST_GEOMETRY_INDEX_DETAIL_PRIORITY_DEQUEUE_HPP

107
include/boost/geometry/index/detail/rtree/query_iterators.hpp
View File

@@ -21,8 +21,6 @@
#include <boost/geometry/index/detail/rtree/visitors/distance_query.hpp>
#include <boost/geometry/index/detail/rtree/visitors/spatial_query.hpp>
//#define BOOST_GEOMETRY_INDEX_DETAIL_QUERY_ITERATORS_USE_MOVE
namespace boost { namespace geometry { namespace index { namespace detail { namespace rtree { namespace iterators {
template <typename Value, typename Allocators>
@@ -69,13 +67,8 @@ struct end_query_iterator
template <typename MembersHolder, typename Predicates>
class spatial_query_iterator
{
typedef typename MembersHolder::parameters_type parameters_type;
typedef typename MembersHolder::translator_type translator_type;
typedef typename MembersHolder::allocators_type allocators_type;
typedef visitors::spatial_query_incremental<MembersHolder, Predicates> visitor_type;
typedef typename visitor_type::node_pointer node_pointer;
public:
typedef std::forward_iterator_tag iterator_category;
typedef typename MembersHolder::value_type value_type;
@@ -83,32 +76,31 @@ public:
typedef typename allocators_type::difference_type difference_type;
typedef typename allocators_type::const_pointer pointer;
inline spatial_query_iterator()
spatial_query_iterator() = default;
explicit spatial_query_iterator(Predicates const& pred)
: m_impl(pred)
{}
inline spatial_query_iterator(parameters_type const& par, translator_type const& t, Predicates const& p)
: m_visitor(par, t, p)
{}
inline spatial_query_iterator(node_pointer root, parameters_type const& par, translator_type const& t, Predicates const& p)
: m_visitor(par, t, p)
spatial_query_iterator(MembersHolder const& members, Predicates const& pred)
: m_impl(members, pred)
{
m_visitor.initialize(root);
m_impl.initialize(members);
}
reference operator*() const
{
return m_visitor.dereference();
return m_impl.dereference();
}
const value_type * operator->() const
{
return boost::addressof(m_visitor.dereference());
return boost::addressof(m_impl.dereference());
}
spatial_query_iterator & operator++()
{
m_visitor.increment();
m_impl.increment();
return *this;
}
@@ -121,33 +113,28 @@ public:
friend bool operator==(spatial_query_iterator const& l, spatial_query_iterator const& r)
{
return l.m_visitor == r.m_visitor;
return l.m_impl == r.m_impl;
}
friend bool operator==(spatial_query_iterator const& l, end_query_iterator<value_type, allocators_type> const& /*r*/)
{
return l.m_visitor.is_end();
return l.m_impl.is_end();
}
friend bool operator==(end_query_iterator<value_type, allocators_type> const& /*l*/, spatial_query_iterator const& r)
{
return r.m_visitor.is_end();
return r.m_impl.is_end();
}
private:
visitor_type m_visitor;
visitors::spatial_query_incremental<MembersHolder, Predicates> m_impl;
};
template <typename MembersHolder, typename Predicates, unsigned NearestPredicateIndex>
template <typename MembersHolder, typename Predicates>
class distance_query_iterator
{
typedef typename MembersHolder::parameters_type parameters_type;
typedef typename MembersHolder::translator_type translator_type;
typedef typename MembersHolder::allocators_type allocators_type;
typedef visitors::distance_query_incremental<MembersHolder, Predicates, NearestPredicateIndex> visitor_type;
typedef typename visitor_type::node_pointer node_pointer;
public:
typedef std::forward_iterator_tag iterator_category;
typedef typename MembersHolder::value_type value_type;
@@ -155,32 +142,31 @@ public:
typedef typename allocators_type::difference_type difference_type;
typedef typename allocators_type::const_pointer pointer;
inline distance_query_iterator()
distance_query_iterator() = default;
explicit distance_query_iterator(Predicates const& pred)
: m_impl(pred)
{}
inline distance_query_iterator(parameters_type const& par, translator_type const& t, Predicates const& p)
: m_visitor(par, t, p)
{}
inline distance_query_iterator(node_pointer root, parameters_type const& par, translator_type const& t, Predicates const& p)
: m_visitor(par, t, p)
distance_query_iterator(MembersHolder const& members, Predicates const& pred)
: m_impl(members, pred)
{
m_visitor.initialize(root);
m_impl.initialize(members);
}
reference operator*() const
{
return m_visitor.dereference();
return m_impl.dereference();
}
const value_type * operator->() const
{
return boost::addressof(m_visitor.dereference());
return boost::addressof(m_impl.dereference());
}
distance_query_iterator & operator++()
{
m_visitor.increment();
m_impl.increment();
return *this;
}
@@ -193,21 +179,21 @@ public:
friend bool operator==(distance_query_iterator const& l, distance_query_iterator const& r)
{
return l.m_visitor == r.m_visitor;
return l.m_impl == r.m_impl;
}
friend bool operator==(distance_query_iterator const& l, end_query_iterator<value_type, allocators_type> const& /*r*/)
{
return l.m_visitor.is_end();
return l.m_impl.is_end();
}
friend bool operator==(end_query_iterator<value_type, allocators_type> const& /*l*/, distance_query_iterator const& r)
{
return r.m_visitor.is_end();
return r.m_impl.is_end();
}
private:
visitor_type m_visitor;
visitors::distance_query_incremental<MembersHolder, Predicates> m_impl;
};
@@ -284,8 +270,7 @@ public:
typedef typename Allocators::difference_type difference_type;
typedef typename Allocators::const_pointer pointer;
query_iterator()
{}
query_iterator() = default;
template <typename It>
query_iterator(It const& it)
@@ -300,7 +285,6 @@ public:
: m_ptr(o.m_ptr.get() ? o.m_ptr->clone() : 0)
{}
#ifndef BOOST_GEOMETRY_INDEX_DETAIL_QUERY_ITERATORS_USE_MOVE
query_iterator & operator=(query_iterator const& o)
{
if ( this != boost::addressof(o) )
@@ -309,12 +293,13 @@ public:
}
return *this;
}
#ifndef BOOST_NO_CXX11_RVALUE_REFERENCES
query_iterator(query_iterator && o)
: m_ptr(0)
{
m_ptr.swap(o.m_ptr);
}
query_iterator & operator=(query_iterator && o)
{
if ( this != boost::addressof(o) )
@@ -324,34 +309,6 @@ public:
}
return *this;
}
#endif
#else // !BOOST_GEOMETRY_INDEX_DETAIL_QUERY_ITERATORS_USE_MOVE
private:
BOOST_COPYABLE_AND_MOVABLE(query_iterator)
public:
query_iterator & operator=(BOOST_COPY_ASSIGN_REF(query_iterator) o)
{
if ( this != boost::addressof(o) )
{
m_ptr.reset(o.m_ptr.get() ? o.m_ptr->clone() : 0);
}
return *this;
}
query_iterator(BOOST_RV_REF(query_iterator) o)
: m_ptr(0)
{
m_ptr.swap(o.m_ptr);
}
query_iterator & operator=(BOOST_RV_REF(query_iterator) o)
{
if ( this != boost::addressof(o) )
{
m_ptr.swap(o.m_ptr);
o.m_ptr.reset();
}
return *this;
}
#endif // BOOST_GEOMETRY_INDEX_DETAIL_QUERY_ITERATORS_USE_MOVE
reference operator*() const
{

667
include/boost/geometry/index/detail/rtree/visitors/distance_query.hpp
View File

@@ -15,8 +15,11 @@
#ifndef BOOST_GEOMETRY_INDEX_DETAIL_RTREE_VISITORS_DISTANCE_QUERY_HPP
#define BOOST_GEOMETRY_INDEX_DETAIL_RTREE_VISITORS_DISTANCE_QUERY_HPP
#include <queue>
#include <boost/geometry/index/detail/distance_predicates.hpp>
#include <boost/geometry/index/detail/predicates.hpp>
#include <boost/geometry/index/detail/priority_dequeue.hpp>
#include <boost/geometry/index/detail/rtree/node/node_elements.hpp>
#include <boost/geometry/index/detail/translator.hpp>
#include <boost/geometry/index/parameters.hpp>
@@ -46,90 +49,114 @@ struct pair_first_greater
}
};
template <typename T, typename Comp>
struct priority_dequeue : index::detail::priority_dequeue<T, std::vector<T>, Comp>
{
priority_dequeue() = default;
//void reserve(typename std::vector<T>::size_type n)
//{
// this->c.reserve(n);
//}
//void clear()
//{
// this->c.clear();
//}
};
template <typename Value, typename Translator, typename DistanceType, typename OutIt>
template <typename T, typename Comp>
struct priority_queue : std::priority_queue<T, std::vector<T>, Comp>
{
priority_queue() = default;
//void reserve(typename std::vector<T>::size_type n)
//{
// this->c.reserve(n);
//}
void clear()
{
this->c.clear();
}
};
struct branch_data_comp
{
template <typename BranchData>
bool operator()(BranchData const& b1, BranchData const& b2) const
{
return b1.distance > b2.distance || (b1.distance == b2.distance && b1.reverse_level > b2.reverse_level);
}
};
template <typename DistanceType, typename Value>
class distance_query_result
{
using neighbor_data = std::pair<DistanceType, const Value *>;
using neighbors_type = std::vector<neighbor_data>;
using size_type = typename neighbors_type::size_type;
public:
typedef DistanceType distance_type;
inline explicit distance_query_result(size_t k, OutIt out_it)
: m_count(k), m_out_it(out_it)
inline distance_query_result(size_type k)
: m_count(k)
{
BOOST_GEOMETRY_INDEX_ASSERT(0 < m_count, "Number of neighbors should be greater than 0");
m_neighbors.reserve(m_count);
}
inline void store(Value const& val, distance_type const& curr_comp_dist)
// NOTE: Do not call if max_count() == 0
inline void store(DistanceType const& distance, const Value * value_ptr)
{
if ( m_neighbors.size() < m_count )
if (m_neighbors.size() < m_count)
{
m_neighbors.push_back(std::make_pair(curr_comp_dist, val));
m_neighbors.push_back(std::make_pair(distance, value_ptr));
if ( m_neighbors.size() == m_count )
std::make_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
}
else
{
if ( curr_comp_dist < m_neighbors.front().first )
if (m_neighbors.size() == m_count)
{
std::pop_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
m_neighbors.back().first = curr_comp_dist;
m_neighbors.back().second = val;
std::push_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
std::make_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
}
}
else if (distance < m_neighbors.front().first)
{
std::pop_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
m_neighbors.back().first = distance;
m_neighbors.back().second = value_ptr;
std::push_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
}
}
inline bool has_enough_neighbors() const
// NOTE: Do not call if max_count() == 0
inline bool ignore_branch(DistanceType const& distance) const
{
return m_count <= m_neighbors.size();
return m_neighbors.size() == m_count
&& m_neighbors.front().first <= distance;
}
inline distance_type greatest_comparable_distance() const
template <typename OutIt>
inline size_type finish(OutIt out_it) const
{
// greatest distance is in the first neighbor only
// if there is at least m_count values found
// this is just for safety reasons since is_comparable_distance_valid() is checked earlier
// TODO - may be replaced by ASSERT
return m_neighbors.size() < m_count
? (std::numeric_limits<distance_type>::max)()
: m_neighbors.front().first;
}
inline size_t finish()
{
typedef typename std::vector< std::pair<distance_type, Value> >::const_iterator neighbors_iterator;
for ( neighbors_iterator it = m_neighbors.begin() ; it != m_neighbors.end() ; ++it, ++m_out_it )
*m_out_it = it->second;
for (auto const& p : m_neighbors)
{
*out_it = *(p.second);
++out_it;
}
return m_neighbors.size();
}
private:
size_t m_count;
OutIt m_out_it;
size_type max_count() const
{
return m_count;
}
std::vector< std::pair<distance_type, Value> > m_neighbors;
private:
size_type m_count;
neighbors_type m_neighbors;
};
template
<
typename MembersHolder,
typename Predicates,
std::size_t DistancePredicateIndex,
typename OutIter
>
template <typename MembersHolder, typename Predicates>
class distance_query
: public MembersHolder::visitor_const
{
public:
typedef typename MembersHolder::value_type value_type;
typedef typename MembersHolder::box_type box_type;
typedef typename MembersHolder::parameters_type parameters_type;
typedef typename MembersHolder::translator_type translator_type;
typedef typename MembersHolder::allocators_type allocators_type;
typedef typename index::detail::strategy_type<parameters_type>::type strategy_type;
@@ -137,7 +164,10 @@ public:
typedef typename MembersHolder::internal_node internal_node;
typedef typename MembersHolder::leaf leaf;
typedef index::detail::predicates_element<DistancePredicateIndex, Predicates> nearest_predicate_access;
typedef index::detail::predicates_element
<
index::detail::predicates_find_distance<Predicates>::value, Predicates
> nearest_predicate_access;
typedef typename nearest_predicate_access::type nearest_predicate_type;
typedef typename indexable_type<translator_type>::type indexable_type;
@@ -146,147 +176,141 @@ public:
typedef typename calculate_value_distance::result_type value_distance_type;
typedef typename calculate_node_distance::result_type node_distance_type;
static const std::size_t predicates_len = index::detail::predicates_length<Predicates>::value;
typedef typename MembersHolder::size_type size_type;
typedef typename MembersHolder::node_pointer node_pointer;
inline distance_query(parameters_type const& parameters, translator_type const& translator, Predicates const& pred, OutIter out_it)
: m_parameters(parameters), m_translator(translator)
using neighbor_data = std::pair<value_distance_type, const value_type *>;
using neighbors_type = std::vector<neighbor_data>;
struct branch_data
{
branch_data(node_distance_type d, size_type rl, node_pointer p)
: distance(d), reverse_level(rl), ptr(p)
{}
node_distance_type distance;
size_type reverse_level;
node_pointer ptr;
};
using branches_type = priority_queue<branch_data, branch_data_comp>;
public:
distance_query(MembersHolder const& members, Predicates const& pred)
: m_tr(members.translator())
, m_strategy(index::detail::get_strategy(members.parameters()))
, m_pred(pred)
, m_result(nearest_predicate_access::get(m_pred).count, out_it)
, m_strategy(index::detail::get_strategy(parameters))
{}
inline void operator()(internal_node const& n)
{
typedef typename rtree::elements_type<internal_node>::type elements_type;
// array of active nodes
typedef typename index::detail::rtree::container_from_elements_type<
elements_type,
std::pair<node_distance_type, typename allocators_type::node_pointer>
>::type active_branch_list_type;
active_branch_list_type active_branch_list;
active_branch_list.reserve(m_parameters.get_max_elements());
elements_type const& elements = rtree::elements(n);
// fill array of nodes meeting predicates
for (typename elements_type::const_iterator it = elements.begin();
it != elements.end(); ++it)
{
// if current node meets predicates
// 0 - dummy value
if ( index::detail::predicates_check
<
index::detail::bounds_tag, 0, predicates_len
>(m_pred, 0, it->first, m_strategy) )
{
// calculate node's distance(s) for distance predicate
node_distance_type node_distance;
// if distance isn't ok - move to the next node
if ( !calculate_node_distance::apply(predicate(), it->first,
m_strategy, node_distance) )
{
continue;
}
// if current node is further than found neighbors - don't analyze it
if ( m_result.has_enough_neighbors() &&
is_node_prunable(m_result.greatest_comparable_distance(), node_distance) )
{
continue;
}
// add current node's data into the list
active_branch_list.push_back( std::make_pair(node_distance, it->second) );
}
}
// if there aren't any nodes in ABL - return
if ( active_branch_list.empty() )
return;
// sort array
std::sort(active_branch_list.begin(), active_branch_list.end(), pair_first_less());
// recursively visit nodes
for ( typename active_branch_list_type::const_iterator it = active_branch_list.begin();
it != active_branch_list.end() ; ++it )
{
// if current node is further than furthest neighbor, the rest of nodes also will be further
if ( m_result.has_enough_neighbors() &&
is_node_prunable(m_result.greatest_comparable_distance(), it->first) )
break;
rtree::apply_visitor(*this, *(it->second));
}
// ALTERNATIVE VERSION - use heap instead of sorted container
// It seems to be faster for greater MaxElements and slower otherwise
// CONSIDER: using one global container/heap for active branches
// instead of a sorted container per level
// This would also change the way how branches are traversed!
// The same may be applied to the iterative version which btw suffers
// from the copying of the whole containers on resize of the ABLs container
//// make a heap
//std::make_heap(active_branch_list.begin(), active_branch_list.end(), pair_first_greater());
//// recursively visit nodes
//while ( !active_branch_list.empty() )
//{
// //if current node is further than furthest neighbor, the rest of nodes also will be further
// if ( m_result.has_enough_neighbors()
// && is_node_prunable(m_result.greatest_comparable_distance(), active_branch_list.front().first) )
// {
// break;
// }
// rtree::apply_visitor(*this, *(active_branch_list.front().second));
// std::pop_heap(active_branch_list.begin(), active_branch_list.end(), pair_first_greater());
// active_branch_list.pop_back();
//}
m_neighbors.reserve((std::min)(members.values_count, size_type(max_count())));
//m_branches.reserve(members.parameters().get_min_elements() * members.leafs_level); ?
// min, max or average?
}
inline void operator()(leaf const& n)
template <typename OutIter>
size_type apply(MembersHolder const& members, OutIter out_it)
{
typedef typename rtree::elements_type<leaf>::type elements_type;
elements_type const& elements = rtree::elements(n);
// search leaf for closest value meeting predicates
for (typename elements_type::const_iterator it = elements.begin();
it != elements.end(); ++it)
{
// if value meets predicates
if ( index::detail::predicates_check
<
index::detail::value_tag, 0, predicates_len
>(m_pred, *it, m_translator(*it), m_strategy) )
{
// calculate values distance for distance predicate
value_distance_type value_distance;
// if distance is ok
if ( calculate_value_distance::apply(predicate(), m_translator(*it),
m_strategy, value_distance) )
{
// store value
m_result.store(*it, value_distance);
}
}
}
}
inline size_t finish()
{
return m_result.finish();
return apply(members.root, members.leafs_level, out_it);
}
private:
template <typename Distance>
static inline bool is_node_prunable(Distance const& greatest_dist, node_distance_type const& d)
template <typename OutIter>
size_type apply(node_pointer ptr, size_type reverse_level, OutIter out_it)
{
return greatest_dist <= d;
namespace id = index::detail;
if (max_count() <= 0)
{
return 0;
}
for (;;)
{
if (reverse_level > 0)
{
internal_node& n = rtree::get<internal_node>(*ptr);
// fill array of nodes meeting predicates
for (auto const& p : rtree::elements(n))
{
node_distance_type node_distance; // for distance predicate
// if current node meets predicates (0 is dummy value)
if (id::predicates_check<id::bounds_tag>(m_pred, 0, p.first, m_strategy)
// and if distance is ok
&& calculate_node_distance::apply(predicate(), p.first, m_strategy, node_distance)
// and if current node is closer than the furthest neighbor
&& ! ignore_branch(node_distance))
{
// add current node's data into the list
m_branches.push(branch_data(node_distance, reverse_level - 1, p.second));
}
}
}
else
{
leaf& n = rtree::get<leaf>(*ptr);
// search leaf for closest value meeting predicates
for (auto const& v : rtree::elements(n))
{
value_distance_type value_distance; // for distance predicate
// if value meets predicates
if (id::predicates_check<id::value_tag>(m_pred, v, m_tr(v), m_strategy)
// and if distance is ok
&& calculate_value_distance::apply(predicate(), m_tr(v), m_strategy, value_distance))
{
// store value
store_value(value_distance, boost::addressof(v));
}
}
}
if (m_branches.empty()
|| ignore_branch(m_branches.top().distance))
{
break;
}
ptr = m_branches.top().ptr;
reverse_level = m_branches.top().reverse_level;
m_branches.pop();
}
for (auto const& p : m_neighbors)
{
*out_it = *(p.second);
++out_it;
}
return m_neighbors.size();
}
bool ignore_branch(node_distance_type const& node_distance) const
{
return m_neighbors.size() == max_count()
&& m_neighbors.front().first <= node_distance;
}
void store_value(value_distance_type value_distance, const value_type * value_ptr)
{
if (m_neighbors.size() < max_count())
{
m_neighbors.push_back(std::make_pair(value_distance, value_ptr));
if (m_neighbors.size() == max_count())
{
std::make_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
}
}
else if (value_distance < m_neighbors.front().first)
{
std::pop_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
m_neighbors.back() = std::make_pair(value_distance, value_ptr);
std::push_heap(m_neighbors.begin(), m_neighbors.end(), pair_first_less());
}
}
std::size_t max_count() const
{
return nearest_predicate_access::get(m_pred).count;
}
nearest_predicate_type const& predicate() const
@@ -294,24 +318,18 @@ private:
return nearest_predicate_access::get(m_pred);
}
parameters_type const& m_parameters;
translator_type const& m_translator;
Predicates m_pred;
distance_query_result<value_type, translator_type, value_distance_type, OutIter> m_result;
translator_type const& m_tr;
strategy_type m_strategy;
Predicates const& m_pred;
branches_type m_branches;
neighbors_type m_neighbors;
};
template <
typename MembersHolder,
typename Predicates,
std::size_t DistancePredicateIndex
>
template <typename MembersHolder, typename Predicates>
class distance_query_incremental
: public MembersHolder::visitor_const
{
public:
typedef typename MembersHolder::value_type value_type;
typedef typename MembersHolder::box_type box_type;
typedef typename MembersHolder::parameters_type parameters_type;
@@ -324,7 +342,10 @@ public:
typedef typename MembersHolder::internal_node internal_node;
typedef typename MembersHolder::leaf leaf;
typedef index::detail::predicates_element<DistancePredicateIndex, Predicates> nearest_predicate_access;
typedef index::detail::predicates_element
<
index::detail::predicates_find_distance<Predicates>::value, Predicates
> nearest_predicate_access;
typedef typename nearest_predicate_access::type nearest_predicate_type;
typedef typename indexable_type<translator_type>::type indexable_type;
@@ -337,90 +358,114 @@ public:
typedef typename allocators_type::const_reference const_reference;
typedef typename allocators_type::node_pointer node_pointer;
static const std::size_t predicates_len = index::detail::predicates_length<Predicates>::value;
typedef typename rtree::elements_type<internal_node>::type internal_elements;
typedef typename internal_elements::const_iterator internal_iterator;
typedef typename rtree::elements_type<leaf>::type leaf_elements;
typedef std::pair<node_distance_type, node_pointer> branch_data;
typedef std::vector<branch_data> internal_heap_type;
using neighbor_data = std::pair<value_distance_type, const value_type *>;
using neighbors_type = priority_dequeue<neighbor_data, pair_first_greater>;
struct branch_data
{
branch_data(node_distance_type d, size_type rl, node_pointer p)
: distance(d), reverse_level(rl), ptr(p)
{}
node_distance_type distance;
size_type reverse_level;
node_pointer ptr;
};
using branches_type = priority_queue<branch_data, branch_data_comp>;
public:
inline distance_query_incremental()
: m_translator(NULL)
: m_tr(nullptr)
// , m_strategy()
// , m_pred()
, current_neighbor((std::numeric_limits<size_type>::max)())
// , m_strategy_type()
, m_neighbors_count(0)
, m_neighbor_ptr(nullptr)
{}
inline distance_query_incremental(parameters_type const& params, translator_type const& translator, Predicates const& pred)
: m_translator(::boost::addressof(translator))
inline distance_query_incremental(Predicates const& pred)
: m_tr(nullptr)
// , m_strategy()
, m_pred(pred)
, current_neighbor((std::numeric_limits<size_type>::max)())
, m_strategy(index::detail::get_strategy(params))
{
BOOST_GEOMETRY_INDEX_ASSERT(0 < max_count(), "k must be greather than 0");
}
, m_neighbors_count(0)
, m_neighbor_ptr(nullptr)
{}
inline distance_query_incremental(MembersHolder const& members, Predicates const& pred)
: m_tr(::boost::addressof(members.translator()))
, m_strategy(index::detail::get_strategy(members.parameters()))
, m_pred(pred)
, m_neighbors_count(0)
, m_neighbor_ptr(nullptr)
{}
const_reference dereference() const
{
return *(neighbors[current_neighbor].second);
return *m_neighbor_ptr;
}
void initialize(node_pointer root)
void initialize(MembersHolder const& members)
{
rtree::apply_visitor(*this, *root);
increment();
if (0 < max_count())
{
apply(members.root, members.leafs_level);
increment();
}
}
void increment()
{
for (;;)
{
size_type new_neighbor = current_neighbor == (std::numeric_limits<size_type>::max)() ? 0 : current_neighbor + 1;
if ( internal_heap.empty() )
if (m_branches.empty())
{
if ( new_neighbor < neighbors.size() )
current_neighbor = new_neighbor;
// there exists a next closest neighbor so we can increment
if (! m_neighbors.empty())
{
m_neighbor_ptr = m_neighbors.top().second;
++m_neighbors_count;
m_neighbors.pop_top();
}
else
{
current_neighbor = (std::numeric_limits<size_type>::max)();
// clear() is used to disable the condition above
neighbors.clear();
// there aren't any neighbors left, end
m_neighbor_ptr = nullptr;
m_neighbors_count = max_count();
}
return;
}
else
{
branch_data const& closest_branch = internal_heap.front();
node_distance_type const& closest_distance = closest_branch.first;
branch_data const& closest_branch = m_branches.top();
// if there are no nodes which can have closer values, set new value
if ( new_neighbor < neighbors.size() &&
// NOTE: In order to use <= current neighbor can't be sorted again
neighbors[new_neighbor].first <= closest_distance )
// if next neighbor is closer or as close as the closest branch, set next neighbor
if (! m_neighbors.empty() && m_neighbors.top().first <= closest_branch.distance )
{
current_neighbor = new_neighbor;
m_neighbor_ptr = m_neighbors.top().second;
++m_neighbors_count;
m_neighbors.pop_top();
return;
}
// if node is further than the furthest neighbor, following nodes will also be further
BOOST_GEOMETRY_INDEX_ASSERT(neighbors.size() <= max_count(), "unexpected neighbors count");
if ( max_count() <= neighbors.size() &&
neighbors.back().first <= closest_distance )
BOOST_GEOMETRY_INDEX_ASSERT(m_neighbors_count + m_neighbors.size() <= max_count(), "unexpected neighbors count");
// if there is enough neighbors and there is no closer branch
if (ignore_branch_or_value(closest_branch.distance))
{
internal_heap.clear();
m_branches.clear();
continue;
}
else
{
node_pointer ptr = closest_branch.second;
std::pop_heap(internal_heap.begin(), internal_heap.end(), pair_first_greater());
internal_heap.pop_back();
node_pointer ptr = closest_branch.ptr;
size_type reverse_level = closest_branch.reverse_level;
m_branches.pop();
rtree::apply_visitor(*this, *ptr);
apply(ptr, reverse_level);
}
}
}
@@ -428,112 +473,80 @@ public:
bool is_end() const
{
return (std::numeric_limits<size_type>::max)() == current_neighbor;
return m_neighbor_ptr == nullptr;
}
friend bool operator==(distance_query_incremental const& l, distance_query_incremental const& r)
{
BOOST_GEOMETRY_INDEX_ASSERT(l.current_neighbor != r.current_neighbor ||
(std::numeric_limits<size_type>::max)() == l.current_neighbor ||
(std::numeric_limits<size_type>::max)() == r.current_neighbor ||
l.neighbors[l.current_neighbor].second == r.neighbors[r.current_neighbor].second,
"not corresponding iterators");
return l.current_neighbor == r.current_neighbor;
return l.m_neighbors_count == r.m_neighbors_count;
}
// Put node's elements into the list of active branches if those elements meets predicates
// and distance predicates(currently not used)
// and aren't further than found neighbours (if there is enough neighbours)
inline void operator()(internal_node const& n)
private:
void apply(node_pointer ptr, size_type reverse_level)
{
typedef typename rtree::elements_type<internal_node>::type elements_type;
elements_type const& elements = rtree::elements(n);
// fill active branch list array of nodes meeting predicates
for ( typename elements_type::const_iterator it = elements.begin() ; it != elements.end() ; ++it )
namespace id = index::detail;
// Put node's elements into the list of active branches if those elements meets predicates
// and distance predicates(currently not used)
// and aren't further than found neighbours (if there is enough neighbours)
if (reverse_level > 0)
{
// if current node meets predicates
// 0 - dummy value
if ( index::detail::predicates_check
<
index::detail::bounds_tag, 0, predicates_len
>(m_pred, 0, it->first, m_strategy) )
internal_node& n = rtree::get<internal_node>(*ptr);
// fill active branch list array of nodes meeting predicates
for (auto const& p : rtree::elements(n))
{
// calculate node's distance(s) for distance predicate
node_distance_type node_distance;
// if distance isn't ok - move to the next node
if ( !calculate_node_distance::apply(predicate(), it->first,
m_strategy, node_distance) )
{
continue;
}
node_distance_type node_distance; // for distance predicate
// if current node is further than found neighbors - don't analyze it
if ( max_count() <= neighbors.size() &&
neighbors.back().first <= node_distance )
// if current node meets predicates (0 is dummy value)
if (id::predicates_check<id::bounds_tag>(m_pred, 0, p.first, m_strategy)
// and if distance is ok
&& calculate_node_distance::apply(predicate(), p.first, m_strategy, node_distance)
// and if current node is closer than the furthest neighbor
&& ! ignore_branch_or_value(node_distance))
{
continue;
// add current node into the queue
m_branches.push(branch_data(node_distance, reverse_level - 1, p.second));
}
// add current node's data into the queue
internal_heap.push_back(std::make_pair(node_distance, it->second));
std::push_heap(internal_heap.begin(), internal_heap.end(), pair_first_greater());
}
}
}
// Put values into the list of neighbours if those values meets predicates
// and distance predicates(currently not used)
// and aren't further than already found neighbours (if there is enough neighbours)
inline void operator()(leaf const& n)
{
typedef typename rtree::elements_type<leaf>::type elements_type;
elements_type const& elements = rtree::elements(n);
// store distance to the furthest neighbour
bool not_enough_neighbors = neighbors.size() < max_count();
value_distance_type greatest_distance = !not_enough_neighbors ? neighbors.back().first : (std::numeric_limits<value_distance_type>::max)();
// search leaf for closest value meeting predicates
for ( typename elements_type::const_iterator it = elements.begin() ; it != elements.end() ; ++it)
// Put values into the list of neighbours if those values meets predicates
// and distance predicates(currently not used)
// and aren't further than already found neighbours (if there is enough neighbours)
else
{
// if value meets predicates
if ( index::detail::predicates_check
<
index::detail::value_tag, 0, predicates_len
>(m_pred, *it, (*m_translator)(*it), m_strategy) )
leaf& n = rtree::get<leaf>(*ptr);
// search leaf for closest value meeting predicates
for (auto const& v : rtree::elements(n))
{
// calculate values distance for distance predicate
value_distance_type value_distance;
// if distance is ok
if ( calculate_value_distance::apply(predicate(), (*m_translator)(*it),
m_strategy, value_distance) )
value_distance_type value_distance; // for distance predicate
// if value meets predicates
if (id::predicates_check<id::value_tag>(m_pred, v, (*m_tr)(v), m_strategy)
// and if distance is ok
&& calculate_value_distance::apply(predicate(), (*m_tr)(v), m_strategy, value_distance)
// and if current value is closer than the furthest neighbor
&& ! ignore_branch_or_value(value_distance))
{
// if there is not enough values or current value is closer than furthest neighbour
if ( not_enough_neighbors || value_distance < greatest_distance )
// add current value into the queue
m_neighbors.push(std::make_pair(value_distance, boost::addressof(v)));
// remove unneeded value
if (m_neighbors_count + m_neighbors.size() > max_count())
{
neighbors.push_back(std::make_pair(value_distance, boost::addressof(*it)));
m_neighbors.pop_bottom();
}
}
}
}
// TODO: sort is probably suboptimal.
// An alternative would be std::set, but it'd probably add constant cost.
// Ideally replace this with double-ended priority queue, e.g. min-max heap.
// NOTE: A condition in increment() relies on the fact that current neighbor doesn't
// participate in sorting anymore.
// sort array
size_type sort_first = current_neighbor == (std::numeric_limits<size_type>::max)() ? 0 : current_neighbor + 1;
std::sort(neighbors.begin() + sort_first, neighbors.end(), pair_first_less());
// remove furthest values
if ( max_count() < neighbors.size() )
neighbors.resize(max_count());
}
private:
inline std::size_t max_count() const
template <typename Distance>
bool ignore_branch_or_value(Distance const& distance)
{
return m_neighbors_count + m_neighbors.size() == max_count()
&& (m_neighbors.empty() || m_neighbors.bottom().first <= distance);
}
std::size_t max_count() const
{
return nearest_predicate_access::get(m_pred).count;
}
@@ -543,15 +556,15 @@ private:
return nearest_predicate_access::get(m_pred);
}
const translator_type * m_translator;
const translator_type * m_tr;
strategy_type m_strategy;
Predicates m_pred;
internal_heap_type internal_heap;
std::vector< std::pair<value_distance_type, const value_type *> > neighbors;
size_type current_neighbor;
strategy_type m_strategy;
branches_type m_branches;
neighbors_type m_neighbors;
size_type m_neighbors_count;
const value_type * m_neighbor_ptr;
};
}}} // namespace detail::rtree::visitors

218
include/boost/geometry/index/detail/rtree/visitors/spatial_query.hpp
View File

@@ -25,7 +25,6 @@ namespace detail { namespace rtree { namespace visitors {
template <typename MembersHolder, typename Predicates, typename OutIter>
struct spatial_query
: public MembersHolder::visitor_const
{
typedef typename MembersHolder::parameters_type parameters_type;
typedef typename MembersHolder::translator_type translator_type;
@@ -37,71 +36,69 @@ struct spatial_query
typedef typename MembersHolder::internal_node internal_node;
typedef typename MembersHolder::leaf leaf;
typedef typename allocators_type::node_pointer node_pointer;
typedef typename allocators_type::size_type size_type;
static const std::size_t predicates_len = index::detail::predicates_length<Predicates>::value;
inline spatial_query(parameters_type const& par, translator_type const& t, Predicates const& p, OutIter out_it)
: tr(t), pred(p), out_iter(out_it), found_count(0), strategy(index::detail::get_strategy(par))
spatial_query(MembersHolder const& members, Predicates const& p, OutIter out_it)
: m_tr(members.translator())
, m_strategy(index::detail::get_strategy(members.parameters()))
, m_pred(p)
, m_out_iter(out_it)
, m_found_count(0)
{}
inline void operator()(internal_node const& n)
size_type apply(node_pointer ptr, size_type reverse_level)
{
typedef typename rtree::elements_type<internal_node>::type elements_type;
elements_type const& elements = rtree::elements(n);
// traverse nodes meeting predicates
for (typename elements_type::const_iterator it = elements.begin();
it != elements.end(); ++it)
namespace id = index::detail;
if (reverse_level > 0)
{
// if node meets predicates
// 0 - dummy value
if ( index::detail::predicates_check
<
index::detail::bounds_tag, 0, predicates_len
>(pred, 0, it->first, strategy) )
internal_node& n = rtree::get<internal_node>(*ptr);
// traverse nodes meeting predicates
for (auto const& p : rtree::elements(n))
{
rtree::apply_visitor(*this, *it->second);
// if node meets predicates (0 is dummy value)
if (id::predicates_check<id::bounds_tag>(m_pred, 0, p.first, m_strategy))
{
apply(p.second, reverse_level - 1);
}
}
}
}
inline void operator()(leaf const& n)
{
typedef typename rtree::elements_type<leaf>::type elements_type;
elements_type const& elements = rtree::elements(n);
// get all values meeting predicates
for (typename elements_type::const_iterator it = elements.begin();
it != elements.end(); ++it)
else
{
// if value meets predicates
if ( index::detail::predicates_check
<
index::detail::value_tag, 0, predicates_len
>(pred, *it, tr(*it), strategy) )
leaf& n = rtree::get<leaf>(*ptr);
// get all values meeting predicates
for (auto const& v : rtree::elements(n))
{
*out_iter = *it;
++out_iter;
++found_count;
// if value meets predicates
if (id::predicates_check<id::value_tag>(m_pred, v, m_tr(v), m_strategy))
{
*m_out_iter = v;
++m_out_iter;
++m_found_count;
}
}
}
return m_found_count;
}
translator_type const& tr;
size_type apply(MembersHolder const& members)
{
return apply(members.root, members.leafs_level);
}
Predicates pred;
private:
translator_type const& m_tr;
strategy_type m_strategy;
OutIter out_iter;
size_type found_count;
Predicates const& m_pred;
OutIter m_out_iter;
strategy_type strategy;
size_type m_found_count;
};
template <typename MembersHolder, typename Predicates>
class spatial_query_incremental
: public MembersHolder::visitor_const
{
typedef typename MembersHolder::value_type value_type;
typedef typename MembersHolder::parameters_type parameters_type;
@@ -110,7 +107,6 @@ class spatial_query_incremental
typedef typename index::detail::strategy_type<parameters_type>::type strategy_type;
public:
typedef typename MembersHolder::node node;
typedef typename MembersHolder::internal_node internal_node;
typedef typename MembersHolder::leaf leaf;
@@ -123,37 +119,40 @@ public:
typedef typename rtree::elements_type<leaf>::type leaf_elements;
typedef typename rtree::elements_type<leaf>::type::const_iterator leaf_iterator;
static const std::size_t predicates_len = index::detail::predicates_length<Predicates>::value;
struct internal_data
{
internal_data(internal_iterator f, internal_iterator l, size_type rl)
: first(f), last(l), reverse_level(rl)
{}
internal_iterator first;
internal_iterator last;
size_type reverse_level;
};
inline spatial_query_incremental()
: m_translator(NULL)
// , m_pred()
, m_values(NULL)
, m_current()
public:
spatial_query_incremental()
: m_translator(nullptr)
// , m_strategy()
{}
inline spatial_query_incremental(parameters_type const& params, translator_type const& t, Predicates const& p)
: m_translator(::boost::addressof(t))
, m_pred(p)
, m_values(NULL)
// , m_pred()
, m_values(nullptr)
, m_current()
, m_strategy(index::detail::get_strategy(params))
{}
inline void operator()(internal_node const& n)
{
typedef typename rtree::elements_type<internal_node>::type elements_type;
elements_type const& elements = rtree::elements(n);
spatial_query_incremental(Predicates const& p)
: m_translator(nullptr)
// , m_strategy()
, m_pred(p)
, m_values(nullptr)
, m_current()
{}
m_internal_stack.push_back(std::make_pair(elements.begin(), elements.end()));
}
inline void operator()(leaf const& n)
{
m_values = ::boost::addressof(rtree::elements(n));
m_current = rtree::elements(n).begin();
}
spatial_query_incremental(MembersHolder const& members, Predicates const& p)
: m_translator(::boost::addressof(members.translator()))
, m_strategy(index::detail::get_strategy(members.parameters()))
, m_pred(p)
, m_values(nullptr)
, m_current()
{}
const_reference dereference() const
{
@@ -161,9 +160,9 @@ public:
return *m_current;
}
void initialize(node_pointer root)
void initialize(MembersHolder const& members)
{
rtree::apply_visitor(*this, *root);
apply(members.root, members.leafs_level);
search_value();
}
@@ -173,8 +172,38 @@ public:
search_value();
}
bool is_end() const
{
return 0 == m_values;
}
friend bool operator==(spatial_query_incremental const& l, spatial_query_incremental const& r)
{
return (l.m_values == r.m_values) && (0 == l.m_values || l.m_current == r.m_current);
}
private:
void apply(node_pointer ptr, size_type reverse_level)
{
namespace id = index::detail;
if (reverse_level > 0)
{
internal_node& n = rtree::get<internal_node>(*ptr);
auto const& elements = rtree::elements(n);
m_internal_stack.push_back(internal_data(elements.begin(), elements.end(), reverse_level - 1));
}
else
{
leaf& n = rtree::get<leaf>(*ptr);
m_values = ::boost::addressof(rtree::elements(n));
m_current = rtree::elements(n).begin();
}
}
void search_value()
{
namespace id = index::detail;
for (;;)
{
// if leaf is choosen, move to the next value in leaf
@@ -184,10 +213,7 @@ public:
{
// return if next value is found
value_type const& v = *m_current;
if (index::detail::predicates_check
<
index::detail::value_tag, 0, predicates_len
>(m_pred, v, (*m_translator)(v), m_strategy))
if (id::predicates_check<id::value_tag>(m_pred, v, (*m_translator)(v), m_strategy))
{
return;
}
@@ -204,52 +230,40 @@ public:
else
{
// return if there is no more nodes to traverse
if ( m_internal_stack.empty() )
if (m_internal_stack.empty())
{
return;
}
// no more children in current node, remove it from stack
if ( m_internal_stack.back().first == m_internal_stack.back().second )
internal_data& current_data = m_internal_stack.back();
// no more children in current node, remove it from stack
if (current_data.first == current_data.last)
{
m_internal_stack.pop_back();
continue;
}
internal_iterator it = m_internal_stack.back().first;
++m_internal_stack.back().first;
internal_iterator it = current_data.first;
++current_data.first;
// next node is found, push it to the stack
if (index::detail::predicates_check
<
index::detail::bounds_tag, 0, predicates_len
>(m_pred, 0, it->first, m_strategy))
if (id::predicates_check<id::bounds_tag>(m_pred, 0, it->first, m_strategy))
{
rtree::apply_visitor(*this, *(it->second));
apply(it->second, current_data.reverse_level);
}
}
}
}
bool is_end() const
{
return 0 == m_values;
}
friend bool operator==(spatial_query_incremental const& l, spatial_query_incremental const& r)
{
return (l.m_values == r.m_values) && (0 == l.m_values || l.m_current == r.m_current );
}
private:
const translator_type * m_translator;
strategy_type m_strategy;
Predicates m_pred;
std::vector< std::pair<internal_iterator, internal_iterator> > m_internal_stack;
std::vector<internal_data> m_internal_stack;
const leaf_elements * m_values;
leaf_iterator m_current;
strategy_type m_strategy;
};
}}} // namespace detail::rtree::visitors

131
include/boost/geometry/index/rtree.hpp
View File

@@ -1079,17 +1079,9 @@ public:
template <typename Predicates, typename OutIter>
size_type query(Predicates const& predicates, OutIter out_it) const
{
if ( !m_members.root )
return 0;
static const std::size_t distance_predicates_count = detail::predicates_count_distance<Predicates>::value;
static const bool is_distance_predicate = 0 < distance_predicates_count;
BOOST_GEOMETRY_STATIC_ASSERT((distance_predicates_count <= 1),
"Only one distance predicate can be passed.",
Predicates);
return query_dispatch(predicates, out_it,
std::integral_constant<bool, is_distance_predicate>());
return m_members.root
? query_dispatch(predicates, out_it)
: 0;
}
/*!
@@ -1183,6 +1175,15 @@ public:
return const_query_iterator();
}
private:
template <typename Predicates>
using query_iterator_t = std::conditional_t
<
detail::predicates_count_distance<Predicates>::value == 0,
detail::rtree::iterators::spatial_query_iterator<members_holder, Predicates>,
detail::rtree::iterators::distance_query_iterator<members_holder, Predicates>
>;
#ifndef BOOST_GEOMETRY_INDEX_DETAIL_EXPERIMENTAL
private:
#endif
@@ -1240,38 +1241,15 @@ private:
\return The iterator pointing at the begin of the query range.
*/
template <typename Predicates>
std::conditional_t
<
detail::predicates_count_distance<Predicates>::value == 0,
detail::rtree::iterators::spatial_query_iterator<members_holder, Predicates>,
detail::rtree::iterators::distance_query_iterator
<
members_holder, Predicates,
detail::predicates_find_distance<Predicates>::value
>
>
qbegin_(Predicates const& predicates) const
query_iterator_t<Predicates> qbegin_(Predicates const& predicates) const
{
static const std::size_t distance_predicates_count = detail::predicates_count_distance<Predicates>::value;
BOOST_GEOMETRY_STATIC_ASSERT((distance_predicates_count <= 1),
BOOST_GEOMETRY_STATIC_ASSERT((detail::predicates_count_distance<Predicates>::value <= 1),
"Only one distance predicate can be passed.",
Predicates);
typedef std::conditional_t
<
detail::predicates_count_distance<Predicates>::value == 0,
detail::rtree::iterators::spatial_query_iterator<members_holder, Predicates>,
detail::rtree::iterators::distance_query_iterator
<
members_holder, Predicates,
detail::predicates_find_distance<Predicates>::value
>
> iterator_type;
if ( !m_members.root )
return iterator_type(m_members.parameters(), m_members.translator(), predicates);
return iterator_type(m_members.root, m_members.parameters(), m_members.translator(), predicates);
return m_members.root
? query_iterator_t<Predicates>(m_members, predicates)
: query_iterator_t<Predicates>(predicates);
}
/*!
@@ -1307,35 +1285,13 @@ private:
\return The iterator pointing at the end of the query range.
*/
template <typename Predicates>
std::conditional_t
<
detail::predicates_count_distance<Predicates>::value == 0,
detail::rtree::iterators::spatial_query_iterator<members_holder, Predicates>,
detail::rtree::iterators::distance_query_iterator
<
members_holder, Predicates,
detail::predicates_find_distance<Predicates>::value
>
>
qend_(Predicates const& predicates) const
query_iterator_t<Predicates> qend_(Predicates const& predicates) const
{
static const std::size_t distance_predicates_count = detail::predicates_count_distance<Predicates>::value;
BOOST_GEOMETRY_STATIC_ASSERT((distance_predicates_count <= 1),
BOOST_GEOMETRY_STATIC_ASSERT((detail::predicates_count_distance<Predicates>::value <= 1),
"Only one distance predicate can be passed.",
Predicates);
typedef std::conditional_t
<
detail::predicates_count_distance<Predicates>::value == 0,
detail::rtree::iterators::spatial_query_iterator<members_holder, Predicates>,
detail::rtree::iterators::distance_query_iterator
<
members_holder, Predicates,
detail::predicates_find_distance<Predicates>::value
>
> iterator_type;
return iterator_type(m_members.parameters(), m_members.translator(), predicates);
return query_iterator_t<Predicates>(m_members.parameters(), m_members.translator(), predicates);
}
/*!
@@ -1436,10 +1392,9 @@ public:
*/
const_iterator begin() const
{
if ( !m_members.root )
return const_iterator();
return const_iterator(m_members.root);
return m_members.root
? const_iterator(m_members.root)
: const_iterator();
}
/*!
@@ -1894,15 +1849,16 @@ private:
\par Exception-safety
strong
*/
template <typename Predicates, typename OutIter>
size_type query_dispatch(Predicates const& predicates, OutIter out_it, std::false_type /*is_distance_predicate*/) const
template
<
typename Predicates, typename OutIter,
std::enable_if_t<(detail::predicates_count_distance<Predicates>::value == 0), int> = 0
>
size_type query_dispatch(Predicates const& predicates, OutIter out_it) const
{
detail::rtree::visitors::spatial_query<members_holder, Predicates, OutIter>
find_v(m_members.parameters(), m_members.translator(), predicates, out_it);
detail::rtree::apply_visitor(find_v, *m_members.root);
return find_v.found_count;
query(m_members, predicates, out_it);
return query.apply(m_members);
}
/*!
@@ -1911,22 +1867,21 @@ private:
\par Exception-safety
strong
*/
template <typename Predicates, typename OutIter>
size_type query_dispatch(Predicates const& predicates, OutIter out_it, std::true_type /*is_distance_predicate*/) const
template
<
typename Predicates, typename OutIter,
std::enable_if_t<(detail::predicates_count_distance<Predicates>::value > 0), int> = 0
>
size_type query_dispatch(Predicates const& predicates, OutIter out_it) const
{
BOOST_GEOMETRY_INDEX_ASSERT(m_members.root, "The root must exist");
BOOST_GEOMETRY_STATIC_ASSERT((detail::predicates_count_distance<Predicates>::value == 1),
"Only one distance predicate can be passed.",
Predicates);
static const std::size_t distance_predicate_index = detail::predicates_find_distance<Predicates>::value;
detail::rtree::visitors::distance_query<
members_holder,
Predicates,
distance_predicate_index,
OutIter
> distance_v(m_members.parameters(), m_members.translator(), predicates, out_it);
detail::rtree::visitors::distance_query<members_holder, Predicates>
distance_v(m_members, predicates);
detail::rtree::apply_visitor(distance_v, *m_members.root);
return distance_v.finish();
return distance_v.apply(m_members, out_it);
}
/*!

2
include/boost/geometry/policies/robustness/get_rescale_policy.hpp
View File

@@ -24,6 +24,7 @@
#include <boost/geometry/core/assert.hpp>
#include <boost/geometry/core/config.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/tag_cast.hpp>
#include <boost/geometry/algorithms/envelope.hpp>
@@ -37,7 +38,6 @@
#include <boost/geometry/policies/robustness/no_rescale_policy.hpp>
#include <boost/geometry/policies/robustness/rescale_policy.hpp>
#include <boost/geometry/policies/robustness/robust_type.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/type_traits.hpp>
// TEMP

66
include/boost/geometry/policies/robustness/segment_ratio.hpp
View File

@@ -19,8 +19,8 @@
#include <boost/rational.hpp>
#include <boost/geometry/core/assert.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
namespace boost { namespace geometry
{
@@ -135,11 +135,15 @@ struct possibly_collinear<Type, false>
template <typename Type>
class segment_ratio
{
public :
typedef Type numeric_type;
// Type used for the approximation (a helper value)
// and for the edge value (0..1) (a helper function).
using floating_point_type =
typename detail::promoted_to_floating_point<Type>::type;
// Type-alias for the type itself
typedef segment_ratio<Type> thistype;
using thistype = segment_ratio<Type>;
public :
inline segment_ratio()
: m_numerator(0)
@@ -221,8 +225,8 @@ public :
m_approximation =
m_denominator == 0 ? 0
: (
boost::numeric_cast<fp_type>(m_numerator) * scale()
/ boost::numeric_cast<fp_type>(m_denominator)
boost::numeric_cast<floating_point_type>(m_numerator) * scale()
/ boost::numeric_cast<floating_point_type>(m_denominator)
);
}
@@ -254,21 +258,16 @@ public :
return m_numerator > m_denominator;
}
inline bool near_end() const
//! Returns a value between 0.0 and 1.0
//! 0.0 means: exactly in the middle
//! 1.0 means: exactly on one of the edges (or even over it)
inline floating_point_type edge_value() const
{
if (left() || right())
{
return false;
}
static fp_type const small_part_of_scale = scale() / 100;
return m_approximation < small_part_of_scale
|| m_approximation > scale() - small_part_of_scale;
}
inline bool close_to(thistype const& other) const
{
return geometry::math::abs(m_approximation - other.m_approximation) < 50;
using fp = floating_point_type;
fp const one{1.0};
floating_point_type const result
= fp(2) * geometry::math::abs(fp(0.5) - m_approximation / scale());
return result > one ? one : result;
}
template <typename Threshold>
@@ -313,19 +312,7 @@ public :
}
#endif
private :
// NOTE: if this typedef is used then fp_type is non-fundamental type
// if Type is non-fundamental type
//typedef typename promote_floating_point<Type>::type fp_type;
// TODO: What with user-defined numeric types?
// Shouldn't here is_integral be checked?
typedef std::conditional_t
<
std::is_floating_point<Type>::value, Type, double
> fp_type;
Type m_numerator;
Type m_denominator;
@@ -335,12 +322,19 @@ private :
// Boost.Rational is used if the approximations are close.
// Reason: performance, Boost.Rational does a GCD by default and also the
// comparisons contain while-loops.
fp_type m_approximation;
floating_point_type m_approximation;
static inline fp_type scale()
inline bool close_to(thistype const& other) const
{
return 1000000.0;
static floating_point_type const threshold{50.0};
return geometry::math::abs(m_approximation - other.m_approximation)
< threshold;
}
static inline floating_point_type scale()
{
static floating_point_type const fp_scale{1000000.0};
return fp_scale;
}
};

10
include/boost/geometry/strategies/agnostic/point_in_box_by_side.hpp
View File

@@ -20,12 +20,16 @@
#include <boost/array.hpp>
#include <boost/core/ignore_unused.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_dimension.hpp>
#include <boost/geometry/algorithms/assign.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/covered_by.hpp>
#include <boost/geometry/strategies/within.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategies/geographic/side.hpp>
#include <boost/geometry/strategies/spherical/ssf.hpp>
@@ -118,7 +122,7 @@ struct cartesian_point_box_by_side
<
within::detail::decide_within
>(point, box,
strategy::side::side_by_triangle<CalculationType>());
strategy::side::side_robust<CalculationType>());
}
};
@@ -189,7 +193,7 @@ struct cartesian_point_box_by_side
<
within::detail::decide_covered_by
>(point, box,
strategy::side::side_by_triangle<CalculationType>());
strategy::side::side_robust<CalculationType>());
}
};

2
include/boost/geometry/strategies/cartesian/azimuth.hpp
View File

@@ -14,10 +14,10 @@
#include <cmath>
#include <boost/geometry/core/tags.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/strategies/azimuth.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_most_precise.hpp>
namespace boost { namespace geometry

1
include/boost/geometry/strategies/cartesian/distance_segment_box.hpp
View File

@@ -17,7 +17,6 @@
#include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
#include <boost/geometry/strategies/cartesian/distance_pythagoras_point_box.hpp>
#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
namespace boost { namespace geometry
{

89
include/boost/geometry/strategies/cartesian/intersection.hpp
View File

@@ -38,13 +38,13 @@
#include <boost/geometry/strategy/cartesian/envelope.hpp>
#include <boost/geometry/strategy/cartesian/expand_box.hpp>
#include <boost/geometry/strategy/cartesian/expand_segment.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/cartesian/disjoint_box_box.hpp>
#include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp>
#include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategies/covered_by.hpp>
#include <boost/geometry/strategies/intersection.hpp>
#include <boost/geometry/strategies/intersection_result.hpp>
@@ -68,6 +68,61 @@ namespace boost { namespace geometry
namespace strategy { namespace intersection
{
namespace detail_usage
{
// When calculating the intersection, the information of "a" or "b" can be used.
// Theoretically this gives equal results, but due to floating point precision
// there might be tiny differences. These are edge cases.
// This structure is to determine if "a" or "b" should be used.
// Prefer the segment closer to the endpoint.
// If both are about equally close, then prefer the longer segment
// To avoid hard thresholds, behavior is made fluent.
// Calculate comparable length indications,
// the longer the segment (relatively), the lower the value
// such that the shorter lengths are evaluated higher and will
// be preferred.
template <bool IsArithmetic>
struct use_a
{
template <typename Ct, typename Ev>
static bool apply(Ct const& cla, Ct const& clb, Ev const& eva, Ev const& evb)
{
auto const clm = (std::max)(cla, clb);
if (clm <= 0)
{
return true;
}
// Relative comparible length
auto const rcla = Ct(1.0) - cla / clm;
auto const rclb = Ct(1.0) - clb / clm;
// Multipliers for edgevalue (ev) and relative comparible length (rcl)
// They determine the balance between edge value (should be larger)
// and segment length. In 99.9xx% of the cases there is no difference
// at all (if either a or b is used). Therefore the values of the
// constants are not sensitive for the majority of the situations.
// One known case is #mysql_23023665_6 (difference) which needs mev >= 2
Ev const mev = 5;
Ev const mrcl = 1;
return mev * eva + mrcl * rcla > mev * evb + mrcl * rclb;
}
};
// Specialization for non arithmetic types. They will always use "a"
template <>
struct use_a<false>
{
template <typename Ct, typename Ev>
static bool apply(Ct const& , Ct const& , Ev const& , Ev const& )
{
return true;
}
};
}
/*!
\see http://mathworld.wolfram.com/Line-LineIntersection.html
@@ -180,30 +235,12 @@ struct cartesian_segments
template <typename Point, typename Segment1, typename Segment2>
void calculate(Point& point, Segment1 const& a, Segment2 const& b) const
{
bool use_a = true;
// Prefer one segment if one is on or near an endpoint
bool const a_near_end = robust_ra.near_end();
bool const b_near_end = robust_rb.near_end();
if (a_near_end && ! b_near_end)
{
use_a = true;
}
else if (b_near_end && ! a_near_end)
{
use_a = false;
}
else
{
// Prefer shorter segment
promoted_type const len_a = comparable_length_a();
promoted_type const len_b = comparable_length_b();
if (len_b < len_a)
{
use_a = false;
}
// else use_a is true but was already assigned like that
}
bool const use_a
= detail_usage::use_a
<
std::is_arithmetic<CoordinateType>::value
>::apply(comparable_length_a(), comparable_length_b(),
robust_ra.edge_value(), robust_rb.edge_value());
if (use_a)
{
@@ -402,7 +439,7 @@ struct cartesian_segments
return Policy::disjoint();
}
typedef side::side_by_triangle<CalculationType> side_strategy_type;
typedef side::side_robust<CalculationType> side_strategy_type;
side_info sides;
sides.set<0>(side_strategy_type::apply(q1, q2, p1),
side_strategy_type::apply(q1, q2, p2));

4
include/boost/geometry/strategies/cartesian/point_in_poly_winding.hpp
View File

@@ -24,10 +24,10 @@
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/strategy/cartesian/expand_point.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/cartesian/point_in_box.hpp>
#include <boost/geometry/strategies/cartesian/disjoint_box_box.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategies/covered_by.hpp>
#include <boost/geometry/strategies/within.hpp>
@@ -116,7 +116,7 @@ public:
else // count == 2 || count == -2
{
// 1 left, -1 right
typedef side::side_by_triangle<CalculationType> side_strategy_type;
typedef side::side_robust<CalculationType> side_strategy_type;
side = side_strategy_type::apply(s1, s2, point);
}

270
include/boost/geometry/strategies/cartesian/side_by_triangle.hpp
View File

@@ -1,275 +1,21 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Boost.Geometry
// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
// Copyright (c) 2021, Oracle and/or its affiliates.
// This file was modified by Oracle on 2015-2021.
// Modifications copyright (c) 2015-2021, Oracle and/or its affiliates.
// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
// Licensed under the Boost Software License version 1.0.
// http://www.boost.org/users/license.html
#ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_BY_TRIANGLE_HPP
#define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_BY_TRIANGLE_HPP
#include <type_traits>
#include <boost/geometry/arithmetic/determinant.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp>
#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
#include <boost/geometry/strategies/compare.hpp>
#include <boost/geometry/strategies/side.hpp>
#include <boost/geometry/strategy/cartesian/envelope.hpp>
#include <boost/geometry/util/select_most_precise.hpp>
#include <boost/config/pragma_message.hpp>
BOOST_PRAGMA_MESSAGE("This include file is deprecated and will be removed in the future.")
namespace boost { namespace geometry
{
namespace strategy { namespace side
{
/*!
\brief Check at which side of a segment a point lies:
left of segment (> 0), right of segment (< 0), on segment (0)
\ingroup strategies
\tparam CalculationType \tparam_calculation
*/
template <typename CalculationType = void>
class side_by_triangle
{
template <typename Policy>
struct eps_policy
{
eps_policy() {}
template <typename Type>
eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
: policy(a, b, c, d)
{}
Policy policy;
};
struct eps_empty
{
eps_empty() {}
template <typename Type>
eps_empty(Type const&, Type const&, Type const&, Type const&) {}
};
public :
typedef cartesian_tag cs_tag;
// Template member function, because it is not always trivial
// or convenient to explicitly mention the typenames in the
// strategy-struct itself.
// Types can be all three different. Therefore it is
// not implemented (anymore) as "segment"
template
<
typename CoordinateType,
typename PromotedType,
typename P1,
typename P2,
typename P,
typename EpsPolicy
>
static inline
PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
{
CoordinateType const x = get<0>(p);
CoordinateType const y = get<1>(p);
CoordinateType const sx1 = get<0>(p1);
CoordinateType const sy1 = get<1>(p1);
CoordinateType const sx2 = get<0>(p2);
CoordinateType const sy2 = get<1>(p2);
PromotedType const dx = sx2 - sx1;
PromotedType const dy = sy2 - sy1;
PromotedType const dpx = x - sx1;
PromotedType const dpy = y - sy1;
eps_policy = EpsPolicy(dx, dy, dpx, dpy);
return geometry::detail::determinant<PromotedType>
(
dx, dy,
dpx, dpy
);
}
template
<
typename CoordinateType,
typename PromotedType,
typename P1,
typename P2,
typename P
>
static inline
PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
{
eps_empty dummy;
return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
}
template
<
typename CoordinateType,
typename PromotedType,
bool AreAllIntegralCoordinates
>
struct compute_side_value
{
template <typename P1, typename P2, typename P, typename EpsPolicy>
static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
{
return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
}
};
template <typename CoordinateType, typename PromotedType>
struct compute_side_value<CoordinateType, PromotedType, false>
{
template <typename P1, typename P2, typename P, typename EpsPolicy>
static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
{
// For robustness purposes, first check if any two points are
// the same; in this case simply return that the points are
// collinear
if (equals_point_point(p1, p2)
|| equals_point_point(p1, p)
|| equals_point_point(p2, p))
{
return PromotedType(0);
}
// The side_by_triangle strategy computes the signed area of
// the point triplet (p1, p2, p); as such it is (in theory)
// invariant under cyclic permutations of its three arguments.
//
// In the context of numerical errors that arise in
// floating-point computations, and in order to make the strategy
// consistent with respect to cyclic permutations of its three
// arguments, we cyclically permute them so that the first
// argument is always the lexicographically smallest point.
typedef compare::cartesian<compare::less> less;
if (less::apply(p, p1))
{
if (less::apply(p, p2))
{
// p is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
}
else
{
// p2 is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
}
}
if (less::apply(p1, p2))
{
// p1 is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
}
else
{
// p2 is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
}
}
};
template <typename P1, typename P2, typename P>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
typedef typename coordinate_type<P1>::type coordinate_type1;
typedef typename coordinate_type<P2>::type coordinate_type2;
typedef typename coordinate_type<P>::type coordinate_type3;
typedef std::conditional_t
<
std::is_void<CalculationType>::value,
typename select_most_precise
<
coordinate_type1,
coordinate_type2,
coordinate_type3
>::type,
CalculationType
> coordinate_type;
// Promote float->double, small int->int
typedef typename select_most_precise
<
coordinate_type,
double
>::type promoted_type;
bool const are_all_integral_coordinates =
std::is_integral<coordinate_type1>::value
&& std::is_integral<coordinate_type2>::value
&& std::is_integral<coordinate_type3>::value;
eps_policy< math::detail::equals_factor_policy<promoted_type> > epsp;
promoted_type s = compute_side_value
<
coordinate_type, promoted_type, are_all_integral_coordinates
>::apply(p1, p2, p, epsp);
promoted_type const zero = promoted_type();
return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
: s > zero ? 1
: -1;
}
private:
template <typename P1, typename P2>
static inline bool equals_point_point(P1 const& p1, P2 const& p2)
{
return strategy::within::cartesian_point_point::apply(p1, p2);
}
};
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
namespace services
{
template <typename CalculationType>
struct default_strategy<cartesian_tag, CalculationType>
{
typedef side_by_triangle<CalculationType> type;
};
}
#endif
}} // namespace strategy::side
}} // namespace boost::geometry
#include <boost/geometry/strategy/cartesian/side_by_triangle.hpp>
#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_SIDE_BY_TRIANGLE_HPP

2
include/boost/geometry/strategies/geographic/distance.hpp
View File

@@ -17,6 +17,7 @@
#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_HPP
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/radius.hpp>
@@ -32,7 +33,6 @@
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/geometries/point_xy.hpp>

2
include/boost/geometry/strategies/geographic/distance_cross_track.hpp
View File

@@ -21,6 +21,7 @@
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/tags.hpp>
@@ -37,7 +38,6 @@
#include <boost/geometry/formulas/vincenty_direct.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>

3
include/boost/geometry/strategies/geographic/distance_segment_box.hpp
View File

@@ -13,6 +13,8 @@
#include <boost/geometry/srs/spheroid.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/algorithms/detail/distance/segment_to_box.hpp>
#include <boost/geometry/strategies/distance.hpp>
@@ -26,7 +28,6 @@
#include <boost/geometry/strategies/spherical/distance_segment_box.hpp>
#include <boost/geometry/strategies/spherical/point_in_point.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
namespace boost { namespace geometry

2
include/boost/geometry/strategies/geographic/mapping_ssf.hpp
View File

@@ -17,10 +17,10 @@
#include <boost/core/ignore_unused.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/radius.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/strategies/side.hpp>

2
include/boost/geometry/strategies/geographic/side.hpp
View File

@@ -16,6 +16,7 @@
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/radius.hpp>
@@ -32,7 +33,6 @@
#include <boost/geometry/strategy/geographic/envelope.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
namespace boost { namespace geometry

4
include/boost/geometry/strategies/intersection_strategies.hpp
View File

@@ -19,12 +19,12 @@
#include <boost/geometry/policies/relate/intersection_policy.hpp>
#include <boost/geometry/policies/robustness/segment_ratio_type.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategies/intersection.hpp>
#include <boost/geometry/strategies/intersection_result.hpp>
#include <boost/geometry/strategies/side.hpp>
#include <boost/geometry/strategies/cartesian/intersection.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategies/spherical/intersection.hpp>
#include <boost/geometry/strategies/spherical/ssf.hpp>

2
include/boost/geometry/strategies/is_convex/cartesian.hpp
View File

@@ -11,10 +11,10 @@
#define BOOST_GEOMETRY_STRATEGIES_IS_CONVEX_CARTESIAN_HPP
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategies/convex_hull/cartesian.hpp>
#include <boost/geometry/strategies/detail.hpp>
#include <boost/geometry/strategies/is_convex/services.hpp>
#include <boost/geometry/strategy/cartesian/side_by_triangle.hpp>
namespace boost { namespace geometry

13
include/boost/geometry/strategies/relate/cartesian.hpp
View File

@@ -26,6 +26,8 @@
#include <boost/geometry/strategies/detail.hpp>
#include <boost/geometry/strategy/cartesian/area.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategy/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/area_box.hpp>
#include <boost/geometry/util/type_traits.hpp>
@@ -156,7 +158,7 @@ public:
static auto side()
{
return strategy::side::side_by_triangle<CalculationType>();
return strategy::side::side_robust<CalculationType>();
}
// within
@@ -381,6 +383,15 @@ struct strategy_converter<strategy::side::side_by_triangle<CalculationType>>
}
};
template <typename CalculationType>
struct strategy_converter<strategy::side::side_robust<CalculationType>>
{
static auto get(strategy::side::side_robust<CalculationType> const&)
{
return strategies::relate::cartesian<CalculationType>();
}
};
} // namespace services

2
include/boost/geometry/strategies/spherical/distance_cross_track.hpp
View File

@@ -24,6 +24,7 @@
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/tags.hpp>
@@ -36,7 +37,6 @@
#include <boost/geometry/strategies/spherical/intersection.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK

2
include/boost/geometry/strategies/spherical/distance_haversine.hpp
View File

@@ -17,6 +17,7 @@
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/radian_access.hpp>
@@ -26,7 +27,6 @@
#include <boost/geometry/strategies/spherical/get_radius.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>

1
include/boost/geometry/strategies/spherical/point_in_point.hpp
View File

@@ -28,6 +28,7 @@
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/coordinate_dimension.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/coordinate_system.hpp>
#include <boost/geometry/core/coordinate_type.hpp>
#include <boost/geometry/core/cs.hpp>

2
include/boost/geometry/strategies/spherical/side_by_cross_track.hpp
View File

@@ -16,6 +16,7 @@
#define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_SIDE_BY_CROSS_TRACK_HPP
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/radian_access.hpp>
@@ -26,7 +27,6 @@
#include <boost/geometry/strategies/spherical/point_in_point.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
namespace boost { namespace geometry

2
include/boost/geometry/strategies/spherical/ssf.hpp
View File

@@ -16,10 +16,10 @@
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/strategy/spherical/envelope.hpp>

3
include/boost/geometry/strategies/strategies.hpp
View File

@@ -67,7 +67,6 @@
#include <boost/geometry/strategies/cartesian/point_in_poly_crossings_multiply.hpp>
#include <boost/geometry/strategies/cartesian/point_in_poly_winding.hpp>
#include <boost/geometry/strategies/cartesian/line_interpolate.hpp>
#include <boost/geometry/strategies/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategies/spherical/azimuth.hpp>
#include <boost/geometry/strategies/spherical/densify.hpp>
@@ -135,6 +134,8 @@
#include <boost/geometry/strategy/cartesian/expand_box.hpp>
#include <boost/geometry/strategy/cartesian/expand_point.hpp>
#include <boost/geometry/strategy/cartesian/expand_segment.hpp>
#include <boost/geometry/strategy/cartesian/side_by_triangle.hpp>
#include <boost/geometry/strategy/cartesian/side_robust.hpp>
#include <boost/geometry/strategy/geographic/area.hpp>
#include <boost/geometry/strategy/geographic/envelope.hpp>

2
include/boost/geometry/strategies/strategy_transform.hpp
View File

@@ -30,10 +30,10 @@
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/radian_access.hpp>
#include <boost/geometry/core/coordinate_dimension.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/strategies/transform.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_coordinate_type.hpp>
namespace boost { namespace geometry

2
include/boost/geometry/strategies/transform/matrix_transformers.hpp
View File

@@ -33,9 +33,9 @@
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/core/coordinate_dimension.hpp>
#include <boost/geometry/core/coordinate_promotion.hpp>
#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/util/math.hpp>
#include <boost/geometry/util/promote_floating_point.hpp>
#include <boost/geometry/util/select_coordinate_type.hpp>
#include <boost/geometry/util/select_most_precise.hpp>

262
include/boost/geometry/strategy/cartesian/side_by_triangle.hpp Normal file
View File

@@ -0,0 +1,262 @@
// Boost.Geometry (aka GGL, Generic Geometry Library)
// Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
// This file was modified by Oracle on 2015-2021.
// Modifications copyright (c) 2015-2021, Oracle and/or its affiliates.
// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
// Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
// Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
// (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
#ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP
#define BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP
#include <type_traits>
#include <boost/geometry/arithmetic/determinant.hpp>
#include <boost/geometry/core/access.hpp>
#include <boost/geometry/strategies/cartesian/disjoint_segment_box.hpp>
#include <boost/geometry/strategies/cartesian/point_in_point.hpp>
#include <boost/geometry/strategies/compare.hpp>
#include <boost/geometry/strategies/side.hpp>
#include <boost/geometry/strategy/cartesian/envelope.hpp>
#include <boost/geometry/util/select_most_precise.hpp>
namespace boost { namespace geometry
{
namespace strategy { namespace side
{
/*!
\brief Check at which side of a segment a point lies:
left of segment (> 0), right of segment (< 0), on segment (0)
\ingroup strategies
\tparam CalculationType \tparam_calculation
*/
template <typename CalculationType = void>
class side_by_triangle
{
template <typename Policy>
struct eps_policy
{
eps_policy() {}
template <typename Type>
eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
: policy(a, b, c, d)
{}
Policy policy;
};
struct eps_empty
{
eps_empty() {}
template <typename Type>
eps_empty(Type const&, Type const&, Type const&, Type const&) {}
};
public :
typedef cartesian_tag cs_tag;
// Template member function, because it is not always trivial
// or convenient to explicitly mention the typenames in the
// strategy-struct itself.
// Types can be all three different. Therefore it is
// not implemented (anymore) as "segment"
template
<
typename CoordinateType,
typename PromotedType,
typename P1,
typename P2,
typename P,
typename EpsPolicy
>
static inline
PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
{
CoordinateType const x = get<0>(p);
CoordinateType const y = get<1>(p);
CoordinateType const sx1 = get<0>(p1);
CoordinateType const sy1 = get<1>(p1);
CoordinateType const sx2 = get<0>(p2);
CoordinateType const sy2 = get<1>(p2);
PromotedType const dx = sx2 - sx1;
PromotedType const dy = sy2 - sy1;
PromotedType const dpx = x - sx1;
PromotedType const dpy = y - sy1;
eps_policy = EpsPolicy(dx, dy, dpx, dpy);
return geometry::detail::determinant<PromotedType>
(
dx, dy,
dpx, dpy
);
}
template
<
typename CoordinateType,
typename PromotedType,
typename P1,
typename P2,
typename P
>
static inline
PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
{
eps_empty dummy;
return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
}
template
<
typename CoordinateType,
typename PromotedType,
bool AreAllIntegralCoordinates
>
struct compute_side_value
{
template <typename P1, typename P2, typename P, typename EpsPolicy>
static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
{
return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
}
};
template <typename CoordinateType, typename PromotedType>
struct compute_side_value<CoordinateType, PromotedType, false>
{
template <typename P1, typename P2, typename P, typename EpsPolicy>
static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
{
// For robustness purposes, first check if any two points are
// the same; in this case simply return that the points are
// collinear
if (equals_point_point(p1, p2)
|| equals_point_point(p1, p)
|| equals_point_point(p2, p))
{
return PromotedType(0);
}
// The side_by_triangle strategy computes the signed area of
// the point triplet (p1, p2, p); as such it is (in theory)
// invariant under cyclic permutations of its three arguments.
//
// In the context of numerical errors that arise in
// floating-point computations, and in order to make the strategy
// consistent with respect to cyclic permutations of its three
// arguments, we cyclically permute them so that the first
// argument is always the lexicographically smallest point.
typedef compare::cartesian<compare::less> less;
if (less::apply(p, p1))
{
if (less::apply(p, p2))
{
// p is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
}
else
{
// p2 is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
}
}
if (less::apply(p1, p2))
{
// p1 is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
}
else
{
// p2 is the lexicographically smallest
return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
}
}
};
template <typename P1, typename P2, typename P>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
typedef typename coordinate_type<P1>::type coordinate_type1;
typedef typename coordinate_type<P2>::type coordinate_type2;
typedef typename coordinate_type<P>::type coordinate_type3;
typedef std::conditional_t
<
std::is_void<CalculationType>::value,
typename select_most_precise
<
coordinate_type1,
coordinate_type2,
coordinate_type3
>::type,
CalculationType
> coordinate_type;
// Promote float->double, small int->int
typedef typename select_most_precise
<
coordinate_type,
double
>::type promoted_type;
bool const are_all_integral_coordinates =
std::is_integral<coordinate_type1>::value
&& std::is_integral<coordinate_type2>::value
&& std::is_integral<coordinate_type3>::value;
eps_policy< math::detail::equals_factor_policy<promoted_type> > epsp;
promoted_type s = compute_side_value
<
coordinate_type, promoted_type, are_all_integral_coordinates
>::apply(p1, p2, p, epsp);
promoted_type const zero = promoted_type();
return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
: s > zero ? 1
: -1;
}
private:
template <typename P1, typename P2>
static inline bool equals_point_point(P1 const& p1, P2 const& p2)
{
return strategy::within::cartesian_point_point::apply(p1, p2);
}
};
}} // namespace strategy::side
}} // namespace boost::geometry
#endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP

49
include/boost/geometry/strategy/cartesian/side_non_robust.hpp
View File

@@ -1,6 +1,6 @@
// Boost.Geometry
// Copyright (c) 2020, Oracle and/or its affiliates.
// Copyright (c) 2020-2021, Oracle and/or its affiliates.
// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
@@ -14,6 +14,8 @@
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/util/precise_math.hpp>
#include <boost/geometry/arithmetic/determinant.hpp>
namespace boost { namespace geometry
{
@@ -21,7 +23,7 @@ namespace strategy { namespace side
{
/*!
\brief Adaptive precision predicate to check at which side of a segment a point lies:
\brief Predicate to check at which side of a segment a point lies:
left of segment (>0), right of segment (< 0), on segment (0).
\ingroup strategies
\tparam CalculationType \tparam_calculation
@@ -35,8 +37,6 @@ struct side_non_robust
{
public:
//! \brief Computes double the signed area of the CCW triangle p1, p2, p
#ifndef DOXYGEN_SHOULD_SKIP_THIS
template
<
typename P1,
@@ -51,21 +51,44 @@ public:
P1,
P2,
P
>::type coordinate_type;
>::type CoordinateType;
typedef typename select_most_precise
<
coordinate_type,
CoordinateType,
double
>::type promoted_type;
>::type PromotedType;
auto detleft = (promoted_type(get<0>(p1)) - promoted_type(get<0>(p)))
* (promoted_type(get<1>(p2)) - promoted_type(get<1>(p)));
auto detright = (promoted_type(get<1>(p1)) - promoted_type(get<1>(p)))
* (promoted_type(get<0>(p2)) - promoted_type(get<0>(p)));
return detleft > detright ? 1 : (detleft < detright ? -1 : 0 );
CoordinateType const x = get<0>(p);
CoordinateType const y = get<1>(p);
CoordinateType const sx1 = get<0>(p1);
CoordinateType const sy1 = get<1>(p1);
CoordinateType const sx2 = get<0>(p2);
CoordinateType const sy2 = get<1>(p2);
//non-robust 1
//the following is 2x slower in some generic cases when compiled with g++
//(tested versions 9 and 10)
//
//auto detleft = (sx1 - x) * (sy2 - y);
//auto detright = (sy1 - y) * (sx2 - x);
//return detleft > detright ? 1 : (detleft < detright ? -1 : 0 );
//non-robust 2
PromotedType const dx = sx2 - sx1;
PromotedType const dy = sy2 - sy1;
PromotedType const dpx = x - sx1;
PromotedType const dpy = y - sy1;
PromotedType sv = geometry::detail::determinant<PromotedType>
(
dx, dy,
dpx, dpy
);
PromotedType const zero = PromotedType();
return sv == 0 ? 0 : sv > zero ? 1 : -1;
}
#endif
};

119
include/boost/geometry/strategy/cartesian/side_robust.hpp
View File

@@ -7,6 +7,7 @@
// This file was modified by Oracle on 2021.
// Modifications copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -17,9 +18,14 @@
#ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_ROBUST_HPP
#define BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_ROBUST_HPP
#include <boost/geometry/strategy/cartesian/side_non_robust.hpp>
#include <boost/geometry/strategies/side.hpp>
#include <boost/geometry/util/select_most_precise.hpp>
#include <boost/geometry/util/select_calculation_type.hpp>
#include <boost/geometry/util/precise_math.hpp>
#include <boost/geometry/util/math.hpp>
namespace boost { namespace geometry
{
@@ -27,6 +33,27 @@ namespace boost { namespace geometry
namespace strategy { namespace side
{
struct epsilon_equals_policy
{
public:
template <typename Policy, typename T1, typename T2>
static bool apply(T1 const& a, T2 const& b, Policy const& policy)
{
return boost::geometry::math::detail::equals_by_policy(a, b, policy);
}
};
struct fp_equals_policy
{
public:
template <typename Policy, typename T1, typename T2>
static bool apply(T1 const& a, T2 const& b, Policy const&)
{
return a == b;
}
};
/*!
\brief Adaptive precision predicate to check at which side of a segment a point lies:
left of segment (>0), right of segment (< 0), on segment (0).
@@ -38,21 +65,52 @@ namespace strategy { namespace side
template
<
typename CalculationType = void,
typename EqualsPolicy = epsilon_equals_policy,
std::size_t Robustness = 3
>
struct side_robust
{
template <typename CT>
struct epsilon_policy
{
using Policy = boost::geometry::math::detail::equals_factor_policy<CT>;
epsilon_policy() {}
template <typename Type>
epsilon_policy(Type const& a, Type const& b, Type const& c, Type const& d)
: m_policy(a, b, c, d)
{}
Policy m_policy;
public:
template <typename T1, typename T2>
bool apply(T1 a, T2 b) const
{
return EqualsPolicy::apply(a, b, m_policy);
}
};
public:
//! \brief Computes double the signed area of the CCW triangle p1, p2, p
typedef cartesian_tag cs_tag;
//! \brief Computes the sign of the CCW triangle p1, p2, p
template
<
typename PromotedType,
typename P1,
typename P2,
typename P
typename P,
typename EpsPolicyInternal,
std::enable_if_t<std::is_fundamental<PromotedType>::value, int> = 0
>
static inline PromotedType side_value(P1 const& p1, P2 const& p2,
P const& p)
static inline PromotedType side_value(P1 const& p1,
P2 const& p2,
P const& p,
EpsPolicyInternal& eps_policy)
{
using vec2d = ::boost::geometry::detail::precise_math::vec2d<PromotedType>;
vec2d pa;
@@ -65,7 +123,22 @@ public:
pc.x = get<0>(p);
pc.y = get<1>(p);
return ::boost::geometry::detail::precise_math::orient2d
<PromotedType, Robustness>(pa, pb, pc);
<PromotedType, Robustness>(pa, pb, pc, eps_policy);
}
template
<
typename PromotedType,
typename P1,
typename P2,
typename P,
typename EpsPolicyInternal,
std::enable_if_t<!std::is_fundamental<PromotedType>::value, int> = 0
>
static inline auto side_value(P1 const& p1, P2 const& p2, P const& p,
EpsPolicyInternal&)
{
return side_non_robust<>::apply(p1, p2, p);
}
#ifndef DOXYGEN_SHOULD_SKIP_THIS
@@ -77,28 +150,48 @@ public:
>
static inline int apply(P1 const& p1, P2 const& p2, P const& p)
{
typedef typename select_calculation_type_alt
using coordinate_type = typename select_calculation_type_alt
<
CalculationType,
P1,
P2,
P
>::type coordinate_type;
typedef typename select_most_precise
>::type;
using promoted_type = typename select_most_precise
<
coordinate_type,
double
>::type promoted_type;
>::type;
promoted_type sv = side_value<promoted_type>(p1, p2, p);
return sv > 0 ? 1
: sv < 0 ? -1
: 0;
epsilon_policy<promoted_type> epsp;
promoted_type sv = side_value<promoted_type>(p1, p2, p, epsp);
promoted_type const zero = promoted_type();
return epsp.apply(sv, zero) ? 0
: sv > zero ? 1
: -1;
}
#endif
};
#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
namespace services
{
template <typename CalculationType>
struct default_strategy<cartesian_tag, CalculationType>
{
typedef side_robust<CalculationType> type;
};
}
#endif
}} // namespace strategy::side
}} // namespace boost::geometry

41
include/boost/geometry/util/precise_math.hpp
View File

@@ -5,6 +5,10 @@
// Contributed and/or modified by Tinko Bartels,
// as part of Google Summer of Code 2019 program.
// This file was modified by Oracle on 2021.
// Modifications copyright (c) 2021, Oracle and/or its affiliates.
// Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle
// Use, modification and distribution is subject to the Boost Software License,
// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
// http://www.boost.org/LICENSE_1_0.txt)
@@ -157,7 +161,9 @@ inline int fast_expansion_sum_zeroelim(
int n = InSize2)
{
std::array<RealNumber, 2> Qh;
int i_e = 0, i_f = 0, i_h = 0;
int i_e = 0;
int i_f = 0;
int i_h = 0;
if (std::abs(f[0]) > std::abs(e[0]))
{
Qh[0] = e[i_e++];
@@ -282,14 +288,16 @@ inline RealNumber orient2dtail(vec2d<RealNumber> const& p1,
std::array<RealNumber, 2>& t4,
std::array<RealNumber, 2>& t5_01,
std::array<RealNumber, 2>& t6_01,
RealNumber const& magnitude
)
RealNumber const& magnitude)
{
t5_01[1] = two_product_tail(t1[0], t2[0], t5_01[0]);
t6_01[1] = two_product_tail(t3[0], t4[0], t6_01[0]);
std::array<RealNumber, 4> tA_03 = two_two_expansion_diff(t5_01, t6_01);
RealNumber det = std::accumulate(tA_03.begin(), tA_03.end(), static_cast<RealNumber>(0));
if(Robustness == 1) return det;
if (Robustness == 1)
{
return det;
}
// see p.39, mind the different definition of epsilon for error bound
RealNumber B_relative_bound =
(1 + 3 * std::numeric_limits<RealNumber>::epsilon())
@@ -347,29 +355,37 @@ inline RealNumber orient2dtail(vec2d<RealNumber> const& p1,
return(D[D_nz - 1]);
}
// see page 38, Figure 21 for the calculations, notation follows the notation in the figure.
// see page 38, Figure 21 for the calculations, notation follows the notation
// in the figure.
template
<
typename RealNumber,
std::size_t Robustness = 3
std::size_t Robustness = 3,
typename EpsPolicy
>
inline RealNumber orient2d(vec2d<RealNumber> const& p1,
vec2d<RealNumber> const& p2,
vec2d<RealNumber> const& p3)
vec2d<RealNumber> const& p3,
EpsPolicy& eps_policy)
{
if(Robustness == 0)
{
return (p1.x - p3.x) * (p2.y - p3.y) - (p1.y - p3.y) * (p2.x - p3.x);
}
std::array<RealNumber, 2> t1, t2, t3, t4;
t1[0] = p1.x - p3.x;
t2[0] = p2.y - p3.y;
t3[0] = p1.y - p3.y;
t4[0] = p2.x - p3.x;
eps_policy = EpsPolicy(t1[0], t2[0], t3[0], t4[0]);
std::array<RealNumber, 2> t5_01, t6_01;
t5_01[0] = t1[0] * t2[0];
t6_01[0] = t3[0] * t4[0];
RealNumber det = t5_01[0] - t6_01[0];
if (Robustness == 0)
{
return det;
}
RealNumber const magnitude = std::abs(t5_01[0]) + std::abs(t6_01[0]);
// see p.39, mind the different definition of epsilon for error bound
@@ -388,7 +404,8 @@ inline RealNumber orient2d(vec2d<RealNumber> const& p1,
//obvious
return det;
}
return orient2dtail<RealNumber, Robustness>(p1, p2, p3, t1, t2, t3, t4, t5_01, t6_01, magnitude);
return orient2dtail<RealNumber, Robustness>(p1, p2, p3, t1, t2, t3, t4,
t5_01, t6_01, magnitude);
}
// This method adaptively computes increasingly precise approximations of the following

32
include/boost/geometry/util/promote_floating_point.hpp
View File

@@ -18,35 +18,9 @@
#ifndef BOOST_GEOMETRY_UTIL_PROMOTE_FLOATING_POINT_HPP
#define BOOST_GEOMETRY_UTIL_PROMOTE_FLOATING_POINT_HPP
#include <boost/config/header_deprecated.hpp>
BOOST_HEADER_DEPRECATED("boost/geometry/core/coordinate_promotion.hpp")
#include <type_traits>
namespace boost { namespace geometry
{
/*!
\brief Meta-function converting, if necessary, to "a floating point" type
\details
- if input type is integer, type is double
- else type is input type
\ingroup utility
*/
template <typename T, typename PromoteIntegerTo = double>
struct promote_floating_point
{
typedef std::conditional_t
<
std::is_integral<T>::value,
PromoteIntegerTo,
T
> type;
};
}} // namespace boost::geometry
#include <boost/geometry/core/coordinate_promotion.hpp>
#endif // BOOST_GEOMETRY_UTIL_PROMOTE_FLOATING_POINT_HPP