strategy/cartesian/side_robust.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2019 Tinko Bartels, Berlin, Germany.
0004
0005 // Contributed and/or modified by Tinko Bartels,
0006 //   as part of Google Summer of Code 2019 program.
0007
0008 // This file was modified by Oracle on 2021.
0009 // Modifications copyright (c) 2021, Oracle and/or its affiliates.
0010
0011 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0012 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0013
0014 // Use, modification and distribution is subject to the Boost Software License,
0015 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0016 // http://www.boost.org/LICENSE_1_0.txt)
0017
0018 #ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_ROBUST_HPP
0019 #define BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_ROBUST_HPP
0020
0021 #include <boost/geometry/core/config.hpp>
0022 #include <boost/geometry/strategy/cartesian/side_non_robust.hpp>
0023
0024 #include <boost/geometry/strategies/side.hpp>
0025
0026 #include <boost/geometry/util/select_most_precise.hpp>
0027 #include <boost/geometry/util/select_calculation_type.hpp>
0028 #include <boost/geometry/util/precise_math.hpp>
0029 #include <boost/geometry/util/math.hpp>
0030
0031 namespace boost { namespace geometry
0032 {
0033
0034 namespace strategy { namespace side
0035 {
0036
0037 struct epsilon_equals_policy
0038 {
0039 public:
0040     template <typename Policy, typename T1, typename T2>
0041     static bool apply(T1 const& a, T2 const& b, Policy const& policy)
0042     {
0043         return boost::geometry::math::detail::equals_by_policy(a, b, policy);
0044     }
0045 };
0046
0047 struct fp_equals_policy
0048 {
0049 public:
0050     template <typename Policy, typename T1, typename T2>
0051     static bool apply(T1 const& a, T2 const& b, Policy const&)
0052     {
0053         return a == b;
0054     }
0055 };
0056
0057
0058 /*!
0059 \brief Adaptive precision predicate to check at which side of a segment a point lies:
0060     left of segment (>0), right of segment (< 0), on segment (0).
0061 \ingroup strategies
0062 \tparam CalculationType \tparam_calculation (numeric_limits<ct>::epsilon() and numeric_limits<ct>::digits must be supported for calculation type ct)
0063 \tparam Robustness std::size_t value from 0 (fastest) to 3 (default, guarantees correct results).
0064 \details This predicate determines at which side of a segment a point lies using an algorithm that is adapted from orient2d as described in "Adaptive Precision Floating-Point Arithmetic and Fast Robust Geometric Predicates" by Jonathan Richard Shewchuk ( https://dl.acm.org/citation.cfm?doid=237218.237337 ). More information and copies of the paper can also be found at https://www.cs.cmu.edu/~quake/robust.html . It is designed to be adaptive in the sense that it should be fast for inputs that lead to correct results with plain float operations but robust for inputs that require higher precision arithmetics.
0065  */
0066 template
0067 <
0068     typename CalculationType = void,
0069     typename EqualsPolicy = epsilon_equals_policy,
0070     std::size_t Robustness = 3
0071 >
0072 struct side_robust
0073 {
0074
0075     template <typename CT>
0076     struct epsilon_policy
0077     {
0078         using Policy = boost::geometry::math::detail::equals_factor_policy<CT>;
0079
0080         epsilon_policy() {}
0081
0082         template <typename Type>
0083         epsilon_policy(Type const& a, Type const& b, Type const& c, Type const& d)
0084             : m_policy(a, b, c, d)
0085         {}
0086         Policy m_policy;
0087
0088     public:
0089
0090         template <typename T1, typename T2>
0091         bool apply(T1 a, T2 b) const
0092         {
0093             return EqualsPolicy::apply(a, b, m_policy);
0094         }
0095     };
0096
0097 public:
0098
0099     typedef cartesian_tag cs_tag;
0100
0101     //! \brief Computes the sign of the CCW triangle p1, p2, p
0102     template
0103     <
0104         typename PromotedType,
0105         typename P1,
0106         typename P2,
0107         typename P,
0108         typename EpsPolicyInternal,
0109         std::enable_if_t<std::is_fundamental<PromotedType>::value, int> = 0
0110     >
0111     static inline PromotedType side_value(P1 const& p1,
0112                                           P2 const& p2,
0113                                           P const& p,
0114                                           EpsPolicyInternal& eps_policy)
0115     {
0116         using vec2d = ::boost::geometry::detail::precise_math::vec2d<PromotedType>;
0117         vec2d pa;
0118         pa.x = get<0>(p1);
0119         pa.y = get<1>(p1);
0120         vec2d pb;
0121         pb.x = get<0>(p2);
0122         pb.y = get<1>(p2);
0123         vec2d pc;
0124         pc.x = get<0>(p);
0125         pc.y = get<1>(p);
0126         return ::boost::geometry::detail::precise_math::orient2d
0127             <PromotedType, Robustness>(pa, pb, pc, eps_policy);
0128     }
0129
0130     template
0131     <
0132         typename PromotedType,
0133         typename P1,
0134         typename P2,
0135         typename P,
0136         typename EpsPolicyInternal,
0137         std::enable_if_t<!std::is_fundamental<PromotedType>::value, int> = 0
0138     >
0139     static inline auto side_value(P1 const& p1, P2 const& p2, P const& p,
0140                                   EpsPolicyInternal&)
0141     {
0142         return side_non_robust<>::apply(p1, p2, p);
0143     }
0144
0145 #ifndef DOXYGEN_SHOULD_SKIP_THIS
0146     template
0147     <
0148         typename P1,
0149         typename P2,
0150         typename P
0151     >
0152     static inline int apply(P1 const& p1, P2 const& p2, P const& p)
0153     {
0154         using coordinate_type = typename select_calculation_type_alt
0155             <
0156                 CalculationType,
0157                 P1,
0158                 P2,
0159                 P
0160             >::type;
0161
0162         using promoted_type = typename select_most_precise
0163             <
0164                 coordinate_type,
0165                 double
0166             >::type;
0167
0168         epsilon_policy<promoted_type> epsp;
0169         promoted_type sv = side_value<promoted_type>(p1, p2, p, epsp);
0170         promoted_type const zero = promoted_type();
0171
0172         return epsp.apply(sv, zero) ? 0
0173             : sv > zero ? 1
0174             : -1;
0175     }
0176
0177 #endif
0178
0179 };
0180
0181 }} // namespace strategy::side
0182
0183 }} // namespace boost::geometry
0184
0185 #endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_ROBUST_HPP