diff options
Diffstat (limited to 'boost/geometry/strategies/geographic')
8 files changed, 693 insertions, 201 deletions
diff --git a/boost/geometry/strategies/geographic/buffer_point_circle.hpp b/boost/geometry/strategies/geographic/buffer_point_circle.hpp index 7e7c8fe81d..51aacefa82 100644 --- a/boost/geometry/strategies/geographic/buffer_point_circle.hpp +++ b/boost/geometry/strategies/geographic/buffer_point_circle.hpp @@ -88,6 +88,7 @@ public : > direct_t; calculation_type const two_pi = geometry::math::two_pi<calculation_type>(); + calculation_type const pi = geometry::math::pi<calculation_type>(); calculation_type const diff = two_pi / calculation_type(m_count); // TODO: after calculation of some angles is corrected, @@ -96,6 +97,11 @@ public : for (std::size_t i = 0; i < m_count; i++, angle += diff) { + if (angle > pi) + { + angle -= two_pi; + } + typename direct_t::result_type dir_r = direct_t::apply(get_as_radian<0>(point), get_as_radian<1>(point), buffer_distance, angle, diff --git a/boost/geometry/strategies/geographic/disjoint_segment_box.hpp b/boost/geometry/strategies/geographic/disjoint_segment_box.hpp index 49c80a5e78..c4b8430720 100644 --- a/boost/geometry/strategies/geographic/disjoint_segment_box.hpp +++ b/boost/geometry/strategies/geographic/disjoint_segment_box.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2017-2018 Oracle and/or its affiliates. +// Copyright (c) 2017-2019 Oracle and/or its affiliates. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -31,11 +31,12 @@ #include <boost/geometry/srs/spheroid.hpp> +// TODO: spherical_point_box currently defined in the same file as cartesian +#include <boost/geometry/strategies/cartesian/point_in_box.hpp> #include <boost/geometry/strategies/disjoint.hpp> #include <boost/geometry/strategies/geographic/azimuth.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/strategies/normalize.hpp> -#include <boost/geometry/strategies/cartesian/point_in_box.hpp> #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp> @@ -65,22 +66,11 @@ public: : m_spheroid(spheroid) {} - template <typename Segment, typename Box> - struct point_in_geometry_strategy - : services::default_strategy - < - typename point_type<Segment>::type, - Box - > - {}; - - template <typename Segment, typename Box> - static inline typename point_in_geometry_strategy<Segment, Box>::type - get_point_in_geometry_strategy() + typedef covered_by::spherical_point_box disjoint_point_box_strategy_type; + + static inline disjoint_point_box_strategy_type get_disjoint_point_box_strategy() { - typedef typename point_in_geometry_strategy<Segment, Box>::type strategy_type; - - return strategy_type(); + return disjoint_point_box_strategy_type(); } template <typename Segment, typename Box> diff --git a/boost/geometry/strategies/geographic/distance_cross_track.hpp b/boost/geometry/strategies/geographic/distance_cross_track.hpp index 53942529bf..f0ad2b16a1 100644 --- a/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -1,6 +1,6 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) -// Copyright (c) 2016-2018, Oracle and/or its affiliates. +// Copyright (c) 2016-2019, Oracle and/or its affiliates. // Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -34,6 +34,7 @@ #include <boost/geometry/strategies/geographic/azimuth.hpp> #include <boost/geometry/strategies/geographic/distance.hpp> #include <boost/geometry/strategies/geographic/parameters.hpp> +#include <boost/geometry/strategies/geographic/intersection.hpp> #include <boost/geometry/formulas/vincenty_direct.hpp> @@ -44,6 +45,7 @@ #include <boost/geometry/formulas/result_direct.hpp> #include <boost/geometry/formulas/mean_radius.hpp> +#include <boost/geometry/formulas/spherical.hpp> #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK #include <boost/geometry/io/dsv/write.hpp> @@ -63,6 +65,8 @@ namespace boost { namespace geometry namespace strategy { namespace distance { +namespace detail +{ /*! \brief Strategy functor for distance point to segment calculation on ellipsoid Algorithm uses direct and inverse geodesic problems as subroutines. @@ -82,6 +86,7 @@ template typename FormulaPolicy = strategy::andoyer, typename Spheroid = srs::spheroid<double>, typename CalculationType = void, + bool Bisection = false, bool EnableClosestPoint = false > class geographic_cross_track @@ -89,6 +94,29 @@ class geographic_cross_track public : typedef within::spherical_point_point equals_point_point_strategy_type; + typedef intersection::geographic_segments + < + FormulaPolicy, + strategy::default_order<FormulaPolicy>::value, + Spheroid, + CalculationType + > relate_segment_segment_strategy_type; + + inline relate_segment_segment_strategy_type get_relate_segment_segment_strategy() const + { + return relate_segment_segment_strategy_type(m_spheroid); + } + + typedef within::geographic_winding + < + void, void, FormulaPolicy, Spheroid, CalculationType + > point_in_geometry_strategy_type; + + inline point_in_geometry_strategy_type get_point_in_geometry_strategy() const + { + return point_in_geometry_strategy_type(m_spheroid); + } + template <typename Point, typename PointOfSegment> struct return_type : promote_floating_point @@ -112,15 +140,15 @@ public : { typedef typename geometry::detail::cs_angular_units<Point>::type units_type; - return (apply<units_type>(get<0>(sp1), get<1>(sp1), - get<0>(sp2), get<1>(sp2), - get<0>(p), get<1>(p), + return (apply<units_type>(get_as_radian<0>(sp1), get_as_radian<1>(sp1), + get_as_radian<0>(sp2), get_as_radian<1>(sp2), + get_as_radian<0>(p), get_as_radian<1>(p), m_spheroid)).distance; } // points on a meridian not crossing poles template <typename CT> - inline CT vertical_or_meridian(CT lat1, CT lat2) const + inline CT vertical_or_meridian(CT const& lat1, CT const& lat2) const { typedef typename formula::meridian_inverse < @@ -132,6 +160,11 @@ public : m_spheroid); } + Spheroid const& model() const + { + return m_spheroid; + } + private : template <typename CT> @@ -150,7 +183,7 @@ private : template <typename CT> result_distance_point_segment<CT> - static inline non_iterative_case(CT lon, CT lat, CT distance) + static inline non_iterative_case(CT const& lon, CT const& lat, CT const& distance) { result_distance_point_segment<CT> result; result.distance = distance; @@ -165,8 +198,8 @@ private : template <typename CT> result_distance_point_segment<CT> - static inline non_iterative_case(CT lon1, CT lat1, //p1 - CT lon2, CT lat2, //p2 + static inline non_iterative_case(CT const& lon1, CT const& lat1, //p1 + CT const& lon2, CT const& lat2, //p2 Spheroid const& spheroid) { CT distance = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> @@ -176,7 +209,7 @@ private : } template <typename CT> - CT static inline normalize(CT g4, CT& der) + CT static inline normalize(CT const& g4, CT& der) { CT const pi = math::pi<CT>(); if (g4 < -1.25*pi)//close to -270 @@ -205,22 +238,236 @@ private : return g4 - pi/2; } - template <typename Units, typename CT> - result_distance_point_segment<CT> - static inline apply(CT lon1, CT lat1, //p1 - CT lon2, CT lat2, //p2 - CT lon3, CT lat3, //query point p3 - Spheroid const& spheroid) + template <typename CT> + static void bisection(CT const& lon1, CT const& lat1, //p1 + CT const& lon2, CT const& lat2, //p2 + CT const& lon3, CT const& lat3, //query point p3 + Spheroid const& spheroid, + CT const& s14_start, CT const& a12, + result_distance_point_segment<CT>& result) + { + typedef typename FormulaPolicy::template direct<CT, true, false, false, false> + direct_distance_type; + typedef typename FormulaPolicy::template inverse<CT, true, false, false, false, false> + inverse_distance_type; + + geometry::formula::result_direct<CT> res14; + + int counter = 0; // robustness + bool dist_improve = true; + + CT pl_lon = lon1; + CT pl_lat = lat1; + CT pr_lon = lon2; + CT pr_lat = lat2; + + CT s14 = s14_start; + + do{ + // Solve the direct problem to find p4 (GEO) + res14 = direct_distance_type::apply(lon1, lat1, s14, a12, spheroid); + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "dist(pl,p3)=" + << inverse_distance_type::apply(lon3, lat3, + pr_lon, pr_lat, + spheroid).distance + << std::endl; + std::cout << "dist(pr,p3)=" + << inverse_distance_type::apply(lon3, lat3, + pr_lon, pr_lat, + spheroid).distance + << std::endl; +#endif + if (inverse_distance_type::apply(lon3, lat3, pl_lon, pl_lat, spheroid).distance + < inverse_distance_type::apply(lon3, lat3, pr_lon, pr_lat, spheroid).distance) + { + s14 -= inverse_distance_type::apply(res14.lon2, res14.lat2, pl_lon, pl_lat, + spheroid).distance/2; + pr_lon = res14.lon2; + pr_lat = res14.lat2; + } + else + { + s14 += inverse_distance_type::apply(res14.lon2, res14.lat2, pr_lon, pr_lat, + spheroid).distance/2; + pl_lon = res14.lon2; + pl_lat = res14.lat2; + } + + CT new_distance = inverse_distance_type::apply(lon3, lat3, + res14.lon2, res14.lat2, + spheroid).distance; + + dist_improve = new_distance != result.distance; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "p4=" << res14.lon2 * math::r2d<CT>() << + "," << res14.lat2 * math::r2d<CT>() << std::endl; + std::cout << "pl=" << pl_lon * math::r2d<CT>() << "," << pl_lat * math::r2d<CT>()<< std::endl; + std::cout << "pr=" << pr_lon * math::r2d<CT>() << "," << pr_lat * math::r2d<CT>() << std::endl; + std::cout << "new_s14=" << s14 << std::endl; + std::cout << std::setprecision(16) << "result.distance =" << result.distance << std::endl; + std::cout << std::setprecision(16) << "new_distance =" << new_distance << std::endl; + std::cout << "---------end of step " << counter << std::endl<< std::endl; + if (!dist_improve) + { + std::cout << "Stop msg: res34.distance >= prev_distance" << std::endl; + } + if (counter == BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS) + { + std::cout << "Stop msg: counter" << std::endl; + } +#endif + + result.distance = new_distance; + + } while (dist_improve + && counter++ < BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS); + } + + template <typename CT> + static void newton(CT const& lon1, CT const& lat1, //p1 + CT const& lon2, CT const& lat2, //p2 + CT const& lon3, CT const& lat3, //query point p3 + Spheroid const& spheroid, + CT const& s14_start, CT const& a12, + result_distance_point_segment<CT>& result) { typedef typename FormulaPolicy::template inverse<CT, true, true, false, true, true> inverse_distance_azimuth_quantities_type; typedef typename FormulaPolicy::template inverse<CT, false, true, false, false, false> - inverse_azimuth_type; - typedef typename FormulaPolicy::template inverse<CT, false, true, true, false, false> - inverse_azimuth_reverse_type; + inverse_dist_azimuth_type; typedef typename FormulaPolicy::template direct<CT, true, false, false, false> direct_distance_type; + CT const half_pi = math::pi<CT>() / CT(2); + CT prev_distance; + geometry::formula::result_direct<CT> res14; + geometry::formula::result_inverse<CT> res34; + res34.distance = -1; + + int counter = 0; // robustness + CT g4; + CT delta_g4 = 0; + bool dist_improve = true; + CT s14 = s14_start; + + do{ + prev_distance = res34.distance; + + // Solve the direct problem to find p4 (GEO) + res14 = direct_distance_type::apply(lon1, lat1, s14, a12, spheroid); + + // Solve an inverse problem to find g4 + // g4 is the angle between segment (p1,p2) and segment (p3,p4) that meet on p4 (GEO) + + CT a4 = inverse_dist_azimuth_type::apply(res14.lon2, res14.lat2, + lon2, lat2, spheroid).azimuth; + res34 = inverse_distance_azimuth_quantities_type::apply(res14.lon2, res14.lat2, + lon3, lat3, spheroid); + g4 = res34.azimuth - a4; + + CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit + CT m34 = res34.reduced_length; + + if (m34 != 0) + { + CT der = (M43 / m34) * sin(g4); + delta_g4 = normalize(g4, der); + s14 -= der != 0 ? delta_g4 / der : 0; + } + + result.distance = res34.distance; + + dist_improve = prev_distance > res34.distance || prev_distance == -1; + if (!dist_improve) + { + result.distance = prev_distance; + } + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "p4=" << res14.lon2 * math::r2d<CT>() << + "," << res14.lat2 * math::r2d<CT>() << std::endl; + std::cout << "a34=" << res34.azimuth * math::r2d<CT>() << std::endl; + std::cout << "a4=" << a4 * math::r2d<CT>() << std::endl; + std::cout << "g4(normalized)=" << g4 * math::r2d<CT>() << std::endl; + std::cout << "delta_g4=" << delta_g4 * math::r2d<CT>() << std::endl; + std::cout << "der=" << der << std::endl; + std::cout << "M43=" << M43 << std::endl; + std::cout << "m34=" << m34 << std::endl; + std::cout << "new_s14=" << s14 << std::endl; + std::cout << std::setprecision(16) << "dist =" << res34.distance << std::endl; + std::cout << "---------end of step " << counter << std::endl<< std::endl; + if (g4 == half_pi) + { + std::cout << "Stop msg: g4 == half_pi" << std::endl; + } + if (!dist_improve) + { + std::cout << "Stop msg: res34.distance >= prev_distance" << std::endl; + } + if (delta_g4 == 0) + { + std::cout << "Stop msg: delta_g4 == 0" << std::endl; + } + if (counter == BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS) + { + std::cout << "Stop msg: counter" << std::endl; + } +#endif + + } while (g4 != half_pi + && dist_improve + && delta_g4 != 0 + && counter++ < BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS); +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "distance=" << res34.distance << std::endl; + + std::cout << "s34(geo) =" + << inverse_distance_azimuth_quantities_type::apply(res14.lon2, res14.lat2, lon3, lat3, spheroid).distance + << ", p4=(" << res14.lon2 * math::r2d<double>() << "," + << res14.lat2 * math::r2d<double>() << ")" + << std::endl; + + CT s31 = inverse_distance_azimuth_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; + CT s32 = inverse_distance_azimuth_quantities_type::apply(lon3, lat3, lon2, lat2, spheroid).distance; + + CT a4 = inverse_dist_azimuth_type::apply(res14.lon2, res14.lat2, lon2, lat2, spheroid).azimuth; + geometry::formula::result_direct<CT> res4 = direct_distance_type::apply(res14.lon2, res14.lat2, .04, a4, spheroid); + CT p4_plus = inverse_distance_azimuth_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; + + geometry::formula::result_direct<CT> res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a12, spheroid); + CT p4_minus = inverse_distance_azimuth_quantities_type::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; + + std::cout << "s31=" << s31 << "\ns32=" << s32 + << "\np4_plus=" << p4_plus << ", p4=(" << res4.lon2 * math::r2d<double>() << "," << res4.lat2 * math::r2d<double>() << ")" + << "\np4_minus=" << p4_minus << ", p4=(" << res1.lon2 * math::r2d<double>() << "," << res1.lat2 * math::r2d<double>() << ")" + << std::endl; + + if (res34.distance <= p4_plus && res34.distance <= p4_minus) + { + std::cout << "Closest point computed" << std::endl; + } + else + { + std::cout << "There is a closer point nearby" << std::endl; + } +#endif + } + + template <typename Units, typename CT> + result_distance_point_segment<CT> + static inline apply(CT const& lo1, CT const& la1, //p1 + CT const& lo2, CT const& la2, //p2 + CT const& lo3, CT const& la3, //query point p3 + Spheroid const& spheroid) + { + typedef typename FormulaPolicy::template inverse<CT, true, true, false, false, false> + inverse_dist_azimuth_type; + typedef typename FormulaPolicy::template inverse<CT, true, true, true, false, false> + inverse_dist_azimuth_reverse_type; + CT const earth_radius = geometry::formula::mean_radius<CT>(spheroid); result_distance_point_segment<CT> result; @@ -231,13 +478,12 @@ private : CT const half_pi = pi / CT(2); CT const c0 = CT(0); - // Convert to radians - lon1 = math::as_radian<Units>(lon1); - lat1 = math::as_radian<Units>(lat1); - lon2 = math::as_radian<Units>(lon2); - lat2 = math::as_radian<Units>(lat2); - lon3 = math::as_radian<Units>(lon3); - lat3 = math::as_radian<Units>(lat3); + CT lon1 = lo1; + CT lat1 = la1; + CT lon2 = lo2; + CT lat2 = la2; + CT lon3 = lo3; + CT lat3 = la3; if (lon1 > lon2) { @@ -306,49 +552,44 @@ private : std::cout << "Meridian segment crossing pole" << std::endl; #endif CT sign_non_zero = lat3 >= c0 ? 1 : -1; - result_distance_point_segment<CT> d1 = + result_distance_point_segment<CT> res13 = apply<geometry::radian>(lon1, lat1, lon1, half_pi * sign_non_zero, lon3, lat3, spheroid); - result_distance_point_segment<CT> d2 = + result_distance_point_segment<CT> res23 = apply<geometry::radian>(lon2, lat2, lon2, half_pi * sign_non_zero, lon3, lat3, spheroid); - if (d1.distance < d2.distance) + if (res13.distance < res23.distance) { - return d1; + return res13; } else { - return d2; + return res23; } } - CT d1 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> - ::apply(lon1, lat1, lon3, lat3, spheroid); - - CT d3 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> - ::apply(lon1, lat1, lon2, lat2, spheroid); + geometry::formula::result_inverse<CT> res12 = + inverse_dist_azimuth_reverse_type::apply(lon1, lat1, lon2, lat2, spheroid); + geometry::formula::result_inverse<CT> res13 = + inverse_dist_azimuth_type::apply(lon1, lat1, lon3, lat3, spheroid); - if (geometry::math::equals(d3, c0)) + if (geometry::math::equals(res12.distance, c0)) { #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK std::cout << "Degenerate segment" << std::endl; - std::cout << "distance between points=" << d1 << std::endl; + std::cout << "distance between points=" << res13.distance << std::endl; #endif - return non_iterative_case(lon1, lat2, d1); - } + typename meridian_inverse::result res = + meridian_inverse::apply(lon1, lat1, lon3, lat3, spheroid); - CT d2 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> - ::apply(lon2, lat2, lon3, lat3, spheroid); + return non_iterative_case(lon1, lat2, + res.meridian ? res.distance : res13.distance); + } // Compute a12 (GEO) - geometry::formula::result_inverse<CT> res12 = - inverse_azimuth_reverse_type::apply(lon1, lat1, lon2, lat2, spheroid); - CT a12 = res12.azimuth; - CT a13 = inverse_azimuth_type::apply(lon1, lat1, lon3, lat3, spheroid).azimuth; - - CT a312 = a13 - a12; + CT a312 = res13.azimuth - res12.azimuth; // TODO: meridian case optimization if (geometry::math::equals(a312, c0) && meridian_not_crossing_pole) @@ -364,11 +605,11 @@ private : } } - CT projection1 = cos( a312 ) * d1 / d3; + CT projection1 = cos( a312 ) * res13.distance / res12.distance; #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK - std::cout << "a1=" << a12 * math::r2d<CT>() << std::endl; - std::cout << "a13=" << a13 * math::r2d<CT>() << std::endl; + std::cout << "a1=" << res12.azimuth * math::r2d<CT>() << std::endl; + std::cout << "a13=" << res13.azimuth * math::r2d<CT>() << std::endl; std::cout << "a312=" << a312 * math::r2d<CT>() << std::endl; std::cout << "cos(a312)=" << cos(a312) << std::endl; std::cout << "projection 1=" << projection1 << std::endl; @@ -381,19 +622,19 @@ private : #endif // projection of p3 on geodesic spanned by segment (p1,p2) fall // outside of segment on the side of p1 + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } - CT a21 = res12.reverse_azimuth - pi; - CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; - - CT a321 = a23 - a21; + geometry::formula::result_inverse<CT> res23 = + inverse_dist_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid); - CT projection2 = cos( a321 ) * d2 / d3; + CT a321 = res23.azimuth - res12.reverse_azimuth + pi; + CT projection2 = cos( a321 ) * res23.distance / res12.distance; #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK - std::cout << "a21=" << a21 * math::r2d<CT>() << std::endl; - std::cout << "a23=" << a23 * math::r2d<CT>() << std::endl; + std::cout << "a21=" << res12.reverse_azimuth * math::r2d<CT>() << std::endl; + std::cout << "a23=" << res23.azimuth * math::r2d<CT>() << std::endl; std::cout << "a321=" << a321 * math::r2d<CT>() << std::endl; std::cout << "cos(a321)=" << cos(a321) << std::endl; std::cout << "projection 2=" << projection2 << std::endl; @@ -421,139 +662,49 @@ private : point p3 = point(lon3, lat3); geometry::strategy::distance::cross_track<CT> cross_track(earth_radius); - CT s34 = cross_track.apply(p3, p1, p2); + CT s34_sph = cross_track.apply(p3, p1, p2); geometry::strategy::distance::haversine<CT> str(earth_radius); - CT s13 = str.apply(p1, p3); + CT s13_sph = str.apply(p1, p3); //CT s14 = acos( cos(s13/earth_radius) / cos(s34/earth_radius) ) * earth_radius; - CT cos_frac = cos(s13/earth_radius) / cos(s34/earth_radius); - CT s14 = cos_frac >= 1 ? CT(0) - : cos_frac <= -1 ? pi * earth_radius - : acos(cos_frac) * earth_radius; - -#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK - std::cout << "s34=" << s34 << std::endl; - std::cout << "s13=" << s13 << std::endl; - std::cout << "s14=" << s14 << std::endl; - std::cout << "===============" << std::endl; -#endif - - // Update s14 (using Newton method) - CT prev_distance; - geometry::formula::result_direct<CT> res14; - geometry::formula::result_inverse<CT> res34; - res34.distance = -1; + CT cos_frac = cos(s13_sph / earth_radius) / cos(s34_sph / earth_radius); + CT s14_sph = cos_frac >= 1 ? CT(0) + : cos_frac <= -1 ? pi * earth_radius + : acos(cos_frac) * earth_radius; - int counter = 0; // robustness - CT g4; - CT delta_g4; - bool dist_improve = true; + CT a12_sph = geometry::formula::spherical_azimuth<>(lon1, lat1, lon2, lat2); - do{ - prev_distance = res34.distance; + geometry::formula::result_direct<CT> res + = geometry::formula::spherical_direct<true, false> + (lon1, lat1, s14_sph, a12_sph, srs::sphere<CT>(earth_radius)); - // Solve the direct problem to find p4 (GEO) - res14 = direct_distance_type::apply(lon1, lat1, s14, a12, spheroid); - - // Solve an inverse problem to find g4 - // g4 is the angle between segment (p1,p2) and segment (p3,p4) that meet on p4 (GEO) - - CT a4 = inverse_azimuth_type::apply(res14.lon2, res14.lat2, - lon2, lat2, spheroid).azimuth; - res34 = inverse_distance_azimuth_quantities_type::apply(res14.lon2, res14.lat2, - lon3, lat3, spheroid); - g4 = res34.azimuth - a4; - - CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit - CT m34 = res34.reduced_length; - CT der = (M43 / m34) * sin(g4); - - //normalize - delta_g4 = normalize(g4, der); - s14 = s14 - delta_g4 / der; - result.distance = res34.distance; - - dist_improve = prev_distance > res34.distance || prev_distance == -1; - if (!dist_improve) - { - result.distance = prev_distance; - } + // this is what postgis (version 2.5) returns + // geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> + // ::apply(lon3, lat3, res.lon2, res.lat2, spheroid); #ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK - std::cout << "p4=" << res14.lon2 * math::r2d<CT>() << - "," << res14.lat2 * math::r2d<CT>() << std::endl; - std::cout << "a34=" << res34.azimuth * math::r2d<CT>() << std::endl; - std::cout << "a4=" << a4 * math::r2d<CT>() << std::endl; - std::cout << "g4(normalized)=" << g4 * math::r2d<CT>() << std::endl; - std::cout << "delta_g4=" << delta_g4 * math::r2d<CT>() << std::endl; - std::cout << "der=" << der << std::endl; - std::cout << "M43=" << M43 << std::endl; - std::cout << "spherical limit=" << cos(s14/earth_radius) << std::endl; - std::cout << "m34=" << m34 << std::endl; - std::cout << "new_s14=" << s14 << std::endl; - std::cout << std::setprecision(16) << "dist =" << res34.distance << std::endl; - std::cout << "---------end of step " << counter << std::endl<< std::endl; - if (g4 == half_pi) - { - std::cout << "Stop msg: g4 == half_pi" << std::endl; - } - if (!dist_improve) - { - std::cout << "Stop msg: res34.distance >= prev_distance" << std::endl; - } - if (delta_g4 == 0) - { - std::cout << "Stop msg: delta_g4 == 0" << std::endl; - } - if (counter == BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS) - { - std::cout << "Stop msg: counter" << std::endl; - } + std::cout << "s34=" << s34_sph << std::endl; + std::cout << "s13=" << res13.distance << std::endl; + std::cout << "s14=" << s14_sph << std::endl; + std::cout << "===============" << std::endl; #endif - } while (g4 != half_pi - && dist_improve - && delta_g4 != 0 - && counter++ < BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS); - -#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK - std::cout << "distance=" << res34.distance << std::endl; - - point p4(res14.lon2, res14.lat2); - CT s34_sph = str.apply(p4, p3); - - std::cout << "s34(sph) =" << s34_sph << std::endl; - std::cout << "s34(geo) =" - << inverse_distance_azimuth_quantities_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance - << ", p4=(" << get<0>(p4) * math::r2d<double>() << "," - << get<1>(p4) * math::r2d<double>() << ")" - << std::endl; - - CT s31 = inverse_distance_azimuth_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; - CT s32 = inverse_distance_azimuth_quantities_type::apply(lon3, lat3, lon2, lat2, spheroid).distance; - - CT a4 = inverse_azimuth_type::apply(get<0>(p4), get<1>(p4), lon2, lat2, spheroid).azimuth; - geometry::formula::result_direct<CT> res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); - CT p4_plus = inverse_distance_azimuth_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; - - geometry::formula::result_direct<CT> res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a12, spheroid); - CT p4_minus = inverse_distance_azimuth_quantities_type::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; - - std::cout << "s31=" << s31 << "\ns32=" << s32 - << "\np4_plus=" << p4_plus << ", p4=(" << res4.lon2 * math::r2d<double>() << "," << res4.lat2 * math::r2d<double>() << ")" - << "\np4_minus=" << p4_minus << ", p4=(" << res1.lon2 * math::r2d<double>() << "," << res1.lat2 * math::r2d<double>() << ")" - << std::endl; + // Update s14 (using Newton method) - if (res34.distance <= p4_plus && res34.distance <= p4_minus) + if (Bisection) { - std::cout << "Closest point computed" << std::endl; + bisection<CT>(lon1, lat1, lon2, lat2, lon3, lat3, + spheroid, res12.distance/2, res12.azimuth, result); } else { - std::cout << "There is a closer point nearby" << std::endl; + CT s14_start = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT> + ::apply(lon1, lat1, res.lon2, res.lat2, spheroid); + + newton<CT>(lon1, lat1, lon2, lat2, lon3, lat3, + spheroid, s14_start, res12.azimuth, result); } -#endif return result; } @@ -561,7 +712,37 @@ private : Spheroid m_spheroid; }; +} // namespace detail + +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid<double>, + typename CalculationType = void +> +class geographic_cross_track + : public detail::geographic_cross_track + < + FormulaPolicy, + Spheroid, + CalculationType, + false, + false + > +{ +public : + explicit geographic_cross_track(Spheroid const& spheroid = Spheroid()) + : + detail::geographic_cross_track< + FormulaPolicy, + Spheroid, + CalculationType, + false, + false + >(spheroid) + {} +}; #ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS namespace services @@ -595,6 +776,17 @@ struct tag<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> > typedef strategy_tag_distance_point_segment type; }; +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + bool Bisection +> +struct tag<detail::geographic_cross_track<FormulaPolicy, Spheroid, CalculationType, Bisection> > +{ + typedef strategy_tag_distance_point_segment type; +}; //return types template <typename FormulaPolicy, typename P, typename PS> @@ -625,6 +817,19 @@ struct return_type<geographic_cross_track<FormulaPolicy, Spheroid, CalculationTy : geographic_cross_track<FormulaPolicy, Spheroid, CalculationType>::template return_type<P, PS> {}; +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + bool Bisection, + typename P, + typename PS +> +struct return_type<detail::geographic_cross_track<FormulaPolicy, Spheroid, CalculationType, Bisection>, P, PS> + : detail::geographic_cross_track<FormulaPolicy, Spheroid, CalculationType, Bisection>::template return_type<P, PS> +{}; + //comparable types template < @@ -640,6 +845,22 @@ struct comparable_type<geographic_cross_track<FormulaPolicy, Spheroid, Calculati > type; }; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + bool Bisection +> +struct comparable_type<detail::geographic_cross_track<FormulaPolicy, Spheroid, CalculationType, Bisection> > +{ + typedef detail::geographic_cross_track + < + FormulaPolicy, Spheroid, CalculationType, Bisection + > type; +}; + template < typename FormulaPolicy, @@ -656,6 +877,23 @@ public : } }; +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + bool Bisection +> +struct get_comparable<detail::geographic_cross_track<FormulaPolicy, Spheroid, CalculationType, Bisection> > +{ +public : + static inline detail::geographic_cross_track<FormulaPolicy, Spheroid, CalculationType, Bisection> + apply(detail::geographic_cross_track<FormulaPolicy, Spheroid, CalculationType, Bisection> const& strategy) + { + return strategy; + } +}; + template < diff --git a/boost/geometry/strategies/geographic/distance_segment_box.hpp b/boost/geometry/strategies/geographic/distance_segment_box.hpp index 615f938739..99b181fa25 100644 --- a/boost/geometry/strategies/geographic/distance_segment_box.hpp +++ b/boost/geometry/strategies/geographic/distance_segment_box.hpp @@ -16,9 +16,10 @@ #include <boost/geometry/algorithms/detail/distance/segment_to_box.hpp> #include <boost/geometry/strategies/distance.hpp> -#include <boost/geometry/strategies/geographic/parameters.hpp> #include <boost/geometry/strategies/geographic/azimuth.hpp> #include <boost/geometry/strategies/geographic/distance_cross_track.hpp> +#include <boost/geometry/strategies/geographic/parameters.hpp> +#include <boost/geometry/strategies/geographic/side.hpp> #include <boost/geometry/strategies/normalize.hpp> #include <boost/geometry/strategies/spherical/disjoint_box_box.hpp> #include <boost/geometry/strategies/spherical/distance_segment_box.hpp> @@ -55,6 +56,8 @@ struct geographic_segment_box > {}; + typedef geographic_tag cs_tag; + // point-point strategy getters struct distance_pp_strategy { @@ -83,6 +86,33 @@ struct geographic_segment_box return distance_type(m_spheroid); } + struct distance_pb_strategy + { + typedef geographic_cross_track_point_box + < + FormulaPolicy, + Spheroid, + CalculationType + > type; + }; + + inline typename distance_pb_strategy::type get_distance_pb_strategy() const + { + return typename distance_pb_strategy::type(m_spheroid); + } + + typedef side::geographic + < + FormulaPolicy, + Spheroid, + CalculationType + > side_strategy_type; + + inline side_strategy_type get_side_strategy() const + { + return side_strategy_type(m_spheroid); + } + typedef within::spherical_point_point equals_point_point_strategy_type; static inline equals_point_point_strategy_type get_equals_point_point_strategy() diff --git a/boost/geometry/strategies/geographic/index.hpp b/boost/geometry/strategies/geographic/index.hpp new file mode 100644 index 0000000000..202f0620df --- /dev/null +++ b/boost/geometry/strategies/geographic/index.hpp @@ -0,0 +1,216 @@ +// Boost.Geometry Index +// +// R-tree strategies +// +// Copyright (c) 2019, Oracle and/or its affiliates. +// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle +// +// Use, modification and distribution is subject to the Boost Software License, +// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) + +#ifndef BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INDEX_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INDEX_HPP + + +#include <boost/geometry/strategies/geographic/distance.hpp> +#include <boost/geometry/strategies/geographic/distance_andoyer.hpp> // backward compatibility +#include <boost/geometry/strategies/geographic/distance_cross_track.hpp> +#include <boost/geometry/strategies/geographic/distance_cross_track_point_box.hpp> +#include <boost/geometry/strategies/geographic/distance_segment_box.hpp> +#include <boost/geometry/strategies/geographic/distance_thomas.hpp> // backward compatibility +#include <boost/geometry/strategies/geographic/distance_vincenty.hpp> // backward compatibility +#include <boost/geometry/strategies/geographic/envelope_segment.hpp> +#include <boost/geometry/strategies/geographic/expand_segment.hpp> +#include <boost/geometry/strategies/geographic/intersection.hpp> +#include <boost/geometry/strategies/geographic/point_in_poly_winding.hpp> + +#include <boost/geometry/strategies/spherical/index.hpp> + + +namespace boost { namespace geometry { namespace strategy { namespace index +{ + +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = geometry::srs::spheroid<double>, + typename CalculationType = void +> +struct geographic + : spherical<CalculationType> +{ + typedef geographic_tag cs_tag; + + typedef geometry::strategy::envelope::geographic_segment + < + FormulaPolicy, Spheroid, CalculationType + > envelope_segment_strategy_type; + + inline envelope_segment_strategy_type get_envelope_segment_strategy() const + { + return envelope_segment_strategy_type(m_spheroid); + } + + typedef geometry::strategy::expand::geographic_segment + < + FormulaPolicy, Spheroid, CalculationType + > expand_segment_strategy_type; + + inline expand_segment_strategy_type get_expand_segment_strategy() const + { + return expand_segment_strategy_type(m_spheroid); + } + + // used in equals(Seg, Seg) but only to get_point_in_point_strategy() + typedef geometry::strategy::intersection::geographic_segments + < + FormulaPolicy, + // If index::geographic formula is derived from intersection::geographic_segments + // formula with different Order this may cause an inconsistency + strategy::default_order<FormulaPolicy>::value, + Spheroid, + CalculationType + > relate_segment_segment_strategy_type; + + inline relate_segment_segment_strategy_type get_relate_segment_segment_strategy() const + { + return relate_segment_segment_strategy_type(m_spheroid); + } + + typedef geometry::strategy::distance::geographic + < + FormulaPolicy, Spheroid, CalculationType + > comparable_distance_point_point_strategy_type; + + inline comparable_distance_point_point_strategy_type get_comparable_distance_point_point_strategy() const + { + return comparable_distance_point_point_strategy_type(m_spheroid); + } + + typedef geometry::strategy::distance::geographic_cross_track_point_box + < + FormulaPolicy, Spheroid, CalculationType + > comparable_distance_point_box_strategy_type; + + inline comparable_distance_point_box_strategy_type get_comparable_distance_point_box_strategy() const + { + return comparable_distance_point_box_strategy_type(m_spheroid); + } + + typedef geometry::strategy::distance::geographic_cross_track + < + FormulaPolicy, Spheroid, CalculationType + > comparable_distance_point_segment_strategy_type; + + inline comparable_distance_point_segment_strategy_type get_comparable_distance_point_segment_strategy() const + { + return comparable_distance_point_segment_strategy_type(m_spheroid); + } + + typedef geometry::strategy::distance::geographic_segment_box + < + FormulaPolicy, Spheroid, CalculationType + > comparable_distance_segment_box_strategy_type; + + inline comparable_distance_segment_box_strategy_type get_comparable_distance_segment_box_strategy() const + { + return comparable_distance_segment_box_strategy_type(m_spheroid); + } + + geographic() + : m_spheroid() + {} + + explicit geographic(Spheroid const& spheroid) + : m_spheroid(spheroid) + {} + +public: + Spheroid m_spheroid; +}; + + +namespace services +{ + +template <typename Geometry> +struct default_strategy<Geometry, geographic_tag> +{ + typedef geographic<> type; +}; + + +// within and relate (MPt, Mls/MPoly) +template <typename Point1, typename Point2, typename Formula, typename Spheroid, typename CalculationType> +struct from_strategy<within::geographic_winding<Point1, Point2, Formula, Spheroid, CalculationType> > +{ + typedef strategy::index::geographic<Formula, Spheroid, CalculationType> type; + + static inline type get(within::geographic_winding<Point1, Point2, Formula, Spheroid, CalculationType> const& strategy) + { + return type(strategy.model()); + } +}; + +// distance (MPt, MPt) +template <typename Formula, typename Spheroid, typename CalculationType> +struct from_strategy<distance::geographic<Formula, Spheroid, CalculationType> > +{ + typedef strategy::index::geographic<Formula, Spheroid, CalculationType> type; + + static inline type get(distance::geographic<Formula, Spheroid, CalculationType> const& strategy) + { + return type(strategy.model()); + } +}; +template <typename Spheroid, typename CalculationType> +struct from_strategy<distance::andoyer<Spheroid, CalculationType> > +{ + typedef strategy::index::geographic<strategy::andoyer, Spheroid, CalculationType> type; + + static inline type get(distance::andoyer<Spheroid, CalculationType> const& strategy) + { + return type(strategy.model()); + } +}; +template <typename Spheroid, typename CalculationType> +struct from_strategy<distance::thomas<Spheroid, CalculationType> > +{ + typedef strategy::index::geographic<strategy::thomas, Spheroid, CalculationType> type; + + static inline type get(distance::thomas<Spheroid, CalculationType> const& strategy) + { + return type(strategy.model()); + } +}; +template <typename Spheroid, typename CalculationType> +struct from_strategy<distance::vincenty<Spheroid, CalculationType> > +{ + typedef strategy::index::geographic<strategy::vincenty, Spheroid, CalculationType> type; + + static inline type get(distance::vincenty<Spheroid, CalculationType> const& strategy) + { + return type(strategy.model()); + } +}; + +// distance (MPt, Linear/Areal) +template <typename Formula, typename Spheroid, typename CalculationType> +struct from_strategy<distance::geographic_cross_track<Formula, Spheroid, CalculationType> > +{ + typedef strategy::index::geographic<Formula, Spheroid, CalculationType> type; + + static inline type get(distance::geographic_cross_track<Formula, Spheroid, CalculationType> const& strategy) + { + return type(strategy.model()); + } +}; + + +} // namespace services + + +}}}} // namespace boost::geometry::strategy::index + +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_INDEX_HPP diff --git a/boost/geometry/strategies/geographic/intersection.hpp b/boost/geometry/strategies/geographic/intersection.hpp index 37f79e5895..1d69933502 100644 --- a/boost/geometry/strategies/geographic/intersection.hpp +++ b/boost/geometry/strategies/geographic/intersection.hpp @@ -74,6 +74,8 @@ template > struct geographic_segments { + typedef geographic_tag cs_tag; + typedef side::geographic < FormulaPolicy, Spheroid, CalculationType @@ -193,6 +195,9 @@ struct geographic_segments } typedef covered_by::spherical_point_box disjoint_point_box_strategy_type; + typedef covered_by::spherical_point_box covered_by_point_box_strategy_type; + typedef within::spherical_point_box within_point_box_strategy_type; + typedef envelope::spherical_box envelope_box_strategy_type; typedef expand::spherical_box expand_box_strategy_type; enum intersection_point_flag { ipi_inters = 0, ipi_at_a1, ipi_at_a2, ipi_at_b1, ipi_at_b2 }; @@ -494,7 +499,6 @@ private: } // NOTE: this is probably not needed - calc_t const c0 = 0; int a1_on_b = position_value(c0, dist_a1_b1, dist_a1_b2); int a2_on_b = position_value(dist_a1_a2, dist_a1_b1, dist_a1_b2); int b1_on_a = position_value(c0, dist_b1_a1, dist_b1_a2); diff --git a/boost/geometry/strategies/geographic/point_in_poly_winding.hpp b/boost/geometry/strategies/geographic/point_in_poly_winding.hpp index 95a1961474..fe45f8db7a 100644 --- a/boost/geometry/strategies/geographic/point_in_poly_winding.hpp +++ b/boost/geometry/strategies/geographic/point_in_poly_winding.hpp @@ -1,6 +1,6 @@ // Boost.Geometry -// Copyright (c) 2017 Oracle and/or its affiliates. +// Copyright (c) 2017, 2019 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, @@ -38,8 +38,8 @@ namespace strategy { namespace within */ template < - typename Point, - typename PointOfSegment = Point, + typename Point = void, // for backward compatibility + typename PointOfSegment = Point, // for backward compatibility typename FormulaPolicy = strategy::andoyer, typename Spheroid = srs::spheroid<double>, typename CalculationType = void @@ -47,16 +47,12 @@ template class geographic_winding : public within::detail::spherical_winding_base < - Point, - PointOfSegment, side::geographic<FormulaPolicy, Spheroid, CalculationType>, CalculationType > { typedef within::detail::spherical_winding_base < - Point, - PointOfSegment, side::geographic<FormulaPolicy, Spheroid, CalculationType>, CalculationType > base_t; @@ -68,6 +64,11 @@ public: explicit geographic_winding(Spheroid const& model) : base_t(model) {} + + Spheroid const& model() const + { + return base_t::m_side_strategy.model(); + } }; diff --git a/boost/geometry/strategies/geographic/side.hpp b/boost/geometry/strategies/geographic/side.hpp index e22718f32b..1a9648e74c 100644 --- a/boost/geometry/strategies/geographic/side.hpp +++ b/boost/geometry/strategies/geographic/side.hpp @@ -2,8 +2,8 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014-2018. -// Modifications copyright (c) 2014-2018 Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2019. +// Modifications copyright (c) 2014-2019 Oracle and/or its affiliates. // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle @@ -65,6 +65,8 @@ template class geographic { public: + typedef geographic_tag cs_tag; + typedef strategy::envelope::geographic < FormulaPolicy, @@ -123,6 +125,11 @@ public: return formula::azimuth_side_value(a1p, a12); } + Spheroid const& model() const + { + return m_model; + } + private: template <typename ResultType, typename InverseFormulaType, |