side_by_triangle.hpp 8.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258
  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/select_calculation_type.hpp>
  25. #include <boost/geometry/util/select_most_precise.hpp>
  26. namespace boost { namespace geometry
  27. {
  28. namespace strategy { namespace side
  29. {
  30. /*!
  31. \brief Check at which side of a segment a point lies:
  32. left of segment (> 0), right of segment (< 0), on segment (0)
  33. \ingroup strategies
  34. \tparam CalculationType \tparam_calculation
  35. */
  36. template <typename CalculationType = void>
  37. class side_by_triangle
  38. {
  39. template <typename Policy>
  40. struct eps_policy
  41. {
  42. eps_policy() {}
  43. template <typename Type>
  44. eps_policy(Type const& a, Type const& b, Type const& c, Type const& d)
  45. : policy(a, b, c, d)
  46. {}
  47. Policy policy;
  48. };
  49. struct eps_empty
  50. {
  51. eps_empty() {}
  52. template <typename Type>
  53. eps_empty(Type const&, Type const&, Type const&, Type const&) {}
  54. };
  55. public :
  56. using cs_tag = cartesian_tag;
  57. // Template member function, because it is not always trivial
  58. // or convenient to explicitly mention the typenames in the
  59. // strategy-struct itself.
  60. // Types can be all three different. Therefore it is
  61. // not implemented (anymore) as "segment"
  62. template
  63. <
  64. typename CoordinateType,
  65. typename PromotedType,
  66. typename P1,
  67. typename P2,
  68. typename P,
  69. typename EpsPolicy
  70. >
  71. static inline
  72. PromotedType side_value(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & eps_policy)
  73. {
  74. CoordinateType const x = get<0>(p);
  75. CoordinateType const y = get<1>(p);
  76. CoordinateType const sx1 = get<0>(p1);
  77. CoordinateType const sy1 = get<1>(p1);
  78. CoordinateType const sx2 = get<0>(p2);
  79. CoordinateType const sy2 = get<1>(p2);
  80. PromotedType const dx = sx2 - sx1;
  81. PromotedType const dy = sy2 - sy1;
  82. PromotedType const dpx = x - sx1;
  83. PromotedType const dpy = y - sy1;
  84. eps_policy = EpsPolicy(dx, dy, dpx, dpy);
  85. return geometry::detail::determinant<PromotedType>
  86. (
  87. dx, dy,
  88. dpx, dpy
  89. );
  90. }
  91. template
  92. <
  93. typename CoordinateType,
  94. typename PromotedType,
  95. typename P1,
  96. typename P2,
  97. typename P
  98. >
  99. static inline
  100. PromotedType side_value(P1 const& p1, P2 const& p2, P const& p)
  101. {
  102. eps_empty dummy;
  103. return side_value<CoordinateType, PromotedType>(p1, p2, p, dummy);
  104. }
  105. template
  106. <
  107. typename CoordinateType,
  108. typename PromotedType,
  109. bool AreAllIntegralCoordinates
  110. >
  111. struct compute_side_value
  112. {
  113. template <typename P1, typename P2, typename P, typename EpsPolicy>
  114. static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
  115. {
  116. return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
  117. }
  118. };
  119. template <typename CoordinateType, typename PromotedType>
  120. struct compute_side_value<CoordinateType, PromotedType, false>
  121. {
  122. template <typename P1, typename P2, typename P, typename EpsPolicy>
  123. static inline PromotedType apply(P1 const& p1, P2 const& p2, P const& p, EpsPolicy & epsp)
  124. {
  125. // For robustness purposes, first check if any two points are
  126. // the same; in this case simply return that the points are
  127. // collinear
  128. if (equals_point_point(p1, p2)
  129. || equals_point_point(p1, p)
  130. || equals_point_point(p2, p))
  131. {
  132. return PromotedType(0);
  133. }
  134. // The side_by_triangle strategy computes the signed area of
  135. // the point triplet (p1, p2, p); as such it is (in theory)
  136. // invariant under cyclic permutations of its three arguments.
  137. //
  138. // In the context of numerical errors that arise in
  139. // floating-point computations, and in order to make the strategy
  140. // consistent with respect to cyclic permutations of its three
  141. // arguments, we cyclically permute them so that the first
  142. // argument is always the lexicographically smallest point.
  143. using less = compare::cartesian<compare::less, compare::equals_epsilon>;
  144. if (less::apply(p, p1))
  145. {
  146. if (less::apply(p, p2))
  147. {
  148. // p is the lexicographically smallest
  149. return side_value<CoordinateType, PromotedType>(p, p1, p2, epsp);
  150. }
  151. else
  152. {
  153. // p2 is the lexicographically smallest
  154. return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
  155. }
  156. }
  157. if (less::apply(p1, p2))
  158. {
  159. // p1 is the lexicographically smallest
  160. return side_value<CoordinateType, PromotedType>(p1, p2, p, epsp);
  161. }
  162. else
  163. {
  164. // p2 is the lexicographically smallest
  165. return side_value<CoordinateType, PromotedType>(p2, p, p1, epsp);
  166. }
  167. }
  168. };
  169. template <typename P1, typename P2, typename P>
  170. static inline int apply(P1 const& p1, P2 const& p2, P const& p)
  171. {
  172. using coor_t = typename select_calculation_type_alt<CalculationType, P1, P2, P>::type;
  173. // Promote float->double, small int->int
  174. using promoted_t = typename select_most_precise<coor_t, double>::type;
  175. bool const are_all_integral_coordinates =
  176. std::is_integral<typename coordinate_type<P1>::type>::value
  177. && std::is_integral<typename coordinate_type<P2>::type>::value
  178. && std::is_integral<typename coordinate_type<P>::type>::value;
  179. eps_policy< math::detail::equals_factor_policy<promoted_t> > epsp;
  180. promoted_t s = compute_side_value
  181. <
  182. coor_t, promoted_t, are_all_integral_coordinates
  183. >::apply(p1, p2, p, epsp);
  184. promoted_t const zero = promoted_t();
  185. return math::detail::equals_by_policy(s, zero, epsp.policy) ? 0
  186. : s > zero ? 1
  187. : -1;
  188. }
  189. private:
  190. template <typename P1, typename P2>
  191. static inline bool equals_point_point(P1 const& p1, P2 const& p2)
  192. {
  193. return strategy::within::cartesian_point_point::apply(p1, p2);
  194. }
  195. };
  196. #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
  197. namespace services
  198. {
  199. template <typename CalculationType>
  200. struct default_strategy<cartesian_tag, CalculationType>
  201. {
  202. using type = side_by_triangle<CalculationType>;
  203. };
  204. }
  205. #endif
  206. }} // namespace strategy::side
  207. }} // namespace boost::geometry
  208. #endif // BOOST_GEOMETRY_STRATEGY_CARTESIAN_SIDE_BY_TRIANGLE_HPP