strategies/spherical/closest_points_pt_seg.hpp

0001 // Boost.Geometry
0002
0003 // Copyright (c) 2024 Adam Wulkiewicz, Lodz, Poland.
0004
0005 // Copyright (c) 2021, Oracle and/or its affiliates.
0006
0007 // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
0008
0009 // Licensed under the Boost Software License version 1.0.
0010 // http://www.boost.org/users/license.html
0011
0012 #ifndef BOOST_GEOMETRY_STRATEGIES_SPHERICAL_CLOSEST_POINTS_CROSS_TRACK_HPP
0013 #define BOOST_GEOMETRY_STRATEGIES_SPHERICAL_CLOSEST_POINTS_CROSS_TRACK_HPP
0014
0015 #include <algorithm>
0016 #include <type_traits>
0017
0018 #include <boost/config.hpp>
0019 #include <boost/concept_check.hpp>
0020
0021 #include <boost/geometry/algorithms/detail/convert_point_to_point.hpp>
0022
0023 #include <boost/geometry/core/cs.hpp>
0024 #include <boost/geometry/core/access.hpp>
0025 #include <boost/geometry/core/coordinate_promotion.hpp>
0026 #include <boost/geometry/core/radian_access.hpp>
0027 #include <boost/geometry/core/tags.hpp>
0028
0029 #include <boost/geometry/formulas/spherical.hpp>
0030
0031 #include <boost/geometry/strategies/distance.hpp>
0032 #include <boost/geometry/strategies/concepts/distance_concept.hpp>
0033 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
0034 #include <boost/geometry/strategies/spherical/distance_cross_track.hpp>
0035 #include <boost/geometry/strategies/spherical/point_in_point.hpp>
0036 #include <boost/geometry/strategies/spherical/intersection.hpp>
0037
0038 #include <boost/geometry/util/math.hpp>
0039 #include <boost/geometry/util/select_calculation_type.hpp>
0040
0041 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0042 #  include <boost/geometry/io/dsv/write.hpp>
0043 #endif
0044
0045
0046 namespace boost { namespace geometry
0047 {
0048
0049 namespace strategy { namespace closest_points
0050 {
0051
0052 template
0053 <
0054     typename CalculationType = void,
0055     typename Strategy = distance::comparable::haversine<double, CalculationType>
0056 >
0057 class cross_track
0058 {
0059 public:
0060     template <typename Point, typename PointOfSegment>
0061     struct calculation_type
0062         : promote_floating_point
0063           <
0064               typename select_calculation_type
0065                   <
0066                       Point,
0067                       PointOfSegment,
0068                       CalculationType
0069                   >::type
0070           >
0071     {};
0072
0073     using radius_type = typename Strategy::radius_type;
0074
0075     cross_track() = default;
0076
0077     explicit inline cross_track(typename Strategy::radius_type const& r)
0078         : m_strategy(r)
0079     {}
0080
0081     inline cross_track(Strategy const& s)
0082         : m_strategy(s)
0083     {}
0084
0085     template <typename Point, typename PointOfSegment>
0086     inline auto apply(Point const& p,
0087                       PointOfSegment const& sp1,
0088                       PointOfSegment const& sp2) const
0089     {
0090         using CT = typename calculation_type<Point, PointOfSegment>::type;
0091
0092         model::point
0093             <
0094                 CT,
0095                 dimension<PointOfSegment>::value,
0096                 coordinate_system_t<PointOfSegment>
0097             > result;
0098
0099         // http://williams.best.vwh.net/avform.htm#XTE
0100         CT d3 = m_strategy.apply(sp1, sp2);
0101
0102         if (geometry::math::equals(d3, 0.0))
0103         {
0104             // "Degenerate" segment, return either d1 or d2
0105             geometry::detail::conversion::convert_point_to_point(sp1, result);
0106             return result;
0107         }
0108
0109         CT d1 = m_strategy.apply(sp1, p);
0110         CT d2 = m_strategy.apply(sp2, p);
0111
0112         auto d_crs_pair = distance::detail::compute_cross_track_pair<CT>::apply(
0113             p, sp1, sp2);
0114
0115         // d1, d2, d3 are in principle not needed, only the sign matters
0116         CT projection1 = cos(d_crs_pair.first) * d1 / d3;
0117         CT projection2 = cos(d_crs_pair.second) * d2 / d3;
0118
0119 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0120         std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " "
0121                   << crs_AD * geometry::math::r2d<CT>() << std::endl;
0122         std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " "
0123                   << crs_AB * geometry::math::r2d<CT>() << std::endl;
0124         std::cout << "Course " << dsv(sp2) << " to " << dsv(sp1) << " "
0125                   << crs_BA * geometry::math::r2d<CT>() << std::endl;
0126         std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " "
0127                   << crs_BD * geometry::math::r2d<CT>() << std::endl;
0128         std::cout << "Projection AD-AB " << projection1 << " : "
0129                   << d_crs1 * geometry::math::r2d<CT>() << std::endl;
0130         std::cout << "Projection BD-BA " << projection2 << " : "
0131                   << d_crs2 * geometry::math::r2d<CT>() << std::endl;
0132         std::cout << " d1: " << (d1 )
0133                   << " d2: " << (d2 )
0134                   << std::endl;
0135 #endif
0136
0137         if (projection1 > 0.0 && projection2 > 0.0)
0138         {
0139 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0140             CT XTD = radius() * geometry::math::abs( asin( sin( d1 ) * sin( d_crs1 ) ));
0141
0142             std::cout << "Projection ON the segment" << std::endl;
0143             std::cout << "XTD: " << XTD
0144                       << " d1: " << (d1 * radius())
0145                       << " d2: " << (d2 * radius())
0146                       << std::endl;
0147 #endif
0148             auto distance = distance::detail::compute_cross_track_distance::apply(
0149                 d_crs_pair.first, d1);
0150
0151             CT lon1 = geometry::get_as_radian<0>(sp1);
0152             CT lat1 = geometry::get_as_radian<1>(sp1);
0153             CT lon2 = geometry::get_as_radian<0>(sp2);
0154             CT lat2 = geometry::get_as_radian<1>(sp2);
0155
0156             CT dist = CT(2) * asin(math::sqrt(distance)) * m_strategy.radius();
0157             CT dist_d1 = CT(2) * asin(math::sqrt(d1)) * m_strategy.radius();
0158
0159             // Note: this is similar to spherical computation in geographic
0160             // point_segment_distance formula
0161             CT earth_radius = m_strategy.radius();
0162             CT cos_frac = cos(dist_d1 / earth_radius) / cos(dist / earth_radius);
0163             CT s14_sph = cos_frac >= 1
0164                 ? CT(0) : cos_frac <= -1 ? math::pi<CT>() * earth_radius
0165                                          : acos(cos_frac) * earth_radius;
0166
0167             CT a12 = geometry::formula::spherical_azimuth<>(lon1, lat1, lon2, lat2);
0168             auto res_direct = geometry::formula::spherical_direct
0169                 <
0170                     true,
0171                     false
0172                 >(lon1, lat1, s14_sph, a12, srs::sphere<CT>(earth_radius));
0173
0174             geometry::set_from_radian<0>(result, res_direct.lon2);
0175             geometry::set_from_radian<1>(result, res_direct.lat2);
0176
0177             return result;
0178         }
0179         else
0180         {
0181 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
0182             std::cout << "Projection OUTSIDE the segment" << std::endl;
0183 #endif
0184             geometry::detail::conversion::convert_point_to_point(d1 < d2 ? sp1 : sp2, result);
0185             return result;
0186         }
0187     }
0188
0189     template <typename T1, typename T2>
0190     inline radius_type vertical_or_meridian(T1 lat1, T2 lat2) const
0191     {
0192         return m_strategy.radius() * (lat1 - lat2);
0193     }
0194
0195     inline typename Strategy::radius_type radius() const
0196     { return m_strategy.radius(); }
0197
0198 private :
0199     Strategy m_strategy;
0200 };
0201
0202 }} // namespace strategy::closest_points
0203
0204 }} // namespace boost::geometry
0205
0206 #endif // BOOST_GEOMETRY_STRATEGIES_SPHERICAL_CLOSEST_POINTS_CROSS_TRACK_HPP