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/select_calculation_type.hpp>
0037 #include <boost/geometry/util/select_most_precise.hpp>
0038
0039
0040 namespace boost { namespace geometry
0041 {
0042
0043 namespace strategy { namespace side
0044 {
0045
0046 /*!
0047 \brief Check at which side of a segment a point lies:
0048     left of segment (> 0), right of segment (< 0), on segment (0)
0049 \ingroup strategies
0050 \tparam CalculationType \tparam_calculation
0051  */
0052 template <typename CalculationType = void>
0053 class side_by_triangle
0054 {
0055     template <typename Policy>
0056     struct eps_policy
0057     {
0058         eps_policy() {}
0059         template <typename Type>
0060         eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
0061             : policy(a, b, c, d)
0062         {}
0063         Policy policy;
0064     };
0065
0066     struct eps_empty
0067     {
0068         eps_empty() {}
0069         template <typename Type>
0070         eps_empty(Type const&, Type const&, Type const&, Type const&) {}
0071     };
0072
0073 public :
0074     using cs_tag = cartesian_tag;
0075
0076     // Template member function, because it is not always trivial
0077     // or convenient to explicitly mention the typenames in the
0078     // strategy-struct itself.
0079
0080     // Types can be all three different. Therefore it is
0081     // not implemented (anymore) as "segment"
0082
0083     template
0084     <
0085         typename CoordinateType,
0086         typename PromotedType,
0087         typename P1,
0088         typename P2,
0089         typename P,
0090         typename EpsPolicy
0091     >
0092     static inline
0093     PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
0094     {
0095         CoordinateType const x = get<0>(p);
0096         CoordinateType const y = get<1>(p);
0097
0098         CoordinateType const sx1 = get<0>(p1);
0099         CoordinateType const sy1 = get<1>(p1);
0100         CoordinateType const sx2 = get<0>(p2);
0101         CoordinateType const sy2 = get<1>(p2);
0102
0103         PromotedType const dx = sx2 - sx1;
0104         PromotedType const dy = sy2 - sy1;
0105         PromotedType const dpx = x - sx1;
0106         PromotedType const dpy = y - sy1;
0107
0108         eps_policy = EpsPolicy(dx, dy, dpx, dpy);
0109
0110         return geometry::detail::determinant<PromotedType>
0111                 (
0112                     dx, dy,
0113                     dpx, dpy
0114                 );
0115
0116     }
0117
0118     template
0119     <
0120         typename CoordinateType,
0121         typename PromotedType,
0122         typename P1,
0123         typename P2,
0124         typename P
0125     >
0126     static inline
0127     PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
0128     {
0129         eps_empty dummy;
0130         return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
0131     }
0132
0133
0134     template
0135     <
0136         typename CoordinateType,
0137         typename PromotedType,
0138         bool AreAllIntegralCoordinates
0139     >
0140     struct compute_side_value
0141     {
0142         template <typename P1, typename P2, typename P, typename EpsPolicy>
0143         static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
0144         {
0145             return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
0146         }
0147     };
0148
0149     template <typename CoordinateType, typename PromotedType>
0150     struct compute_side_value<CoordinateType, PromotedType, false>
0151     {
0152         template <typename P1, typename P2, typename P, typename EpsPolicy>
0153         static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
0154         {
0155             // For robustness purposes, first check if any two points are
0156             // the same; in this case simply return that the points are
0157             // collinear
0158             if (equals_point_point(p1, p2)
0159                 || equals_point_point(p1, p)
0160                 || equals_point_point(p2, p))
0161             {
0162                 return PromotedType(0);
0163             }
0164
0165             // The side_by_triangle strategy computes the signed area of
0166             // the point triplet (p1, p2, p); as such it is (in theory)
0167             // invariant under cyclic permutations of its three arguments.
0168             //
0169             // In the context of numerical errors that arise in
0170             // floating-point computations, and in order to make the strategy
0171             // consistent with respect to cyclic permutations of its three
0172             // arguments, we cyclically permute them so that the first
0173             // argument is always the lexicographically smallest point.
0174
0175             using less = compare::cartesian<compare::less, compare::equals_epsilon>;
0176
0177             if (less::apply(p, p1))
0178             {
0179                 if (less::apply(p, p2))
0180                 {
0181                     // p is the lexicographically smallest
0182                     return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
0183                 }
0184                 else
0185                 {
0186                     // p2 is the lexicographically smallest
0187                     return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
0188                 }
0189             }
0190
0191             if (less::apply(p1, p2))
0192             {
0193                 // p1 is the lexicographically smallest
0194                 return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
0195             }
0196             else
0197             {
0198                 // p2 is the lexicographically smallest
0199                 return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
0200             }
0201         }
0202     };
0203
0204
0205     template <typename P1, typename P2, typename P>
0206     static inline int apply(P1 const& p1, P2 const& p2, P const& p)
0207     {
0208         using coor_t = typename select_calculation_type_alt<CalculationType, P1, P2, P>::type;
0209
0210         // Promote float->double, small int->int
0211         using promoted_t = typename select_most_precise<coor_t, double>::type;
0212
0213         bool const are_all_integral_coordinates =
0214             std::is_integral<typename coordinate_type<P1>::type>::value
0215             && std::is_integral<typename coordinate_type<P2>::type>::value
0216             && std::is_integral<typename coordinate_type<P>::type>::value;
0217
0218         eps_policy< math::detail::equals_factor_policy<promoted_t> > epsp;
0219         promoted_t s = compute_side_value
0220             <
0221                 coor_t, promoted_t, are_all_integral_coordinates
0222             >::apply(p1, p2, p, epsp);
0223
0224         promoted_t const zero = promoted_t();
0225         return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
0226             : s > zero ? 1
0227             : -1;
0228     }
0229
0230 private:
0231     template <typename P1, typename P2>
0232     static inline bool equals_point_point(P1 const& p1, P2 const& p2)
0233     {
0234         return strategy::within::cartesian_point_point::apply(p1, p2);
0235     }
0236 };
0237
0238 #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
0239
0240 namespace services
0241 {
0242
0243 template <typename CalculationType>
0244 struct default_strategy<cartesian_tag, CalculationType>
0245 {
0246     using type = side_by_triangle<CalculationType>;
0247 };
0248
0249 }
0250
0251 #endif
0252
0253 }} // namespace strategy::side
0254
0255 }} // namespace boost::geometry
0256
0257
0258 #endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP