strategy/cartesian/side_by_triangle.hpp

0001 // Boost.Geometry (aka GGL, Generic Geometry Library)
0002
0003 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
0004 // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
0005 // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
0006
0007 // This file was modified by Oracle on 2015-2023.
0008 // Modifications copyright (c) 2015-2023, Oracle and/or its affiliates.
0009
0010 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0011 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
0012 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0013
0014 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
0015 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
0016
0017 // Use, modification and distribution is subject to the Boost Software License,
0018 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0019 // http://www.boost.org/LICENSE_1_0.txt)
0020
0021 #ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP
0022 #define BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP
0023
0024
0025 #include <type_traits>
0026
0027 #include <boost/geometry/core/config.hpp>
0028 #include <boost/geometry/arithmetic/determinant.hpp>
0029
0030 #include <boost/geometry/core/access.hpp>
0031
0032 #include <boost/geometry/strategies/cartesian/point_in_point.hpp>
0033 #include <boost/geometry/strategies/compare.hpp>
0034 #include <boost/geometry/strategies/side.hpp>
0035
0036 #include <boost/geometry/util/promote_integral.hpp>
0037 #include <boost/geometry/util/select_calculation_type.hpp>
0038 #include <boost/geometry/util/select_most_precise.hpp>
0039
0040
0041 namespace boost { namespace geometry
0042 {
0043
0044 namespace strategy { namespace side
0045 {
0046
0047 /*!
0048 \brief Check at which side of a segment a point lies:
0049     left of segment (> 0), right of segment (< 0), on segment (0)
0050 \ingroup strategies
0051 \tparam CalculationType \tparam_calculation
0052  */
0053 template <typename CalculationType = void>
0054 class side_by_triangle
0055 {
0056     template <typename Policy>
0057     struct eps_policy
0058     {
0059         eps_policy() {}
0060         template <typename Type>
0061         eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
0062             : policy(a, b, c, d)
0063         {}
0064         Policy policy;
0065     };
0066
0067     struct eps_empty
0068     {
0069         eps_empty() {}
0070         template <typename Type>
0071         eps_empty(Type const&, Type const&, Type const&, Type const&) {}
0072     };
0073
0074 public :
0075     using cs_tag = cartesian_tag;
0076
0077     // Template member function, because it is not always trivial
0078     // or convenient to explicitly mention the typenames in the
0079     // strategy-struct itself.
0080
0081     // Types can be all three different. Therefore it is
0082     // not implemented (anymore) as "segment"
0083
0084     template
0085     <
0086         typename CoordinateType,
0087         typename PromotedType,
0088         typename P1,
0089         typename P2,
0090         typename P,
0091         typename EpsPolicy
0092     >
0093     static inline
0094     PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
0095     {
0096         CoordinateType const x = get<0>(p);
0097         CoordinateType const y = get<1>(p);
0098
0099         CoordinateType const sx1 = get<0>(p1);
0100         CoordinateType const sy1 = get<1>(p1);
0101         CoordinateType const sx2 = get<0>(p2);
0102         CoordinateType const sy2 = get<1>(p2);
0103
0104         PromotedType const dx = sx2 - sx1;
0105         PromotedType const dy = sy2 - sy1;
0106         PromotedType const dpx = x - sx1;
0107         PromotedType const dpy = y - sy1;
0108
0109         eps_policy = EpsPolicy(dx, dy, dpx, dpy);
0110
0111         return geometry::detail::determinant<PromotedType>
0112                 (
0113                     dx, dy,
0114                     dpx, dpy
0115                 );
0116
0117     }
0118
0119     template
0120     <
0121         typename CoordinateType,
0122         typename PromotedType,
0123         typename P1,
0124         typename P2,
0125         typename P
0126     >
0127     static inline
0128     PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
0129     {
0130         eps_empty dummy;
0131         return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
0132     }
0133
0134
0135     template
0136     <
0137         typename CoordinateType,
0138         typename PromotedType,
0139         bool AreAllIntegralCoordinates
0140     >
0141     struct compute_side_value
0142     {
0143         template <typename P1, typename P2, typename P, typename EpsPolicy>
0144         static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
0145         {
0146             return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
0147         }
0148     };
0149
0150     template <typename CoordinateType, typename PromotedType>
0151     struct compute_side_value<CoordinateType, PromotedType, false>
0152     {
0153         template <typename P1, typename P2, typename P, typename EpsPolicy>
0154         static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
0155         {
0156             // For robustness purposes, first check if any two points are
0157             // the same; in this case simply return that the points are
0158             // collinear
0159             if (equals_point_point(p1, p2)
0160                 || equals_point_point(p1, p)
0161                 || equals_point_point(p2, p))
0162             {
0163                 return PromotedType(0);
0164             }
0165
0166             // The side_by_triangle strategy computes the signed area of
0167             // the point triplet (p1, p2, p); as such it is (in theory)
0168             // invariant under cyclic permutations of its three arguments.
0169             //
0170             // In the context of numerical errors that arise in
0171             // floating-point computations, and in order to make the strategy
0172             // consistent with respect to cyclic permutations of its three
0173             // arguments, we cyclically permute them so that the first
0174             // argument is always the lexicographically smallest point.
0175
0176             using less = compare::cartesian<compare::less, compare::equals_epsilon>;
0177
0178             if (less::apply(p, p1))
0179             {
0180                 if (less::apply(p, p2))
0181                 {
0182                     // p is the lexicographically smallest
0183                     return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
0184                 }
0185                 else
0186                 {
0187                     // p2 is the lexicographically smallest
0188                     return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
0189                 }
0190             }
0191
0192             if (less::apply(p1, p2))
0193             {
0194                 // p1 is the lexicographically smallest
0195                 return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
0196             }
0197             else
0198             {
0199                 // p2 is the lexicographically smallest
0200                 return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
0201             }
0202         }
0203     };
0204
0205
0206     template <typename P1, typename P2, typename P>
0207     static inline int apply(P1 const& p1, P2 const& p2, P const& p)
0208     {
0209         constexpr bool are_all_integral_coordinates =
0210             std::is_integral<coordinate_type_t<P1>>::value
0211             && std::is_integral<coordinate_type_t<P2>>::value
0212             && std::is_integral<coordinate_type_t<P>>::value;
0213
0214         // Promote float to double
0215         // For integer: short -> int -> long
0216         // For larger integers: long, long long, std::int64_t all stay as they are (on a Mac)
0217         using coor_t = typename select_calculation_type_alt<CalculationType, P1, P2, P>::type;
0218         using promoted_t = std::conditional_t
0219             <
0220                 are_all_integral_coordinates,
0221                 typename promote_integral<coor_t>::type,
0222                 typename select_most_precise<coor_t, double>::type
0223             >;
0224
0225         eps_policy< math::detail::equals_factor_policy<promoted_t> > epsp;
0226         promoted_t const s = compute_side_value
0227             <
0228                 coor_t, promoted_t, are_all_integral_coordinates
0229             >::apply(p1, p2, p, epsp);
0230
0231         static promoted_t const zero = promoted_t();
0232         return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
0233             : s > zero ? 1
0234             : -1;
0235     }
0236
0237 private:
0238     template <typename P1, typename P2>
0239     static inline bool equals_point_point(P1 const& p1, P2 const& p2)
0240     {
0241         return strategy::within::cartesian_point_point::apply(p1, p2);
0242     }
0243 };
0244
0245 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0246
0247 namespace services
0248 {
0249
0250 template <typename CalculationType>
0251 struct default_strategy<cartesian_tag, CalculationType>
0252 {
0253     using type = side_by_triangle<CalculationType>;
0254 };
0255
0256 }
0257
0258 #endif
0259
0260 }} // namespace strategy::side
0261
0262 }} // namespace boost::geometry
0263
0264
0265 #endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP