closest_points_pt_seg.hpp 4.0 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146
  1. // Boost.Geometry
  2. // Copyright (c) 2021-2023, Oracle and/or its affiliates.
  3. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
  4. // Licensed under the Boost Software License version 1.0.
  5. // http://www.boost.org/users/license.html
  6. #ifndef BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CLOSEST_POINTS_PT_SEG_HPP
  7. #define BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CLOSEST_POINTS_PT_SEG_HPP
  8. #include <boost/geometry/algorithms/convert.hpp>
  9. #include <boost/geometry/core/coordinate_promotion.hpp>
  10. #include <boost/geometry/geometries/point.hpp>
  11. #include <boost/geometry/strategies/cartesian/distance_pythagoras.hpp>
  12. #include <boost/geometry/strategies/closest_points/services.hpp>
  13. namespace boost { namespace geometry
  14. {
  15. namespace strategy { namespace closest_points
  16. {
  17. #ifndef DOXYGEN_NO_DETAIL
  18. namespace detail
  19. {
  20. template <typename CalculationType>
  21. struct compute_closest_point_to_segment
  22. {
  23. template <typename Point, typename PointOfSegment>
  24. static inline auto
  25. apply(Point const& p, PointOfSegment const& p1, PointOfSegment const& p2)
  26. {
  27. // A projected point of points in Integer coordinates must be able to be
  28. // represented in FP.
  29. using fp_point_type = model::point
  30. <
  31. CalculationType,
  32. dimension<PointOfSegment>::value,
  33. typename coordinate_system<PointOfSegment>::type
  34. >;
  35. // For convenience
  36. using fp_vector_type = fp_point_type;
  37. /*
  38. Algorithm [p: (px,py), p1: (x1,y1), p2: (x2,y2)]
  39. VECTOR v(x2 - x1, y2 - y1)
  40. VECTOR w(px - x1, py - y1)
  41. c1 = w . v
  42. c2 = v . v
  43. b = c1 / c2
  44. RETURN POINT(x1 + b * vx, y1 + b * vy)
  45. */
  46. // v is multiplied below with a (possibly) FP-value, so should be in FP
  47. // For consistency we define w also in FP
  48. fp_vector_type v, w, projected;
  49. geometry::convert(p2, v);
  50. geometry::convert(p, w);
  51. geometry::convert(p1, projected);
  52. subtract_point(v, projected);
  53. subtract_point(w, projected);
  54. CalculationType const zero = CalculationType();
  55. CalculationType const c1 = dot_product(w, v);
  56. if (c1 <= zero)
  57. {
  58. fp_vector_type fp_p1;
  59. geometry::convert(p1, fp_p1);
  60. return fp_p1;
  61. }
  62. CalculationType const c2 = dot_product(v, v);
  63. if (c2 <= c1)
  64. {
  65. fp_vector_type fp_p2;
  66. geometry::convert(p2, fp_p2);
  67. return fp_p2;
  68. }
  69. // See above, c1 > 0 AND c2 > c1 so: c2 != 0
  70. CalculationType const b = c1 / c2;
  71. multiply_value(v, b);
  72. add_point(projected, v);
  73. return projected;
  74. }
  75. };
  76. }
  77. #endif // DOXYGEN_NO_DETAIL
  78. template
  79. <
  80. typename CalculationType = void
  81. >
  82. class projected_point
  83. {
  84. public:
  85. // The three typedefs below are necessary to calculate distances
  86. // from segments defined in integer coordinates.
  87. // Integer coordinates can still result in FP distances.
  88. // There is a division, which must be represented in FP.
  89. // So promote.
  90. template <typename Point, typename PointOfSegment>
  91. struct calculation_type
  92. : promote_floating_point
  93. <
  94. typename select_most_precise
  95. <
  96. typename coordinate_type<Point>::type,
  97. typename coordinate_type<PointOfSegment>::type,
  98. CalculationType
  99. >::type
  100. >
  101. {};
  102. template <typename Point, typename PointOfSegment>
  103. inline auto
  104. apply(Point const& p, PointOfSegment const& p1, PointOfSegment const& p2) const
  105. {
  106. assert_dimension_equal<Point, PointOfSegment>();
  107. using calculation_type = typename calculation_type<Point, PointOfSegment>::type;
  108. return detail::compute_closest_point_to_segment<calculation_type>::apply(p, p1, p2);
  109. }
  110. };
  111. }} // namespace strategy::closest_points
  112. }} // namespace boost::geometry
  113. #endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_CLOSEST_POINTS_PT_SEG_HPP