diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 1700a8842..c3eafc3a7 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -18,6 +18,10 @@ #define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100 #endif +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG +#include +#endif + /*! \brief Algorithm to compute the distance between a segment and a point using direct and inverse geodesic problems as subroutines. The algorithm @@ -33,19 +37,20 @@ template < typename CT, typename Units, - template class Inverse, - template class Direct, + typename FormulaPolicy, + //template class Inverse, + //template class Direct, bool EnableClosestPoint = false > class distance_point_segment{ public: - typedef Inverse inverse_distance_quantities_type; - typedef Inverse inverse_distance_type; - typedef Inverse inverse_azimuth_type; - typedef Inverse inverse_azimuth_reverse_type; - typedef Direct direct_distance_type; + typedef typename FormulaPolicy::template inverse inverse_distance_quantities_type; + typedef typename FormulaPolicy::template inverse inverse_distance_type; + typedef typename FormulaPolicy::template inverse inverse_azimuth_type; + typedef typename FormulaPolicy::template inverse inverse_azimuth_reverse_type; + typedef typename FormulaPolicy::template direct direct_distance_type; struct result_distance_point_segment { @@ -80,9 +85,16 @@ public: CT lon2, CT lat2, //p2 Spheroid const& spheroid) { - CT distance = inverse_distance_quantities_type::apply(lon1, lat1, - lon2, lat2, - spheroid).distance; + //CT distance = inverse_distance_type::apply(lon1, lat1, + //lon2, lat2, + //spheroid).distance; + + //geometry::model::point > point; + // CT distance = geometry::distance(point(lon1, lat1), point(lon2, lat2), geometry::strategy::distance::andoyer<>()); + + //geometry::strategy::distance::geographic<> str(spheroid); + CT distance = geometry::strategy::distance::geographic::apply(lon1, lat1, lon2, lat2, spheroid); + return non_iterative_case(lon1, lat1, distance); } @@ -153,14 +165,25 @@ public: } return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); } - +/* CT d1 = inverse_distance_type::apply(lon1, lat1, lon3, lat3, spheroid).distance; CT d3 = inverse_distance_type::apply(lon1, lat1, lon2, lat2, spheroid).distance; +*/ + + CT d1 = geometry::strategy::distance::geographic + ::apply(lon1, lat1, lon3, lat3, spheroid); + + CT d3 = geometry::strategy::distance::geographic + ::apply(lon1, lat1, lon2, lat2, spheroid); if (geometry::math::equals(d3, c0)) { +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "Degenerate segment" << std::endl; + std::cout << "d1=" << d1 << std::endl; +#endif return non_iterative_case(lon1, lat2, d1); } diff --git a/include/boost/geometry/formulas/elliptic_arc_length.hpp b/include/boost/geometry/formulas/elliptic_arc_length.hpp index 057085eae..e9d7ea7db 100644 --- a/include/boost/geometry/formulas/elliptic_arc_length.hpp +++ b/include/boost/geometry/formulas/elliptic_arc_length.hpp @@ -22,7 +22,6 @@ #include - namespace boost { namespace geometry { namespace formula { @@ -36,6 +35,53 @@ class elliptic_arc_length public : + struct result + { + result() + : distance(0) + , meridian(false) + {} + + CT distance; + bool meridian; + }; + + template + static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid) + { + result res; + + CT c0 = 0; + CT pi = math::pi(); + CT half_pi = pi/CT(2); + CT diff = math::longitude_distance_signed(lon1, lon2); + + if (math::equals(diff, c0)) + { + // single meridian not crossing pole + if (lat1 > lat2) + { + std::swap(lat1, lat2); + } + res.distance = apply(lat2, spheroid) - apply(lat1, spheroid); + res.meridian = true; + } + + if (math::equals(math::abs(diff), pi)) + { + // meridian crosses pole + CT lat_sign = 1; + if (lat1+lat2 < c0) + { + lat_sign = CT(-1); + } + res.distance = math::abs(lat_sign * CT(2) * apply(half_pi, spheroid) + - apply(lat1, spheroid) - apply(lat2, spheroid)); + res.meridian = true; + } + return res; + } + // Distance computation on meridians using series approximations // to elliptic integrals. Formula to compute distance from lattitude 0 to lat // https://en.wikipedia.org/wiki/Meridian_arc @@ -48,7 +94,7 @@ public : CT n = f / (CT(2) - f); CT M = a/(1+n); CT C0 = 1; - + if (Order == 0) { return M * C0 * lat; diff --git a/include/boost/geometry/strategies/geographic/distance.hpp b/include/boost/geometry/strategies/geographic/distance.hpp index 6b195b65e..98ee46a36 100644 --- a/include/boost/geometry/strategies/geographic/distance.hpp +++ b/include/boost/geometry/strategies/geographic/distance.hpp @@ -73,6 +73,29 @@ public : : m_spheroid(spheroid) {} + template + static inline CT apply(CT lon1, CT lat1, CT lon2, CT lat2, + Spheroid const& spheroid) + { + typedef typename formula::elliptic_arc_length + < + CT, strategy::default_order::value + > elliptic_arc_length; + + typename elliptic_arc_length::result res = + elliptic_arc_length::apply(lon1, lat1, lon2, lat2, spheroid); + + if (res.meridian) + { + return res.distance; + } + + return FormulaPolicy::template inverse + < + CT, true, false, false, false, false + >::apply(lon1, lat1, lon2, lat2, spheroid).distance; + } + template inline typename calculation_type::type apply(Point1 const& point1, Point2 const& point2) const @@ -84,45 +107,7 @@ public : CT lon2 = get_as_radian<0>(point2); CT lat2 = get_as_radian<1>(point2); - CT c0 = 0; - CT pi = math::pi(); - CT half_pi = pi/CT(2); - CT diff = math::longitude_distance_signed(lon1, lon2); - - typedef typename formula::elliptic_arc_length - < - CT, strategy::default_order::value - > elliptic_arc_length; - - if (math::equals(diff, c0)) - { - // single meridian not crossing pole - if (lat1 > lat2) - { - std::swap(lat1, lat2); - } - return elliptic_arc_length::apply(lat2, m_spheroid) - - elliptic_arc_length::apply(lat1, m_spheroid); - } - - if (math::equals(math::abs(diff), pi)) - { - // meridian crosses pole - CT lat_sign = 1; - if (lat1+lat2 < c0) - { - lat_sign = CT(-1); - } - return math::abs(lat_sign * CT(2) * - elliptic_arc_length::apply(half_pi, m_spheroid) - - elliptic_arc_length::apply(lat1, m_spheroid) - - elliptic_arc_length::apply(lat2, m_spheroid)); - } - - return FormulaPolicy::template inverse - < - CT, true, false, false, false, false - >::apply(lon1, lat1, lon2, lat2, m_spheroid).distance; + return apply(lon1, lat1, lon2, lat2, m_spheroid); } inline Spheroid const& model() const diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp index 4a220d1c9..872dad358 100644 --- a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -106,8 +106,7 @@ public : < CT, units_type, - FormulaPolicy::template inverse, - FormulaPolicy::template direct + FormulaPolicy >::apply(get<0>(sp1), get<1>(sp1), get<0>(sp2), get<1>(sp2), get<0>(p), get<1>(p), diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 6da51ca55..d3bfd0851 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -8,6 +8,9 @@ // Licensed under the Boost Software License version 1.0. // http://www.boost.org/users/license.html +#define BOOST_GEOMETRY_TEST_DEBUG +#define BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + #include #ifndef BOOST_TEST_MODULE @@ -318,13 +321,12 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, pp_distance("POINT(90 -89)", "POINT(90 90)", strategy_pp), strategy_ps); // degenerate segment - tester::apply("p-s-14", + tester::apply("p-s-15", "POINT(90 90)", "SEGMENT(0.5 -90,175.5 -90)", - pp_distance("POINT(90 -90)", "POINT(90 90)", strategy_pp), + pp_distance("POINT(0.5 -90)", "POINT(90 90)", strategy_pp), strategy_ps); - } //===========================================================================