side_by_triangle.hpp 8.4 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266
  1. // Boost.Geometry (aka GGL, Generic Geometry Library)
  2. // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
  3. // Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
  4. // Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
  5. // This file was modified by Oracle on 2015-2023.
  6. // Modifications copyright (c) 2015-2023, Oracle and/or its affiliates.
  7. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
  8. // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
  9. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
  10. // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
  11. // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
  12. // Use, modification and distribution is subject to the Boost Software License,
  13. // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
  14. // http://www.boost.org/LICENSE_1_0.txt)
  15. #ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP
  16. #define BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP
  17. #include <type_traits>
  18. #include <boost/geometry/core/config.hpp>
  19. #include <boost/geometry/arithmetic/determinant.hpp>
  20. #include <boost/geometry/core/access.hpp>
  21. #include <boost/geometry/strategies/cartesian/point_in_point.hpp>
  22. #include <boost/geometry/strategies/compare.hpp>
  23. #include <boost/geometry/strategies/side.hpp>
  24. #include <boost/geometry/util/promote_integral.hpp>
  25. #include <boost/geometry/util/select_calculation_type.hpp>
  26. #include <boost/geometry/util/select_most_precise.hpp>
  27. namespace boost { namespace geometry
  28. {
  29. namespace strategy { namespace side
  30. {
  31. /*!
  32. \brief Check at which side of a segment a point lies:
  33. left of segment (> 0), right of segment (< 0), on segment (0)
  34. \ingroup strategies
  35. \tparam CalculationType \tparam_calculation
  36. */
  37. template <typename CalculationType = void>
  38. class side_by_triangle
  39. {
  40. template <typename Policy>
  41. struct eps_policy
  42. {
  43. eps_policy() {}
  44. template <typename Type>
  45. eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
  46. : policy(a, b, c, d)
  47. {}
  48. Policy policy;
  49. };
  50. struct eps_empty
  51. {
  52. eps_empty() {}
  53. template <typename Type>
  54. eps_empty(Type const&, Type const&, Type const&, Type const&) {}
  55. };
  56. public :
  57. using cs_tag = cartesian_tag;
  58. // Template member function, because it is not always trivial
  59. // or convenient to explicitly mention the typenames in the
  60. // strategy-struct itself.
  61. // Types can be all three different. Therefore it is
  62. // not implemented (anymore) as "segment"
  63. template
  64. <
  65. typename CoordinateType,
  66. typename PromotedType,
  67. typename P1,
  68. typename P2,
  69. typename P,
  70. typename EpsPolicy
  71. >
  72. static inline
  73. PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
  74. {
  75. CoordinateType const x = get<0>(p);
  76. CoordinateType const y = get<1>(p);
  77. CoordinateType const sx1 = get<0>(p1);
  78. CoordinateType const sy1 = get<1>(p1);
  79. CoordinateType const sx2 = get<0>(p2);
  80. CoordinateType const sy2 = get<1>(p2);
  81. PromotedType const dx = sx2 - sx1;
  82. PromotedType const dy = sy2 - sy1;
  83. PromotedType const dpx = x - sx1;
  84. PromotedType const dpy = y - sy1;
  85. eps_policy = EpsPolicy(dx, dy, dpx, dpy);
  86. return geometry::detail::determinant<PromotedType>
  87. (
  88. dx, dy,
  89. dpx, dpy
  90. );
  91. }
  92. template
  93. <
  94. typename CoordinateType,
  95. typename PromotedType,
  96. typename P1,
  97. typename P2,
  98. typename P
  99. >
  100. static inline
  101. PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
  102. {
  103. eps_empty dummy;
  104. return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
  105. }
  106. template
  107. <
  108. typename CoordinateType,
  109. typename PromotedType,
  110. bool AreAllIntegralCoordinates
  111. >
  112. struct compute_side_value
  113. {
  114. template <typename P1, typename P2, typename P, typename EpsPolicy>
  115. static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
  116. {
  117. return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
  118. }
  119. };
  120. template <typename CoordinateType, typename PromotedType>
  121. struct compute_side_value<CoordinateType, PromotedType, false>
  122. {
  123. template <typename P1, typename P2, typename P, typename EpsPolicy>
  124. static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
  125. {
  126. // For robustness purposes, first check if any two points are
  127. // the same; in this case simply return that the points are
  128. // collinear
  129. if (equals_point_point(p1, p2)
  130. || equals_point_point(p1, p)
  131. || equals_point_point(p2, p))
  132. {
  133. return PromotedType(0);
  134. }
  135. // The side_by_triangle strategy computes the signed area of
  136. // the point triplet (p1, p2, p); as such it is (in theory)
  137. // invariant under cyclic permutations of its three arguments.
  138. //
  139. // In the context of numerical errors that arise in
  140. // floating-point computations, and in order to make the strategy
  141. // consistent with respect to cyclic permutations of its three
  142. // arguments, we cyclically permute them so that the first
  143. // argument is always the lexicographically smallest point.
  144. using less = compare::cartesian<compare::less, compare::equals_epsilon>;
  145. if (less::apply(p, p1))
  146. {
  147. if (less::apply(p, p2))
  148. {
  149. // p is the lexicographically smallest
  150. return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
  151. }
  152. else
  153. {
  154. // p2 is the lexicographically smallest
  155. return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
  156. }
  157. }
  158. if (less::apply(p1, p2))
  159. {
  160. // p1 is the lexicographically smallest
  161. return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
  162. }
  163. else
  164. {
  165. // p2 is the lexicographically smallest
  166. return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
  167. }
  168. }
  169. };
  170. template <typename P1, typename P2, typename P>
  171. static inline int apply(P1 const& p1, P2 const& p2, P const& p)
  172. {
  173. constexpr bool are_all_integral_coordinates =
  174. std::is_integral<coordinate_type_t<P1>>::value
  175. && std::is_integral<coordinate_type_t<P2>>::value
  176. && std::is_integral<coordinate_type_t<P>>::value;
  177. // Promote float to double
  178. // For integer: short -> int -> long
  179. // For larger integers: long, long long, std::int64_t all stay as they are (on a Mac)
  180. using coor_t = typename select_calculation_type_alt<CalculationType, P1, P2, P>::type;
  181. using promoted_t = std::conditional_t
  182. <
  183. are_all_integral_coordinates,
  184. typename promote_integral<coor_t>::type,
  185. typename select_most_precise<coor_t, double>::type
  186. >;
  187. eps_policy< math::detail::equals_factor_policy<promoted_t> > epsp;
  188. promoted_t const s = compute_side_value
  189. <
  190. coor_t, promoted_t, are_all_integral_coordinates
  191. >::apply(p1, p2, p, epsp);
  192. static promoted_t const zero = promoted_t();
  193. return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
  194. : s > zero ? 1
  195. : -1;
  196. }
  197. private:
  198. template <typename P1, typename P2>
  199. static inline bool equals_point_point(P1 const& p1, P2 const& p2)
  200. {
  201. return strategy::within::cartesian_point_point::apply(p1, p2);
  202. }
  203. };
  204. #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  205. #if defined(BOOST_GEOMETRY_DEFAULT_STRATEGY_SIDE_USE_SIDE_BY_TRIANGLE)
  206. namespace services
  207. {
  208. template <typename CalculationType>
  209. struct default_strategy<cartesian_tag, CalculationType>
  210. {
  211. using type = side_by_triangle<CalculationType>;
  212. };
  213. } // namespace services
  214. #endif
  215. #endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  216. }} // namespace strategy::side
  217. }} // namespace boost::geometry
  218. #endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP