From 8b33899ec7b9b09322f442d8aef6c4bb342b7442 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Thu, 13 Jul 2017 09:59:10 +0300 Subject: [PATCH 01/27] [distance] Point-segment distance formula, strategy and tests --- .../detail/distance/point_to_geometry.hpp | 2 +- .../formulas/distance_point_segment.hpp | 288 ++++++++++++++ .../boost/geometry/strategies/distance.hpp | 3 - .../geographic/distance_cross_track.hpp | 339 ++++++++++++++++ .../strategies/geographic/parameters.hpp | 51 +++ .../boost/geometry/strategies/strategies.hpp | 1 + test/algorithms/distance/Jamfile.v2 | 1 + .../distance/distance_brute_force.hpp | 22 +- .../algorithms/distance/distance_geo_pl_l.cpp | 370 +++++++++++++++++ .../distance/test_distance_geo_common.hpp | 371 ++++++++++++++++++ 10 files changed, 1433 insertions(+), 15 deletions(-) create mode 100644 include/boost/geometry/formulas/distance_point_segment.hpp create mode 100644 include/boost/geometry/strategies/geographic/distance_cross_track.hpp create mode 100644 test/algorithms/distance/distance_geo_pl_l.cpp create mode 100644 test/algorithms/distance/test_distance_geo_common.hpp diff --git a/include/boost/geometry/algorithms/detail/distance/point_to_geometry.hpp b/include/boost/geometry/algorithms/detail/distance/point_to_geometry.hpp index ab5de3d9b..9596799cb 100644 --- a/include/boost/geometry/algorithms/detail/distance/point_to_geometry.hpp +++ b/include/boost/geometry/algorithms/detail/distance/point_to_geometry.hpp @@ -391,7 +391,7 @@ struct distance template struct distance < - Point, Polygon, Strategy, point_tag, polygon_tag, + Point, Polygon, Strategy, point_tag, polygon_tag, strategy_tag_distance_point_segment, false > : detail::distance::point_to_polygon < diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp new file mode 100644 index 000000000..6363aac8b --- /dev/null +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -0,0 +1,288 @@ +// Boost.Geometry + +// Copyright (c) 2016-2017 Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, 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_FORMULAS_CROSS_TRACK_GEO_HPP +#define BOOST_GEOMETRY_FORMULAS_CROSS_TRACK_GEO_HPP + +#include +#include + +namespace boost { namespace geometry { namespace formula +{ + +template +< + typename CT, + typename Units, + typename Inverse_type_azimuth, + typename Inverse_type_distance, + typename Direct_type +> +class cross_track_geo_formula{ + +public: + + CT static inline apply(CT lon1, CT lat1, //p1 + CT lon2, CT lat2, //p2 + CT lon3, CT lat3) //query point p3 + { + + int print = 0; + + // Constants + geometry::srs::spheroid spheroid; + CT const f = flattening(spheroid); + CT const pi = math::pi(); + CT const half_pi = math::pi() / CT(2); + geometry::srs::sphere sphere; + CT const earth_radius = sphere.get_radius<1>(); + //CT const earth_radius = 6372795; + //CT const earth_radius = 6378137; + + // Convert to radians + lon1 = math::as_radian(lon1); + lat1 = math::as_radian(lat1); + lon2 = math::as_radian(lon2); + lat2 = math::as_radian(lat2); + lon3 = math::as_radian(lon3); + lat3 = math::as_radian(lat3); + + if (lon1 > lon2) + { + std::swap(lon1, lon2); + std::swap(lat1, lat2); + } + + //segment on equator + if (math::equals(lat1, 0) && math::equals(lat2, 0)) + { + if (lon3 <= lon1) + { + return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + } + if (lon3 >= lon2) + { + return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + } + return Inverse_type_distance::apply(lon3, lat1, lon3, lat3, spheroid).distance; + } + + //segment on meridian + if (math::equals(lon1, lon2)) + { + if (lat1 > lat2) + { + std::swap(lat1, lat2); + } + if (lat3 <= lat1) + { + return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + } + if (lat3 >= lat2) + { + return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + } + return Inverse_type_distance::apply(lon1, lat3, lon3, lat3, spheroid).distance; + } + + // Easy cases + if (lon3 <= lon1 && lat3 <= lat1) + { + return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + } + if (lon3 >= lon2 && lat3 <= lat2) + { + return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + } + + // Compute a1 (GEO) + CT a1 = Inverse_type_azimuth::apply(lon1, lat1, lon2, lat2, spheroid).azimuth; + CT a13 = Inverse_type_azimuth::apply(lon1, lat1, lon3, lat3, spheroid).azimuth; + + CT a312 = a1 - a13; + + if (print) + { + std::cout << "a1=" << a1 * math::r2d() << std::endl; + std::cout << "a13=" << a13 * math::r2d() << std::endl; + std::cout << "a312=" << a312 * math::r2d() << std::endl; + } + if (a312 > half_pi) + { + // projection of p3 on geodesic spanned by segment (p1,p2) fall + // outside of segment on the side of p1 + return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + } + + CT a2 = pi + Inverse_type_azimuth::apply(lon2, lat2, lon1, lat1, spheroid).azimuth; + CT a23 = Inverse_type_azimuth::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; + + CT a321 = a2 - a23; + + if (print) + { + std::cout << "a2=" << a2 * math::r2d() << std::endl; + std::cout << "a23=" << a13 * math::r2d() << std::endl; + std::cout << "a321=" << a321 * math::r2d() << std::endl; + } + if (a321 < half_pi) + { + // projection of p3 on geodesic spanned by segment (p1,p2) fall + // outside of segment on the side of p2 + return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + } + + // Guess s14 (SPHERICAL) + typedef geometry::model::point > point; + + CT bet1 = atan((1 - f) * tan(lon1)); + CT bet2 = atan((1 - f) * tan(lon2)); + CT bet3 = atan((1 - f) * tan(lon3)); + point p1 = point(bet1, lat1); + point p2 = point(bet2, lat2); + point p3 = point(bet3, lat3); + + geometry::strategy::distance::cross_track cross_track(earth_radius); + CT s34 = cross_track.apply(p3, p1, p2); + + geometry::strategy::distance::haversine str(earth_radius); + CT s13 = str.apply(p1, p3); + /* + // Closest point is one of the segment's endpoints + if (s34 == s13) + { + return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + } + CT s23 = str.apply(p2, p3); + if (s34 == s23) + { + return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + } +*/ + CT s14 = acos( cos(s13/earth_radius) / cos(s34/earth_radius)) * earth_radius; + + if (print) + { + std::cout << "s34=" << s34 << std::endl; + std::cout << "s13=" << s13 << std::endl; + std::cout << "s14=" << s14 << std::endl; + std::cout << "===============" << std::endl; + } + + // Update s14 (using Newton method) + CT prev_distance = 0; + geometry::formula::result_direct res14; + geometry::formula::result_inverse res34; + + for (int i=0; i<20; ++i) + { + // Solve the direct problem to find p4 (GEO) + res14 = Direct_type::apply(lon1, lat1, s14, a1, 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_type_azimuth::apply(res14.lon2, res14.lat2, lon2, lat2, spheroid).azimuth; + res34 = Inverse_type_distance::apply(res14.lon2, res14.lat2, lon3, lat3, spheroid); + CT g4 = res34.azimuth - a4; + + // Normalize g4 + CT delta_g4; + + if (g4 < 0 && g4 < -pi)//close to -270 + { + delta_g4 = g4 + 1.5 * pi; + } + else if (g4 > 0 && g4 > pi)//close to 270 + { + delta_g4 = - g4 + 1.5 * pi; + } + else if (g4 < 0 && g4 > -pi)//close to -90 + { + delta_g4 = -g4 - pi/2; + } + else //close to 90 + { + delta_g4 = g4 - pi/2; + } + + CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit + CT m34 = res34.reduced_length; + CT der = M43 / m34; + s14 = s14 - delta_g4 / der; + + //Print + if (print) + { + std::cout << "p4=" << res14.lon2 * math::r2d() << + "," << res14.lat2 * math::r2d() << std::endl; + std::cout << "delta_g4=" << delta_g4 << std::endl; + std::cout << "g4=" << g4 * math::r2d() << std::endl; + std::cout << "der=" << der << std::endl; + std::cout << "new_s14=" << s14 << std::endl; + std::cout << "dist =" << res34.distance << std::endl; + std::cout << "---------end of step " << i << std::endl<< std::endl; + } + + if (g4 == math::pi()/2 || (res34.distance >= prev_distance && prev_distance!=0) || delta_g4 == 0) + { + //std::cout << "stop" << std::endl; + return prev_distance; + } + prev_distance = res34.distance; + + } + + ////////////////////////////TESTING + + std::cout << "dist =" << 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_type_distance::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance + << ", p4=(" << get<0>(p4) * math::r2d() << "," + << get<1>(p4) * math::r2d() << ")" + << std::endl; + + CT s31 = Inverse_type_distance::apply(lon3, lat3, lon1, lat1, spheroid).distance; + CT s32 = Inverse_type_distance::apply(lon3, lat3, lon2, lat2, spheroid).distance; + + CT a4 = Inverse_type_azimuth::apply(get<0>(p4), get<1>(p4), lon2, lat2, spheroid).azimuth; + geometry::formula::result_direct res4 = Direct_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); + CT p4_plus = Inverse_type_distance::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; + + geometry::formula::result_direct res1 = Direct_type::apply(lon1, lat1, s14-.04, a1, spheroid); + CT p4_minus = Inverse_type_distance::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; + + std::cout << "s31=" << s31 << "\ns32=" << s32 + << "\np4_plus=" << p4_plus << ", p4=(" << res4.lon2 * math::r2d() << "," << res4.lat2 * math::r2d() << ")" + << "\np4_minus=" << p4_minus << ", p4=(" << res1.lon2 * math::r2d() << "," << res1.lat2 * math::r2d() << ")" + << std::endl; + + if (res34.distance <= p4_plus && res34.distance <= p4_minus) + { + std::cout << "Closest point computed" << std::endl << std::endl; + } + else + { + std::cout << "There is a closer point nearby!" << std::endl; + } + + CT min_dist_endpoints = std::min(s31,s32); + return std::min(res34.distance, min_dist_endpoints); + } + +}; + +}}} // namespace boost::geometry::formula +#endif // BOOST_GEOMETRY_FORMULAS_CROSS_TRACK_GEO_HPP diff --git a/include/boost/geometry/strategies/distance.hpp b/include/boost/geometry/strategies/distance.hpp index 98ccb8202..0d5a9e6dd 100644 --- a/include/boost/geometry/strategies/distance.hpp +++ b/include/boost/geometry/strategies/distance.hpp @@ -65,9 +65,6 @@ template struct get_comparable template struct result_from_distance {}; - - - // Default strategy diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp new file mode 100644 index 000000000..b339f8f45 --- /dev/null +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -0,0 +1,339 @@ +// Boost.Geometry (aka GGL, Generic Geometry Library) + +// Copyright (c) 2016-2017, Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, 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_DISTANCE_CROSS_TRACK_HPP +#define BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_HPP + +#include + +#include +#include +#include +#include + +#include +#include +#include +#include + +#include + +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK +# include +#endif + + +namespace boost { namespace geometry +{ + +namespace strategy { namespace distance +{ + + +/*! +\brief Strategy functor for distance point to segment calculation +\ingroup strategies +\details Class which calculates the distance of a point to a segment, for points on a sphere or globe +\see http://williams.best.vwh.net/avform.htm +\tparam CalculationType \tparam_calculation +\tparam Strategy underlying point-point distance strategy +\qbk{ +[heading See also] +[link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] +} +*/ +template +< + typename FormulaPolicy = strategy::andoyer, + typename Spheroid = srs::spheroid, + typename CalculationType = void +> +class cross_track_geo +{ +public : + template + struct return_type + : promote_floating_point + < + typename select_calculation_type + < + Point, + PointOfSegment, + CalculationType + >::type + > + {}; + + inline cross_track_geo() + {} + + typedef geographic + < + FormulaPolicy, Spheroid, CalculationType + > DistanceStrategy; + + /* + template + inline T comparable_to_distance(T& a) const + { + return a; + } +*/ + template + inline typename return_type::type + apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const + { + typedef typename coordinate_system::type::units units_type; + + typedef typename return_type::type CT; + + typedef typename FormulaPolicy::template inverse + < + CT, + false, + true, + false, + false, + false + > inverse_type_azimuth; + + typedef typename FormulaPolicy::template inverse + < + CT, + true, + true, + false, + true, + true + > inverse_type_distance; + + typedef typename FormulaPolicy::template direct + < + CT, + true, + false, + false, + false + > direct_type; + + return geometry::formula::cross_track_geo_formula + < + CT, + units_type, + inverse_type_azimuth, + inverse_type_distance, + direct_type + >::apply(get<0>(sp1), get<1>(sp1), + get<0>(sp2), get<1>(sp2), + get<0>(p), get<1>(p)); + } + + //TODO: apply a more general strategy getter + inline DistanceStrategy get_distance_strategy() const + { + return DistanceStrategy(); + } + +}; + + + +#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS +namespace services +{ + +//tags +template +struct tag > +{ + typedef strategy_tag_distance_point_segment type; +}; + +template +< + typename FormulaPolicy, + typename Spheroid +> +struct tag > +{ + typedef strategy_tag_distance_point_segment type; +}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct tag > +{ + typedef strategy_tag_distance_point_segment type; +}; + + +//return types +template +struct return_type, P, PS> + : cross_track_geo::template return_type +{}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename P, + typename PS +> +struct return_type, P, PS> + : cross_track_geo::template return_type +{}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P, + typename PS +> +struct return_type, P, PS> + : cross_track_geo::template return_type +{}; + +//comparable types +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct comparable_type > +{ + typedef cross_track_geo + < + FormulaPolicy, Spheroid, CalculationType + > type; +}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct get_comparable > +{ + typedef typename comparable_type + < + cross_track_geo + >::type comparable_type; +public : + static inline comparable_type + apply(cross_track_geo const& strategy) + { + return comparable_type(); + } +}; + + +template +< + typename FormulaPolicy, + typename P, + typename PS +> +struct result_from_distance, P, PS> +{ +private : + typedef typename cross_track_geo + < + FormulaPolicy + >::template return_type::type return_type; +public : + template + static inline return_type + apply(cross_track_geo const& , T const& distance) + { + return distance; + } +}; + +/* +template +< + typename FormulaPolicy, + typename P, + typename PS +> +struct result_from_comparable + < + cross_track_geo, P, PS + > +{ +private : + typedef cross_track_geo strategy_type; +// typedef typename return_type::type return_type; +public : + template + static inline T apply(strategy_type const& strategy, + T const& comparable_distance) + { + //T c = T(2.0) * asin(math::sqrt(a)); + //return c * strategy.radius(); + return comparable_distance; + } +}; +*/ + +template +struct default_strategy + < + point_tag, segment_tag, Point, PointOfSegment, + geographic_tag, geographic_tag + > +{ + typedef cross_track_geo<> type; +}; + + +template +struct default_strategy + < + segment_tag, point_tag, PointOfSegment, Point, + geographic_tag, geographic_tag + > +{ + typedef typename default_strategy + < + point_tag, segment_tag, Point, PointOfSegment, + geographic_tag, geographic_tag + >::type type; +}; + +} // namespace services +#endif // DOXYGEN_NO_STRATEGY_SPECIALIZATIONS + +}} // namespace strategy::distance + +}} // namespace boost::geometry +#endif // BOOST_GEOMETRY_STRATEGIES_GEOGRAPHIC_DISTANCE_CROSS_TRACK_HPP diff --git a/include/boost/geometry/strategies/geographic/parameters.hpp b/include/boost/geometry/strategies/geographic/parameters.hpp index 5638db50f..c70c636f0 100644 --- a/include/boost/geometry/strategies/geographic/parameters.hpp +++ b/include/boost/geometry/strategies/geographic/parameters.hpp @@ -14,6 +14,8 @@ #include #include #include +#include +#include #include #include @@ -24,6 +26,23 @@ namespace boost { namespace geometry { namespace strategy struct andoyer { + //TODO: this should be replaced by an andoyer direct formula + template + < + typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct direct + : formula::thomas_direct + < + CT, EnableCoordinates, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; + template < typename CT, @@ -45,6 +64,22 @@ struct andoyer struct thomas { + template + < + typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct direct + : formula::thomas_direct + < + CT, EnableCoordinates, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; + template < typename CT, @@ -66,6 +101,22 @@ struct thomas struct vincenty { + template + < + typename CT, + bool EnableCoordinates = true, + bool EnableReverseAzimuth = false, + bool EnableReducedLength = false, + bool EnableGeodesicScale = false + > + struct direct + : formula::vincenty_direct + < + CT, EnableCoordinates, EnableReverseAzimuth, + EnableReducedLength, EnableGeodesicScale + > + {}; + template < typename CT, diff --git a/include/boost/geometry/strategies/strategies.hpp b/include/boost/geometry/strategies/strategies.hpp index 27a025c7c..05895ddf5 100644 --- a/include/boost/geometry/strategies/strategies.hpp +++ b/include/boost/geometry/strategies/strategies.hpp @@ -83,6 +83,7 @@ #include #include #include +#include #include #include #include diff --git a/test/algorithms/distance/Jamfile.v2 b/test/algorithms/distance/Jamfile.v2 index 656555a57..e19750a66 100644 --- a/test/algorithms/distance/Jamfile.v2 +++ b/test/algorithms/distance/Jamfile.v2 @@ -26,4 +26,5 @@ test-suite boost-geometry-algorithms-distance [ run distance_se_pl_l.cpp : : : : algorithms_distance_se_pl_l ] [ run distance_se_pl_pl.cpp : : : : algorithms_distance_se_pl_pl ] [ run distance_se_point_box.cpp : : : : algorithms_distance_se_point_box ] + [ run distance_geo_pl_l.cpp : : : : algorithms_distance_geo_pl_l ] ; diff --git a/test/algorithms/distance/distance_brute_force.hpp b/test/algorithms/distance/distance_brute_force.hpp index b62c09105..7992ffcae 100644 --- a/test/algorithms/distance/distance_brute_force.hpp +++ b/test/algorithms/distance/distance_brute_force.hpp @@ -170,7 +170,7 @@ struct distance_brute_force }; -template +template < typename Point1, typename Point2, @@ -184,7 +184,7 @@ struct distance_brute_force {}; -template +template < typename Point, typename Segment, @@ -198,7 +198,7 @@ struct distance_brute_force {}; -template +template < typename Point, typename Box, @@ -226,7 +226,7 @@ struct distance_brute_force {}; -template +template < typename Point, typename Linear, @@ -258,7 +258,7 @@ struct distance_brute_force }; -template +template < typename Point, typename MultiPoint, @@ -286,7 +286,7 @@ struct distance_brute_force } }; -template +template < typename MultiPoint1, typename MultiPoint2, @@ -320,7 +320,7 @@ struct distance_brute_force }; -template +template < typename MultiPoint, typename Linear, @@ -354,7 +354,7 @@ struct distance_brute_force }; -template +template < typename Linear, typename MultiPoint, @@ -384,7 +384,7 @@ struct distance_brute_force }; -template +template < typename MultiPoint, typename Segment, @@ -413,7 +413,7 @@ struct distance_brute_force }; -template +template < typename Linear, typename Segment, @@ -445,7 +445,7 @@ struct distance_brute_force }; -template +template < typename Linear1, typename Linear2, diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp new file mode 100644 index 000000000..f30665de8 --- /dev/null +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -0,0 +1,370 @@ + +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2016-2017, Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle + +// Licensed under the Boost Software License version 1.0. +// http://www.boost.org/users/license.html + +#include + +#ifndef BOOST_TEST_MODULE +#define BOOST_TEST_MODULE test_distance_geographic_pl_l +#endif + +#include +#include + +#include "test_distance_geo_common.hpp" + +typedef bg::cs::geographic cs_type; +typedef bg::model::point point_type; +typedef bg::model::segment segment_type; +typedef bg::model::multi_point multi_point_type; +typedef bg::model::segment segment_type; +typedef bg::model::linestring linestring_type; +typedef bg::model::multi_linestring multi_linestring_type; + +namespace services = bg::strategy::distance::services; +typedef bg::default_distance_result::type return_type; + +//typedef bg::strategy::distance::haversine point_point_strategy; +//typedef bg::strategy::distance::cross_track<> point_segment_strategy; + +typedef bg::srs::spheroid stype; + +// Strategies for point-point distance + +typedef bg::strategy::distance::andoyer andoyer_pp; +typedef bg::strategy::distance::thomas thomas_pp; +typedef bg::strategy::distance::vincenty vincenty_pp; + +// Strategies for point-segment distance + +typedef bg::strategy::distance::cross_track_geo + andoyer_strategy; + +typedef bg::strategy::distance::cross_track_geo + thomas_strategy; + +typedef bg::strategy::distance::cross_track_geo + vincenty_strategy; + +//=========================================================================== + +template +inline bg::default_distance_result::type +pp_distance(std::string const& wkt1, + std::string const& wkt2, + Strategy const& strategy) +{ + point_type p1, p2; + bg::read_wkt(wkt1, p1); + bg::read_wkt(wkt2, p2); + return bg::distance(p1, p2, strategy); +} + +template +inline bg::default_comparable_distance_result::type +pp_comparable_distance(std::string const& wkt1, + std::string const& wkt2, + Strategy const&) +{ + point_type p1, p2; + bg::read_wkt(wkt1, p1); + bg::read_wkt(wkt2, p2); + return bg::comparable_distance(p1, p2); +} + +template +inline bg::default_distance_result::type +ps_distance(std::string const& wkt1, + std::string const& wkt2, + Strategy const& strategy) +{ + point_type p; + segment_type s; + bg::read_wkt(wkt1, p); + bg::read_wkt(wkt2, s); + return bg::distance(p, s, strategy); +} + +template +inline bg::default_comparable_distance_result::type +ps_comparable_distance(std::string const& wkt1, + std::string const& wkt2, + Strategy const& strategy) +{ + point_type p; + segment_type s; + bg::read_wkt(wkt1, p); + bg::read_wkt(wkt2, s); + return bg::comparable_distance(p, s, strategy); +} + +template +T to_comparable(Strategy const& strategy, T const& distance) +{ + namespace services = bg::strategy::distance::services; + + typedef typename services::comparable_type + < + Strategy + >::type comparable_strategy; + + typedef typename services::result_from_distance + < + comparable_strategy, + point_type, + bg::point_type::type + > get_comparable_distance; + + comparable_strategy cstrategy = services::get_comparable + < + Strategy + >::apply(strategy); + + return get_comparable_distance::apply(cstrategy, distance); +} + +//=========================================================================== + +template +void test_distance_point_segment(Strategy_pp const& strategy_pp, + Strategy_ps const& strategy_ps) +{ + +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << std::endl; + std::cout << "point/segment distance tests" << std::endl; +#endif + typedef test_distance_of_geometries tester; + + tester::apply("p-s-01", + "POINT(0 0)", + "SEGMENT(2 0,3 0)", + pp_distance("POINT(0 0)", "POINT(2 0)", strategy_pp), + strategy_ps); + tester::apply("p-s-02", + "POINT(2.5 3)", + "SEGMENT(2 0,3 0)", + pp_distance("POINT(2.5 3)", "POINT(2.49777 0)", strategy_pp), + strategy_ps); + tester::apply("p-s-03", + "POINT(2 0)", + "SEGMENT(2 0,3 0)", + 0, + strategy_ps); + tester::apply("p-s-04", + "POINT(3 0)", + "SEGMENT(2 0,3 0)", + 0, + strategy_ps); + tester::apply("p-s-06", + "POINT(3.5 3)", + "SEGMENT(2 0,3 0)", + pp_distance("POINT(3 0)", "POINT(3.5 3)", strategy_pp), + strategy_ps); + tester::apply("p-s-07", + "POINT(15 80)", + "SEGMENT(10 15,30 15)", + 7204174.8264546748, + strategy_ps); + tester::apply("p-s-08", + "POINT(15 10)", + "SEGMENT(10 15,30 15)", + 571412.71247283253, + strategy_ps); + tester::apply("p-s-09", + "POINT(5 10)", + "SEGMENT(10 15,30 15)", + pp_distance("POINT(5 10)", "POINT(10 15)", strategy_pp), + strategy_ps); + tester::apply("p-s-10", + "POINT(35 10)", + "SEGMENT(10 15,30 15)", + pp_distance("POINT(35 10)", "POINT(30 15)", strategy_pp), + strategy_ps); + tester::apply("p-s-11", + "POINT(5 10)", + "SEGMENT(30 15,10 15)", + pp_distance("POINT(5 10)", "POINT(10 15)", strategy_pp), + strategy_ps); + tester::apply("p-s-12", + "POINT(35 10)", + "SEGMENT(30 15,10 15)", + pp_distance("POINT(35 10)", "POINT(30 15)", strategy_pp), + strategy_ps); + + tester::apply("p-s-right-up", + "POINT(3.5 3)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(3 2)", "POINT(3.5 3)", strategy_pp), + strategy_ps); + + tester::apply("p-s-left-up", + "POINT(1.5 3)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(2 2)", "POINT(1.5 3)", strategy_pp), + strategy_ps); + + tester::apply("p-s-up-1", + "POINT(2 3)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(2.0003 2)", "POINT(2 3)", strategy_pp), + strategy_ps); + + tester::apply("p-s-up-2", + "POINT(3 3)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(2.9997 2)", "POINT(3 3)", strategy_pp), + strategy_ps); + + tester::apply("p-s-right-down", + "POINT(3.5 1)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(3 2)", "POINT(3.5 1)", strategy_pp), + strategy_ps); + + tester::apply("p-s-left-down", + "POINT(1.5 1)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(2 2)", "POINT(1.5 1)", strategy_pp), + strategy_ps); + + tester::apply("p-s-down-1", + "POINT(2 1)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(2.0003 2)", "POINT(2 1)", strategy_pp), + strategy_ps); + + tester::apply("p-s-down-2", + "POINT(3 1)", + "SEGMENT(2 2,3 2)", + pp_distance("POINT(2.9997 2)", "POINT(3 1)", strategy_pp), + strategy_ps); + + tester::apply("p-s-south", + "POINT(3 -1)", + "SEGMENT(2 -2,3 -2)", + pp_distance("POINT(2.9997 -2)", "POINT(3 -1)", strategy_pp), + strategy_ps); + + tester::apply("p-s-antimeridian-1", + "POINT(3 1)", + "SEGMENT(220 2,3 2)", + pp_distance("POINT(3 2)", "POINT(3 1)", strategy_pp), + strategy_ps); + + tester::apply("p-s-antimeridian-2", + "POINT(220 1)", + "SEGMENT(220 2,3 2)", + pp_distance("POINT(220 2)", "POINT(220 1)", strategy_pp), + strategy_ps); + + // equator special case + tester::apply("p-s-eq1", + "POINT(2.5 0)", + "SEGMENT(2 0,3 0)", + 0, + strategy_ps); + tester::apply("p-s-eq2", + "POINT(2.5 2)", + "SEGMENT(2 0,3 0)", + pp_distance("POINT(2.5 0)", "POINT(2.5 2)", strategy_pp), + strategy_ps); + tester::apply("p-s-eq3", + "POINT(2 2)", + "SEGMENT(2 0,3 0)", + pp_distance("POINT(2 0)", "POINT(2 2)", strategy_pp), + strategy_ps); + tester::apply("p-s-eq4", + "POINT(3 2)", + "SEGMENT(2 0,3 0)", + pp_distance("POINT(3 0)", "POINT(3 2)", strategy_pp), + strategy_ps); + + // meridian special case + tester::apply("p-s-mer1", + "POINT(2.5 2)", + "SEGMENT(2 2,2 4)", + pp_distance("POINT(2.5 2)", "POINT(2 2)", strategy_pp), + strategy_ps); + tester::apply("p-s-mer2", + "POINT(2.5 3)", + "SEGMENT(2 2,2 4)", + pp_distance("POINT(2.5 3)", "POINT(2 3)", strategy_pp), + strategy_ps); + tester::apply("p-s-mer3", + "POINT(2.5 5)", + "SEGMENT(2 2,2 4)", + pp_distance("POINT(2.5 5)", "POINT(2 4)", strategy_pp), + strategy_ps); + + + //Nearly antipodal points not treated correctly in BG so far +/* + // very small distances to segment + tester::apply("p-s-07", + "POINT(90 1e-3)", + "SEGMENT(0.5 0,175.5 0)", + 110.57427582147321, + //1e-3 * d2r * 6335439.3272858784, + //to_comparable(strategy, 1e-3 * d2r * strategy.radius()), + strategy_ps); + tester::apply("p-s-08", + "POINT(90 1e-4)", + "SEGMENT(0.5 0,175.5 0)", + 6176.4469462110519, + //1e-4 * d2r * 6335439.3272858784, + //1e-4 * d2r * strategy.radius(), + //to_comparable(strategy, 1e-4 * d2r * strategy.radius()), + strategy_ps); + tester::apply("p-s-09", + "POINT(90 1e-5)", + "SEGMENT(0.5 0,175.5 0)", + 4917.5121192812849, + //1e-5 * d2r * strategy.radius(), + //to_comparable(strategy, 1e-5 * d2r * strategy.radius()), + strategy_ps); + tester::apply("p-s-10", + "POINT(90 1e-6)", + "SEGMENT(0.5 0,175.5 0)", + 4997.0354014473833, + //1e-6 * d2r * strategy.radius(), + //to_comparable(strategy, 1e-6 * d2r * strategy.radius()), + strategy_ps); + tester::apply("p-s-11", + "POINT(90 1e-7)", + "SEGMENT(0.5 0,175.5 0)", + 109854.20732250962, + //1e-7 * d2r * strategy.radius(), + //to_comparable(strategy, 1e-7 * d2r * strategy.radius()), + strategy_ps); + tester::apply("p-s-12", + "POINT(90 1e-8)", + "SEGMENT(0.5 0,175.5 0)", + 1372504.9365672076, + //1e-8 * d2r * strategy.radius(), + //to_comparable(strategy, 1e-8 * d2r * strategy.radius()), + strategy_ps); +*/ +} + +//=========================================================================== +//=========================================================================== +//=========================================================================== + +BOOST_AUTO_TEST_CASE( test_all_point_segment ) +{ + //test_distance_point_segment(andoyer_strategy()); + //test_distance_point_segment(thomas_strategy()); + test_distance_point_segment(vincenty_pp(), vincenty_strategy()); + test_distance_point_segment(thomas_pp(), thomas_strategy()); + test_distance_point_segment(andoyer_pp(), andoyer_strategy()); + //test_distance_point_segment(point_segment_strategy(earth_radius_km)); + //test_distance_point_segment(point_segment_strategy(earth_radius_miles)); +} diff --git a/test/algorithms/distance/test_distance_geo_common.hpp b/test/algorithms/distance/test_distance_geo_common.hpp new file mode 100644 index 000000000..5a820dc90 --- /dev/null +++ b/test/algorithms/distance/test_distance_geo_common.hpp @@ -0,0 +1,371 @@ + +// Boost.Geometry (aka GGL, Generic Geometry Library) +// Unit Test + +// Copyright (c) 2016-2017, Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle + +// Licensed under the Boost Software License version 1.0. +// http://www.boost.org/users/license.html + +#ifndef BOOST_GEOMETRY_TEST_DISTANCE_GEO_COMMON_HPP +#define BOOST_GEOMETRY_TEST_DISTANCE_GEO_COMMON_HPP + +#include +#include + +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + +#include +#include +#include + +#include + +#include +#include + +#include "distance_brute_force.hpp" + +namespace bg = ::boost::geometry; + +//======================================================================== + + +template +struct check_equal +{ + template + struct equal_to + { + static inline void apply(Value const& x, Value const& y) + { + BOOST_CHECK(x == y); + } + }; + + template + struct equal_to + { + static inline void apply(double x, double y) + { + BOOST_CHECK_CLOSE(x, y, 0.001); + } + }; + + template + static inline void apply(std::string const& /*case_id*/, + std::string const& /*subcase_id*/, + Geometry1 const& /*geometry1*/, + Geometry2 const& /*geometry2*/, + T const& detected, + T const& expected) + { + equal_to::apply(expected, detected); + /* + TODO: + Ideally we would want the following, but it does not work well + approximate equality test. + BOOST_CHECK_MESSAGE(equal_to::apply(expected, detected), + "case ID: " << case_id << "-" << subcase_id << "; " + << "G1: " << bg::wkt(geometry1) + << " - " + << "G2: " << bg::wkt(geometry2) + << " -> Detected: " << detected + << "; Expected: " << expected); + */ + } +}; + +//======================================================================== + +template +< + typename Geometry1, typename Geometry2, + int id1 = bg::geometry_id::value, + int id2 = bg::geometry_id::value +> +struct test_distance_of_geometries + : public test_distance_of_geometries +{}; + + +template +struct test_distance_of_geometries +{ + template + < + typename DistanceType, + typename ComparableDistanceType, + typename Strategy + > + static inline + void apply(std::string const& case_id, + std::string const& wkt1, + std::string const& wkt2, + DistanceType const& expected_distance, + ComparableDistanceType const& expected_comparable_distance, + Strategy const& strategy, + bool test_reversed = true) + { + Geometry1 geometry1 = from_wkt(wkt1); + Geometry2 geometry2 = from_wkt(wkt2); + + apply(case_id, geometry1, geometry2, + expected_distance, expected_comparable_distance, + strategy, test_reversed); + } + + template + static inline + void apply(std::string const& case_id, + std::string const& wkt1, + std::string const& wkt2, + DistanceType const& expected_distance, + Strategy const& strategy, + bool test_reversed = true) + { + Geometry1 geometry1 = from_wkt(wkt1); + Geometry2 geometry2 = from_wkt(wkt2); + + apply(case_id, geometry1, geometry2, + expected_distance, expected_distance, + strategy, test_reversed); + } + + + template + < + typename DistanceType, + typename ComparableDistanceType, + typename Strategy + > + static inline + void apply(std::string const& case_id, + Geometry1 const& geometry1, + Geometry2 const& geometry2, + DistanceType const& expected_distance, + ComparableDistanceType const& expected_comparable_distance, + Strategy const& strategy, + bool test_reversed = true) + { +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << "case ID: " << case_id << "; " + << "G1: " << bg::wkt(geometry1) + << " - " + << "G2: " << bg::wkt(geometry2) + << std::endl; +#endif + namespace services = bg::strategy::distance::services; + + using bg::unit_test::distance_brute_force; + + typedef typename bg::default_distance_result + < + Geometry1, Geometry2 + >::type default_distance_result; + + typedef typename services::return_type + < + Strategy, Geometry1, Geometry2 + >::type distance_result_from_strategy; + + static const bool same_regular = boost::is_same + < + default_distance_result, + distance_result_from_strategy + >::type::value; + + BOOST_CHECK(same_regular); +/* + typedef typename bg::default_comparable_distance_result + < + Geometry1, Geometry2 + >::type default_comparable_distance_result; + typedef typename services::return_type + < + typename services::comparable_type::type, + Geometry1, + Geometry2 + >::type comparable_distance_result_from_strategy; + static const bool same_comparable = boost::is_same + < + default_comparable_distance_result, + comparable_distance_result_from_strategy + >::type::value; + BOOST_CHECK( same_comparable ); +*/ + // check distance with passed strategy + distance_result_from_strategy dist = + bg::distance(geometry1, geometry2, strategy); + + check_equal + < + distance_result_from_strategy + >::apply(case_id, "a", geometry1, geometry2, + dist, expected_distance); + + // check against the comparable distance computed in a + // brute-force manner + default_distance_result dist_brute_force + = distance_brute_force(geometry1, geometry2, strategy); + + check_equal + < + default_distance_result + >::apply(case_id, "b", geometry1, geometry2, + dist_brute_force, expected_distance); + +/* + // check comparable distance with passed strategy + comparable_distance_result_from_strategy cdist = + bg::comparable_distance(geometry1, geometry2, strategy); + check_equal + < + default_comparable_distance_result + >::apply(case_id, "c", geometry1, geometry2, + cdist, expected_comparable_distance); + // check against the comparable distance computed in a + // brute-force manner + default_comparable_distance_result cdist_brute_force + = distance_brute_force(geometry1, + geometry2, + services::get_comparable + < + Strategy + >::apply(strategy)); + check_equal + < + default_comparable_distance_result + >::apply(case_id, "d", geometry1, geometry2, + cdist_brute_force, expected_comparable_distance); +*/ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << string_from_type::type>::name() + << string_from_type::type>::name() + << " -> " + << string_from_type::name() + //<< string_from_type::name() + << std::endl; + //std::cout << "strategy radius: " << strategy.radius() << std::endl; + std::cout << "expected distance = " + << expected_distance << " ; " + << "expected comp. distance = " + // << expected_comparable_distance + << std::endl; + std::cout << "distance = " + << dist << " ; " + << "comp. distance = " + //<< cdist + << std::endl; + + if ( !test_reversed ) + { + std::cout << std::endl; + } +#endif + + if ( test_reversed ) + { + // check distance with given strategy + dist = bg::distance(geometry2, geometry1, strategy); + + check_equal + < + default_distance_result + >::apply(case_id, "ra", geometry2, geometry1, + dist, expected_distance); +/* + // check comparable distance with given strategy + cdist = bg::comparable_distance(geometry2, geometry1, strategy); + check_equal + < + default_comparable_distance_result + >::apply(case_id, "rc", geometry2, geometry1, + cdist, expected_comparable_distance); +*/ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << "distance[reversed args] = " + << dist << " ; " + << "comp. distance[reversed args] = " + //<< cdist + << std::endl; + std::cout << std::endl; +#endif + } + + } +}; + + +//======================================================================== + + +template +void test_empty_input(Geometry1 const& geometry1, + Geometry2 const& geometry2, + Strategy const& strategy) +{ + try + { + bg::distance(geometry1, geometry2); + } + catch(bg::empty_input_exception const& ) + { + return; + } + BOOST_CHECK_MESSAGE(false, + "A empty_input_exception should have been thrown"); + + try + { + bg::distance(geometry2, geometry1); + } + catch(bg::empty_input_exception const& ) + { + return; + } + BOOST_CHECK_MESSAGE(false, + "A empty_input_exception should have been thrown"); + + try + { + bg::distance(geometry1, geometry2, strategy); + } + catch(bg::empty_input_exception const& ) + { + return; + } + BOOST_CHECK_MESSAGE(false, + "A empty_input_exception should have been thrown"); + + try + { + bg::distance(geometry2, geometry1, strategy); + } + catch(bg::empty_input_exception const& ) + { + return; + } + BOOST_CHECK_MESSAGE(false, + "A empty_input_exception should have been thrown"); +} + +#endif // BOOST_GEOMETRY_TEST_DISTANCE_GEO_COMMON_HPP From 36232a3cdb4e040baf15fbe10ed44054852312ee Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Thu, 13 Jul 2017 12:52:52 +0300 Subject: [PATCH 02/27] [distance] Spheroid passed as parameter, more general strategy getter for distance strategy --- .../formulas/distance_point_segment.hpp | 13 ++- .../geographic/distance_cross_track.hpp | 79 ++++++------------- .../strategies/geographic/intersection.hpp | 6 +- .../algorithms/distance/distance_geo_pl_l.cpp | 6 +- .../distance/test_distance_geo_common.hpp | 1 - 5 files changed, 35 insertions(+), 70 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 6363aac8b..a6ce5c77f 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -25,26 +25,25 @@ template typename Inverse_type_distance, typename Direct_type > -class cross_track_geo_formula{ +class distance_point_segment{ public: + template CT static inline apply(CT lon1, CT lat1, //p1 CT lon2, CT lat2, //p2 - CT lon3, CT lat3) //query point p3 + CT lon3, CT lat3, //query point p3 + Spheroid const& spheroid, + CT earth_radius = + geometry::srs::sphere().get_radius<1>()) { int print = 0; // Constants - geometry::srs::spheroid spheroid; CT const f = flattening(spheroid); CT const pi = math::pi(); CT const half_pi = math::pi() / CT(2); - geometry::srs::sphere sphere; - CT const earth_radius = sphere.get_radius<1>(); - //CT const earth_radius = 6372795; - //CT const earth_radius = 6378137; // Convert to radians lon1 = math::as_radian(lon1); diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp index b339f8f45..64009ee6a 100644 --- a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -49,18 +49,15 @@ namespace boost { namespace geometry namespace strategy { namespace distance { - /*! -\brief Strategy functor for distance point to segment calculation +\brief Strategy functor for distance point to segment calculation on ellipsoid \ingroup strategies -\details Class which calculates the distance of a point to a segment, for points on a sphere or globe -\see http://williams.best.vwh.net/avform.htm +\details Class which calculates the distance of a point to a segment, for points +on the ellipsoid +\see https://arxiv.org/abs/1102.1215 +\tparam FormulaPolicy underlying point-point distance strategy +\tparam Spheroid is the spheroidal model used \tparam CalculationType \tparam_calculation -\tparam Strategy underlying point-point distance strategy -\qbk{ -[heading See also] -[link geometry.reference.algorithms.distance.distance_3_with_strategy distance (with strategy)] -} */ template < @@ -84,21 +81,21 @@ public : > {}; - inline cross_track_geo() + struct distance_strategy + { + typedef geographic type; + }; + + inline typename distance_strategy::type get_distance_strategy() const + { + typedef typename distance_strategy::type distance_type; + return distance_type(m_spheroid); + } + + explicit cross_track_geo(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) {} - typedef geographic - < - FormulaPolicy, Spheroid, CalculationType - > DistanceStrategy; - - /* - template - inline T comparable_to_distance(T& a) const - { - return a; - } -*/ template inline typename return_type::type apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const @@ -136,7 +133,7 @@ public : false > direct_type; - return geometry::formula::cross_track_geo_formula + return geometry::formula::distance_point_segment < CT, units_type, @@ -145,14 +142,12 @@ public : direct_type >::apply(get<0>(sp1), get<1>(sp1), get<0>(sp2), get<1>(sp2), - get<0>(p), get<1>(p)); + get<0>(p), get<1>(p), + m_spheroid); } - //TODO: apply a more general strategy getter - inline DistanceStrategy get_distance_strategy() const - { - return DistanceStrategy(); - } +private : + Spheroid m_spheroid; }; @@ -278,32 +273,6 @@ public : } }; -/* -template -< - typename FormulaPolicy, - typename P, - typename PS -> -struct result_from_comparable - < - cross_track_geo, P, PS - > -{ -private : - typedef cross_track_geo strategy_type; -// typedef typename return_type::type return_type; -public : - template - static inline T apply(strategy_type const& strategy, - T const& comparable_distance) - { - //T c = T(2.0) * asin(math::sqrt(a)); - //return c * strategy.radius(); - return comparable_distance; - } -}; -*/ template struct default_strategy diff --git a/include/boost/geometry/strategies/geographic/intersection.hpp b/include/boost/geometry/strategies/geographic/intersection.hpp index 84acd149c..f4c2c1436 100644 --- a/include/boost/geometry/strategies/geographic/intersection.hpp +++ b/include/boost/geometry/strategies/geographic/intersection.hpp @@ -118,7 +118,6 @@ struct geographic_segments return strategy_type(m_spheroid); } - template struct distance_strategy { typedef distance::geographic @@ -129,10 +128,9 @@ struct geographic_segments > type; }; - template - inline typename distance_strategy::type get_distance_strategy() const + inline typename distance_strategy::type get_distance_strategy() const { - typedef typename distance_strategy::type strategy_type; + typedef typename distance_strategy::type strategy_type; return strategy_type(m_spheroid); } diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index f30665de8..83858f64d 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -44,13 +44,13 @@ typedef bg::strategy::distance::vincenty vincenty_pp; // Strategies for point-segment distance -typedef bg::strategy::distance::cross_track_geo +typedef bg::strategy::distance::cross_track_geo andoyer_strategy; -typedef bg::strategy::distance::cross_track_geo +typedef bg::strategy::distance::cross_track_geo thomas_strategy; -typedef bg::strategy::distance::cross_track_geo +typedef bg::strategy::distance::cross_track_geo vincenty_strategy; //=========================================================================== diff --git a/test/algorithms/distance/test_distance_geo_common.hpp b/test/algorithms/distance/test_distance_geo_common.hpp index 5a820dc90..1f831d447 100644 --- a/test/algorithms/distance/test_distance_geo_common.hpp +++ b/test/algorithms/distance/test_distance_geo_common.hpp @@ -1,4 +1,3 @@ - // Boost.Geometry (aka GGL, Generic Geometry Library) // Unit Test From b8453a3caf9ab12e381297f941a1a3180bb00fa6 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 18 Jul 2017 13:05:15 +0300 Subject: [PATCH 03/27] [distance] [test] Clean code in distance test and point segment formula --- .../formulas/distance_point_segment.hpp | 200 +++++++++++------- .../boost/geometry/strategies/distance.hpp | 3 + .../geographic/distance_cross_track.hpp | 4 +- .../algorithms/distance/distance_geo_pl_l.cpp | 58 ----- .../distance/test_distance_geo_common.hpp | 99 +-------- 5 files changed, 128 insertions(+), 236 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index a6ce5c77f..a1b60993a 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -11,9 +11,12 @@ #ifndef BOOST_GEOMETRY_FORMULAS_CROSS_TRACK_GEO_HPP #define BOOST_GEOMETRY_FORMULAS_CROSS_TRACK_GEO_HPP -#include #include +#ifndef BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS +#define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100 +#endif + namespace boost { namespace geometry { namespace formula { @@ -23,27 +26,59 @@ template typename Units, typename Inverse_type_azimuth, typename Inverse_type_distance, - typename Direct_type + typename Direct_type, + bool EnableClosestPoint = false > class distance_point_segment{ public: + struct result_distance_point_segment + { + result_distance_point_segment() + : distance(0) + , closest_point_lon(0) + , closest_point_lat(0) + {} + + CT distance; + CT closest_point_lon; + CT closest_point_lat; + }; + template - CT static inline apply(CT lon1, CT lat1, //p1 - CT lon2, CT lat2, //p2 - CT lon3, CT lat3, //query point p3 - Spheroid const& spheroid, - CT earth_radius = + result_distance_point_segment + static inline non_iterative_case(CT lon1, CT lat1, //p1 + CT lon2, CT lat2, //p2 + Spheroid const& spheroid) + { + result_distance_point_segment result; + result.distance = Inverse_type_distance::apply(lon1, lat1, + lon2, lat2, + spheroid).distance; + if (EnableClosestPoint) + { + result.closest_point_lon = lon1; + result.closest_point_lat = lat1; + } + return result; + } + + template + result_distance_point_segment + static inline apply(CT lon1, CT lat1, //p1 + CT lon2, CT lat2, //p2 + CT lon3, CT lat3, //query point p3 + Spheroid const& spheroid, + CT earth_radius = geometry::srs::sphere().get_radius<1>()) { - - int print = 0; + result_distance_point_segment result; // Constants CT const f = flattening(spheroid); CT const pi = math::pi(); - CT const half_pi = math::pi() / CT(2); + CT const half_pi = pi / CT(2); // Convert to radians lon1 = math::as_radian(lon1); @@ -64,13 +99,13 @@ public: { if (lon3 <= lon1) { - return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } if (lon3 >= lon2) { - return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); } - return Inverse_type_distance::apply(lon3, lat1, lon3, lat3, spheroid).distance; + return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); } //segment on meridian @@ -82,23 +117,23 @@ public: } if (lat3 <= lat1) { - return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } if (lat3 >= lat2) { - return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); } - return Inverse_type_distance::apply(lon1, lat3, lon3, lat3, spheroid).distance; + return non_iterative_case(lon1, lat3, lon3, lat3, spheroid); } // Easy cases if (lon3 <= lon1 && lat3 <= lat1) { - return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } if (lon3 >= lon2 && lat3 <= lat2) { - return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); } // Compute a1 (GEO) @@ -107,17 +142,11 @@ public: CT a312 = a1 - a13; - if (print) - { - std::cout << "a1=" << a1 * math::r2d() << std::endl; - std::cout << "a13=" << a13 * math::r2d() << std::endl; - std::cout << "a312=" << a312 * math::r2d() << std::endl; - } if (a312 > half_pi) { // projection of p3 on geodesic spanned by segment (p1,p2) fall // outside of segment on the side of p1 - return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } CT a2 = pi + Inverse_type_azimuth::apply(lon2, lat2, lon1, lat1, spheroid).azimuth; @@ -125,17 +154,25 @@ public: CT a321 = a2 - a23; - if (print) - { - std::cout << "a2=" << a2 * math::r2d() << std::endl; - std::cout << "a23=" << a13 * math::r2d() << std::endl; - std::cout << "a321=" << a321 * math::r2d() << std::endl; - } +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "segment=(" << lon1 * math::r2d(); + std::cout << "," << lat1 * math::r2d(); + std::cout << "),(" << lon2 * math::r2d(); + std::cout << "," << lat2 * math::r2d(); + std::cout << ")\np=(" << lon3 * math::r2d(); + std::cout << "," << lat3 * math::r2d(); + std::cout << ")\na1=" << a1 * math::r2d() << std::endl; + std::cout << "a13=" << a13 * math::r2d() << std::endl; + std::cout << "a312=" << a312 * math::r2d() << std::endl; + std::cout << "a2=" << a2 * math::r2d() << std::endl; + std::cout << "a23=" << a13 * math::r2d() << std::endl; + std::cout << "a321=" << a321 * math::r2d() << std::endl; +#endif if (a321 < half_pi) { // projection of p3 on geodesic spanned by segment (p1,p2) fall // outside of segment on the side of p2 - return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); } // Guess s14 (SPHERICAL) @@ -153,35 +190,27 @@ public: geometry::strategy::distance::haversine str(earth_radius); CT s13 = str.apply(p1, p3); - /* - // Closest point is one of the segment's endpoints - if (s34 == s13) - { - return Inverse_type_distance::apply(lon1, lat1, lon3, lat3, spheroid).distance; - } - CT s23 = str.apply(p2, p3); - if (s34 == s23) - { - return Inverse_type_distance::apply(lon2, lat2, lon3, lat3, spheroid).distance; - } -*/ - CT s14 = acos( cos(s13/earth_radius) / cos(s34/earth_radius)) * earth_radius; + CT s14 = acos( cos(s13/earth_radius) / cos(s34/earth_radius) ) * earth_radius; - if (print) - { - std::cout << "s34=" << s34 << std::endl; - std::cout << "s13=" << s13 << std::endl; - std::cout << "s14=" << s14 << std::endl; - std::cout << "===============" << std::endl; - } +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + 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 = 0; geometry::formula::result_direct res14; geometry::formula::result_inverse res34; - for (int i=0; i<20; ++i) - { + int counter = 0; // robustness + CT g4; + CT delta_g4; + + do{ + prev_distance = res34.distance; + // Solve the direct problem to find p4 (GEO) res14 = Direct_type::apply(lon1, lat1, s14, a1, spheroid); @@ -190,11 +219,9 @@ public: CT a4 = Inverse_type_azimuth::apply(res14.lon2, res14.lat2, lon2, lat2, spheroid).azimuth; res34 = Inverse_type_distance::apply(res14.lon2, res14.lat2, lon3, lat3, spheroid); - CT g4 = res34.azimuth - a4; + g4 = res34.azimuth - a4; // Normalize g4 - CT delta_g4; - if (g4 < 0 && g4 < -pi)//close to -270 { delta_g4 = g4 + 1.5 * pi; @@ -217,31 +244,44 @@ public: CT der = M43 / m34; s14 = s14 - delta_g4 / der; - //Print - if (print) +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "p4=" << res14.lon2 * math::r2d() << + "," << res14.lat2 * math::r2d() << std::endl; + std::cout << "delta_g4=" << delta_g4 << std::endl; + std::cout << "g4=" << g4 * math::r2d() << std::endl; + std::cout << "der=" << der << 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; +#endif + result.distance = prev_distance; + +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + if (g4 == half_pi) { - std::cout << "p4=" << res14.lon2 * math::r2d() << - "," << res14.lat2 * math::r2d() << std::endl; - std::cout << "delta_g4=" << delta_g4 << std::endl; - std::cout << "g4=" << g4 * math::r2d() << std::endl; - std::cout << "der=" << der << std::endl; - std::cout << "new_s14=" << s14 << std::endl; - std::cout << "dist =" << res34.distance << std::endl; - std::cout << "---------end of step " << i << std::endl<< std::endl; + std::cout << "Stop msg: g4 == half_pi" << std::endl; } - - if (g4 == math::pi()/2 || (res34.distance >= prev_distance && prev_distance!=0) || delta_g4 == 0) + if (res34.distance >= prev_distance && prev_distance != 0) { - //std::cout << "stop" << std::endl; - return prev_distance; + std::cout << "Stop msg: res34.distance >= prev_distance" << std::endl; } - prev_distance = res34.distance; + if (delta_g4 == 0) + { + std::cout << "Stop msg: delta_g4 == 0" << std::endl; + } + if (counter == 19) + { + std::cout << "Stop msg: counter" << std::endl; + } +#endif - } + } while (g4 != half_pi + && (prev_distance > res34.distance || prev_distance == 0) + && delta_g4 != 0 + && ++counter < BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS ) ; - ////////////////////////////TESTING - - std::cout << "dist =" << res34.distance << std::endl; +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "distance=" << res34.distance << std::endl; point p4(res14.lon2, res14.lat2); CT s34_sph = str.apply(p4, p3); @@ -270,15 +310,15 @@ public: if (res34.distance <= p4_plus && res34.distance <= p4_minus) { - std::cout << "Closest point computed" << std::endl << std::endl; + std::cout << "Closest point computed" << std::endl; } else { - std::cout << "There is a closer point nearby!" << std::endl; + std::cout << "There is a closer point nearby" << std::endl; } +#endif - CT min_dist_endpoints = std::min(s31,s32); - return std::min(res34.distance, min_dist_endpoints); + return result; } }; diff --git a/include/boost/geometry/strategies/distance.hpp b/include/boost/geometry/strategies/distance.hpp index 0d5a9e6dd..98ccb8202 100644 --- a/include/boost/geometry/strategies/distance.hpp +++ b/include/boost/geometry/strategies/distance.hpp @@ -65,6 +65,9 @@ template struct get_comparable template struct result_from_distance {}; + + + // Default strategy diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp index 64009ee6a..60ebd54fd 100644 --- a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -133,7 +133,7 @@ public : false > direct_type; - return geometry::formula::distance_point_segment + return (geometry::formula::distance_point_segment < CT, units_type, @@ -143,7 +143,7 @@ public : >::apply(get<0>(sp1), get<1>(sp1), get<0>(sp2), get<1>(sp2), get<0>(p), get<1>(p), - m_spheroid); + m_spheroid)).distance; } private : diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 83858f64d..58cdb3b8b 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -1,4 +1,3 @@ - // Boost.Geometry (aka GGL, Generic Geometry Library) // Unit Test @@ -31,9 +30,6 @@ typedef bg::model::multi_linestring multi_linestring_type; namespace services = bg::strategy::distance::services; typedef bg::default_distance_result::type return_type; -//typedef bg::strategy::distance::haversine point_point_strategy; -//typedef bg::strategy::distance::cross_track<> point_segment_strategy; - typedef bg::srs::spheroid stype; // Strategies for point-point distance @@ -67,18 +63,6 @@ pp_distance(std::string const& wkt1, return bg::distance(p1, p2, strategy); } -template -inline bg::default_comparable_distance_result::type -pp_comparable_distance(std::string const& wkt1, - std::string const& wkt2, - Strategy const&) -{ - point_type p1, p2; - bg::read_wkt(wkt1, p1); - bg::read_wkt(wkt2, p2); - return bg::comparable_distance(p1, p2); -} - template inline bg::default_distance_result::type ps_distance(std::string const& wkt1, @@ -92,44 +76,6 @@ ps_distance(std::string const& wkt1, return bg::distance(p, s, strategy); } -template -inline bg::default_comparable_distance_result::type -ps_comparable_distance(std::string const& wkt1, - std::string const& wkt2, - Strategy const& strategy) -{ - point_type p; - segment_type s; - bg::read_wkt(wkt1, p); - bg::read_wkt(wkt2, s); - return bg::comparable_distance(p, s, strategy); -} - -template -T to_comparable(Strategy const& strategy, T const& distance) -{ - namespace services = bg::strategy::distance::services; - - typedef typename services::comparable_type - < - Strategy - >::type comparable_strategy; - - typedef typename services::result_from_distance - < - comparable_strategy, - point_type, - bg::point_type::type - > get_comparable_distance; - - comparable_strategy cstrategy = services::get_comparable - < - Strategy - >::apply(strategy); - - return get_comparable_distance::apply(cstrategy, distance); -} - //=========================================================================== template @@ -360,11 +306,7 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, BOOST_AUTO_TEST_CASE( test_all_point_segment ) { - //test_distance_point_segment(andoyer_strategy()); - //test_distance_point_segment(thomas_strategy()); test_distance_point_segment(vincenty_pp(), vincenty_strategy()); test_distance_point_segment(thomas_pp(), thomas_strategy()); test_distance_point_segment(andoyer_pp(), andoyer_strategy()); - //test_distance_point_segment(point_segment_strategy(earth_radius_km)); - //test_distance_point_segment(point_segment_strategy(earth_radius_miles)); } diff --git a/test/algorithms/distance/test_distance_geo_common.hpp b/test/algorithms/distance/test_distance_geo_common.hpp index 1f831d447..f1a06e0a2 100644 --- a/test/algorithms/distance/test_distance_geo_common.hpp +++ b/test/algorithms/distance/test_distance_geo_common.hpp @@ -34,7 +34,6 @@ #include #include -#include #include @@ -78,18 +77,6 @@ struct check_equal T const& expected) { equal_to::apply(expected, detected); - /* - TODO: - Ideally we would want the following, but it does not work well - approximate equality test. - BOOST_CHECK_MESSAGE(equal_to::apply(expected, detected), - "case ID: " << case_id << "-" << subcase_id << "; " - << "G1: " << bg::wkt(geometry1) - << " - " - << "G2: " << bg::wkt(geometry2) - << " -> Detected: " << detected - << "; Expected: " << expected); - */ } }; @@ -109,29 +96,6 @@ struct test_distance_of_geometries template struct test_distance_of_geometries { - template - < - typename DistanceType, - typename ComparableDistanceType, - typename Strategy - > - static inline - void apply(std::string const& case_id, - std::string const& wkt1, - std::string const& wkt2, - DistanceType const& expected_distance, - ComparableDistanceType const& expected_comparable_distance, - Strategy const& strategy, - bool test_reversed = true) - { - Geometry1 geometry1 = from_wkt(wkt1); - Geometry2 geometry2 = from_wkt(wkt2); - - apply(case_id, geometry1, geometry2, - expected_distance, expected_comparable_distance, - strategy, test_reversed); - } - template static inline void apply(std::string const& case_id, @@ -145,7 +109,7 @@ struct test_distance_of_geometries Geometry2 geometry2 = from_wkt(wkt2); apply(case_id, geometry1, geometry2, - expected_distance, expected_distance, + expected_distance, strategy, test_reversed); } @@ -153,7 +117,6 @@ struct test_distance_of_geometries template < typename DistanceType, - typename ComparableDistanceType, typename Strategy > static inline @@ -161,7 +124,6 @@ struct test_distance_of_geometries Geometry1 const& geometry1, Geometry2 const& geometry2, DistanceType const& expected_distance, - ComparableDistanceType const& expected_comparable_distance, Strategy const& strategy, bool test_reversed = true) { @@ -193,24 +155,7 @@ struct test_distance_of_geometries >::type::value; BOOST_CHECK(same_regular); -/* - typedef typename bg::default_comparable_distance_result - < - Geometry1, Geometry2 - >::type default_comparable_distance_result; - typedef typename services::return_type - < - typename services::comparable_type::type, - Geometry1, - Geometry2 - >::type comparable_distance_result_from_strategy; - static const bool same_comparable = boost::is_same - < - default_comparable_distance_result, - comparable_distance_result_from_strategy - >::type::value; - BOOST_CHECK( same_comparable ); -*/ + // check distance with passed strategy distance_result_from_strategy dist = bg::distance(geometry1, geometry2, strategy); @@ -232,47 +177,17 @@ struct test_distance_of_geometries >::apply(case_id, "b", geometry1, geometry2, dist_brute_force, expected_distance); -/* - // check comparable distance with passed strategy - comparable_distance_result_from_strategy cdist = - bg::comparable_distance(geometry1, geometry2, strategy); - check_equal - < - default_comparable_distance_result - >::apply(case_id, "c", geometry1, geometry2, - cdist, expected_comparable_distance); - // check against the comparable distance computed in a - // brute-force manner - default_comparable_distance_result cdist_brute_force - = distance_brute_force(geometry1, - geometry2, - services::get_comparable - < - Strategy - >::apply(strategy)); - check_equal - < - default_comparable_distance_result - >::apply(case_id, "d", geometry1, geometry2, - cdist_brute_force, expected_comparable_distance); -*/ #ifdef BOOST_GEOMETRY_TEST_DEBUG std::cout << string_from_type::type>::name() << string_from_type::type>::name() << " -> " << string_from_type::name() - //<< string_from_type::name() << std::endl; - //std::cout << "strategy radius: " << strategy.radius() << std::endl; std::cout << "expected distance = " << expected_distance << " ; " - << "expected comp. distance = " - // << expected_comparable_distance << std::endl; std::cout << "distance = " << dist << " ; " - << "comp. distance = " - //<< cdist << std::endl; if ( !test_reversed ) @@ -291,15 +206,7 @@ struct test_distance_of_geometries default_distance_result >::apply(case_id, "ra", geometry2, geometry1, dist, expected_distance); -/* - // check comparable distance with given strategy - cdist = bg::comparable_distance(geometry2, geometry1, strategy); - check_equal - < - default_comparable_distance_result - >::apply(case_id, "rc", geometry2, geometry1, - cdist, expected_comparable_distance); -*/ + #ifdef BOOST_GEOMETRY_TEST_DEBUG std::cout << "distance[reversed args] = " << dist << " ; " From ca2ab199cf196fd5c0ea0568a59a148a5011c063 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 18 Jul 2017 14:00:10 +0300 Subject: [PATCH 04/27] [formula] Documentation for distance point segment --- .../boost/geometry/formulas/distance_point_segment.hpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index a1b60993a..83f443ab2 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -17,6 +17,14 @@ #define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100 #endif +/*! +\brief Algorithm to compute the distance between a segment and a point using + direct and inverse geodesic problems as subroutines. The algorithm + approximates the distance by an iterative Newton method. +\see C.F.F.Karney - Geodesics on an ellipsoid of revolution, + https://arxiv.org/abs/1102.1215 +*/ + namespace boost { namespace geometry { namespace formula { From dbe6df7ad3531f944a4e1904d4eee3e7a6ebdaa8 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 30 Aug 2017 15:04:01 +0300 Subject: [PATCH 05/27] [formulas] [strategies] [distance] Change interface and earth radius in distance point segment geographic formula --- .../formulas/distance_point_segment.hpp | 52 +++++++------ .../geographic/distance_cross_track.hpp | 74 ++++++------------- .../strategies/geographic/intersection.hpp | 6 +- .../algorithms/distance/distance_geo_pl_l.cpp | 6 +- 4 files changed, 59 insertions(+), 79 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 83f443ab2..a2f32cd2a 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -25,6 +25,8 @@ https://arxiv.org/abs/1102.1215 */ +#define BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + namespace boost { namespace geometry { namespace formula { @@ -32,15 +34,18 @@ template < typename CT, typename Units, - typename Inverse_type_azimuth, - typename Inverse_type_distance, - typename Direct_type, + template class Inverse, + template class Direct, bool EnableClosestPoint = false > class distance_point_segment{ public: + typedef Inverse inverse_distance_type; + typedef Inverse inverse_azimuth_type; + typedef Direct direct_distance_type; + struct result_distance_point_segment { result_distance_point_segment() @@ -61,7 +66,7 @@ public: Spheroid const& spheroid) { result_distance_point_segment result; - result.distance = Inverse_type_distance::apply(lon1, lat1, + result.distance = inverse_distance_type::apply(lon1, lat1, lon2, lat2, spheroid).distance; if (EnableClosestPoint) @@ -77,10 +82,11 @@ public: static inline apply(CT lon1, CT lat1, //p1 CT lon2, CT lat2, //p2 CT lon3, CT lat3, //query point p3 - Spheroid const& spheroid, - CT earth_radius = - geometry::srs::sphere().get_radius<1>()) + Spheroid const& spheroid) { + CT earth_radius = (CT(2) * spheroid.get_radius<1>() + + spheroid.get_radius<2>()) / CT(3); + result_distance_point_segment result; // Constants @@ -145,8 +151,8 @@ public: } // Compute a1 (GEO) - CT a1 = Inverse_type_azimuth::apply(lon1, lat1, lon2, lat2, spheroid).azimuth; - CT a13 = Inverse_type_azimuth::apply(lon1, lat1, lon3, lat3, spheroid).azimuth; + CT a1 = inverse_azimuth_type::apply(lon1, lat1, lon2, lat2, spheroid).azimuth; + CT a13 = inverse_azimuth_type::apply(lon1, lat1, lon3, lat3, spheroid).azimuth; CT a312 = a1 - a13; @@ -157,8 +163,8 @@ public: return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } - CT a2 = pi + Inverse_type_azimuth::apply(lon2, lat2, lon1, lat1, spheroid).azimuth; - CT a23 = Inverse_type_azimuth::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; + CT a2 = pi + inverse_azimuth_type::apply(lon2, lat2, lon1, lat1, spheroid).azimuth; + CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; CT a321 = a2 - a23; @@ -220,13 +226,15 @@ public: prev_distance = res34.distance; // Solve the direct problem to find p4 (GEO) - res14 = Direct_type::apply(lon1, lat1, s14, a1, spheroid); + res14 = direct_distance_type::apply(lon1, lat1, s14, a1, 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_type_azimuth::apply(res14.lon2, res14.lat2, lon2, lat2, spheroid).azimuth; - res34 = Inverse_type_distance::apply(res14.lon2, res14.lat2, lon3, lat3, spheroid); + CT a4 = inverse_azimuth_type::apply(res14.lon2, res14.lat2, + lon2, lat2, spheroid).azimuth; + res34 = inverse_distance_type::apply(res14.lon2, res14.lat2, + lon3, lat3, spheroid); g4 = res34.azimuth - a4; // Normalize g4 @@ -296,20 +304,20 @@ public: std::cout << "s34(sph) =" << s34_sph << std::endl; std::cout << "s34(geo) =" - << Inverse_type_distance::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance + << inverse_distance_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance << ", p4=(" << get<0>(p4) * math::r2d() << "," << get<1>(p4) * math::r2d() << ")" << std::endl; - CT s31 = Inverse_type_distance::apply(lon3, lat3, lon1, lat1, spheroid).distance; - CT s32 = Inverse_type_distance::apply(lon3, lat3, lon2, lat2, spheroid).distance; + CT s31 = inverse_distance_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; + CT s32 = inverse_distance_type::apply(lon3, lat3, lon2, lat2, spheroid).distance; - CT a4 = Inverse_type_azimuth::apply(get<0>(p4), get<1>(p4), lon2, lat2, spheroid).azimuth; - geometry::formula::result_direct res4 = Direct_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); - CT p4_plus = Inverse_type_distance::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; + CT a4 = inverse_azimuth_type::apply(get<0>(p4), get<1>(p4), lon2, lat2, spheroid).azimuth; + geometry::formula::result_direct res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); + CT p4_plus = inverse_distance_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; - geometry::formula::result_direct res1 = Direct_type::apply(lon1, lat1, s14-.04, a1, spheroid); - CT p4_minus = Inverse_type_distance::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; + geometry::formula::result_direct res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a1, spheroid); + CT p4_minus = inverse_distance_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() << "," << res4.lat2 * math::r2d() << ")" diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp index 60ebd54fd..295858f6e 100644 --- a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -65,7 +65,7 @@ template typename Spheroid = srs::spheroid, typename CalculationType = void > -class cross_track_geo +class geographic_cross_track { public : template @@ -92,7 +92,7 @@ public : return distance_type(m_spheroid); } - explicit cross_track_geo(Spheroid const& spheroid = Spheroid()) + explicit geographic_cross_track(Spheroid const& spheroid = Spheroid()) : m_spheroid(spheroid) {} @@ -104,42 +104,12 @@ public : typedef typename return_type::type CT; - typedef typename FormulaPolicy::template inverse - < - CT, - false, - true, - false, - false, - false - > inverse_type_azimuth; - - typedef typename FormulaPolicy::template inverse - < - CT, - true, - true, - false, - true, - true - > inverse_type_distance; - - typedef typename FormulaPolicy::template direct - < - CT, - true, - false, - false, - false - > direct_type; - return (geometry::formula::distance_point_segment < CT, units_type, - inverse_type_azimuth, - inverse_type_distance, - direct_type + FormulaPolicy::template inverse, + FormulaPolicy::template direct >::apply(get<0>(sp1), get<1>(sp1), get<0>(sp2), get<1>(sp2), get<0>(p), get<1>(p), @@ -159,7 +129,7 @@ namespace services //tags template -struct tag > +struct tag > { typedef strategy_tag_distance_point_segment type; }; @@ -169,7 +139,7 @@ template typename FormulaPolicy, typename Spheroid > -struct tag > +struct tag > { typedef strategy_tag_distance_point_segment type; }; @@ -180,7 +150,7 @@ template typename Spheroid, typename CalculationType > -struct tag > +struct tag > { typedef strategy_tag_distance_point_segment type; }; @@ -188,8 +158,8 @@ struct tag > //return types template -struct return_type, P, PS> - : cross_track_geo::template return_type +struct return_type, P, PS> + : geographic_cross_track::template return_type {}; template @@ -199,8 +169,8 @@ template typename P, typename PS > -struct return_type, P, PS> - : cross_track_geo::template return_type +struct return_type, P, PS> + : geographic_cross_track::template return_type {}; template @@ -211,8 +181,8 @@ template typename P, typename PS > -struct return_type, P, PS> - : cross_track_geo::template return_type +struct return_type, P, PS> + : geographic_cross_track::template return_type {}; //comparable types @@ -222,9 +192,9 @@ template typename Spheroid, typename CalculationType > -struct comparable_type > +struct comparable_type > { - typedef cross_track_geo + typedef geographic_cross_track < FormulaPolicy, Spheroid, CalculationType > type; @@ -236,15 +206,15 @@ template typename Spheroid, typename CalculationType > -struct get_comparable > +struct get_comparable > { typedef typename comparable_type < - cross_track_geo + geographic_cross_track >::type comparable_type; public : static inline comparable_type - apply(cross_track_geo const& strategy) + apply(geographic_cross_track const& strategy) { return comparable_type(); } @@ -257,17 +227,17 @@ template typename P, typename PS > -struct result_from_distance, P, PS> +struct result_from_distance, P, PS> { private : - typedef typename cross_track_geo + typedef typename geographic_cross_track < FormulaPolicy >::template return_type::type return_type; public : template static inline return_type - apply(cross_track_geo const& , T const& distance) + apply(geographic_cross_track const& , T const& distance) { return distance; } @@ -281,7 +251,7 @@ struct default_strategy geographic_tag, geographic_tag > { - typedef cross_track_geo<> type; + typedef geographic_cross_track<> type; }; diff --git a/include/boost/geometry/strategies/geographic/intersection.hpp b/include/boost/geometry/strategies/geographic/intersection.hpp index f4c2c1436..84acd149c 100644 --- a/include/boost/geometry/strategies/geographic/intersection.hpp +++ b/include/boost/geometry/strategies/geographic/intersection.hpp @@ -118,6 +118,7 @@ struct geographic_segments return strategy_type(m_spheroid); } + template struct distance_strategy { typedef distance::geographic @@ -128,9 +129,10 @@ struct geographic_segments > type; }; - inline typename distance_strategy::type get_distance_strategy() const + template + inline typename distance_strategy::type get_distance_strategy() const { - typedef typename distance_strategy::type strategy_type; + typedef typename distance_strategy::type strategy_type; return strategy_type(m_spheroid); } diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 58cdb3b8b..0e611f065 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -40,13 +40,13 @@ typedef bg::strategy::distance::vincenty vincenty_pp; // Strategies for point-segment distance -typedef bg::strategy::distance::cross_track_geo +typedef bg::strategy::distance::geographic_cross_track andoyer_strategy; -typedef bg::strategy::distance::cross_track_geo +typedef bg::strategy::distance::geographic_cross_track thomas_strategy; -typedef bg::strategy::distance::cross_track_geo +typedef bg::strategy::distance::geographic_cross_track vincenty_strategy; //=========================================================================== From cd3f264a275fc72e8f7758573336ef384d4bfd89 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Thu, 31 Aug 2017 11:02:55 +0300 Subject: [PATCH 06/27] [formulas] [distance] Activate differential quantities computation in distance point segment --- include/boost/geometry/formulas/distance_point_segment.hpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index a2f32cd2a..afa4eebd7 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -25,8 +25,6 @@ https://arxiv.org/abs/1102.1215 */ -#define BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - namespace boost { namespace geometry { namespace formula { @@ -42,7 +40,7 @@ class distance_point_segment{ public: - typedef Inverse inverse_distance_type; + typedef Inverse inverse_distance_type; typedef Inverse inverse_azimuth_type; typedef Direct direct_distance_type; @@ -266,6 +264,9 @@ public: std::cout << "delta_g4=" << delta_g4 << std::endl; std::cout << "g4=" << g4 * math::r2d() << 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; From ed80e868c2c79aef02f44019963739a32e566393 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 6 Sep 2017 17:10:32 +0300 Subject: [PATCH 07/27] [formulas] Series expansion formula for distance on meridians --- .../geometry/formulas/andoyer_inverse.hpp | 56 ++++-- .../geometry/formulas/elliptic_arc_length.hpp | 178 ++++++++++++++++++ 2 files changed, 218 insertions(+), 16 deletions(-) create mode 100644 include/boost/geometry/formulas/elliptic_arc_length.hpp diff --git a/include/boost/geometry/formulas/andoyer_inverse.hpp b/include/boost/geometry/formulas/andoyer_inverse.hpp index 902fd7d8f..ef8dca9ab 100644 --- a/include/boost/geometry/formulas/andoyer_inverse.hpp +++ b/include/boost/geometry/formulas/andoyer_inverse.hpp @@ -2,6 +2,7 @@ // Copyright (c) 2015-2017 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 // Use, modification and distribution is subject to the Boost Software License, @@ -73,28 +74,44 @@ public: CT const c0 = CT(0); CT const c1 = CT(1); - CT const pi = math::pi(); CT const f = formula::flattening(spheroid); - CT const dlon = lon2 - lon1; - CT const sin_dlon = sin(dlon); - CT const cos_dlon = cos(dlon); + CT d; //the spherical distance + CT dlon; + CT cos_d; + CT cos_dlon; + CT sin_dlon; + CT cos_lat1; + CT cos_lat2; CT const sin_lat1 = sin(lat1); - CT const cos_lat1 = cos(lat1); CT const sin_lat2 = sin(lat2); - CT const cos_lat2 = cos(lat2); - // H,G,T = infinity if cos_d = 1 or cos_d = -1 - // lat1 == +-90 && lat2 == +-90 - // lat1 == lat2 && lon1 == lon2 - CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; - // on some platforms cos_d may be outside valid range - if (cos_d < -c1) - cos_d = -c1; - else if (cos_d > c1) - cos_d = c1; + if (math::equals(lon1, lon2) ) + { + d = lat1 - lat2; + cos_d = cos(d); + } + else + { + dlon = lon2 - lon1; + sin_dlon = sin(dlon); + cos_dlon = cos(dlon); + cos_lat1 = cos(lat1); + cos_lat2 = cos(lat2); + + // H,G,T = infinity if cos_d = 1 or cos_d = -1 + // lat1 == +-90 && lat2 == +-90 + // lat1 == lat2 && lon1 == lon2 + cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; + // on some platforms cos_d may be outside valid range + if (cos_d < -c1) + cos_d = -c1; + else if (cos_d > c1) + cos_d = c1; + + d = acos(cos_d); // [0, pi] + } - CT const d = acos(cos_d); // [0, pi] CT const sin_d = sin(d); // [-1, 1] if ( BOOST_GEOMETRY_CONDITION(EnableDistance) ) @@ -123,6 +140,8 @@ public: if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) ) { + CT const pi = math::pi(); + // sin_d = 0 <=> antipodal points if (math::equals(sin_d, c0)) { @@ -131,6 +150,11 @@ public: // azimuth = -inf result.azimuth = lat1 <= lat2 ? c0 : pi; } + // meridian segment + else if (math::equals(lon1, lon2)) + { + result.azimuth = c0; + } else { CT const c2 = CT(2); diff --git a/include/boost/geometry/formulas/elliptic_arc_length.hpp b/include/boost/geometry/formulas/elliptic_arc_length.hpp new file mode 100644 index 000000000..3d26aa14d --- /dev/null +++ b/include/boost/geometry/formulas/elliptic_arc_length.hpp @@ -0,0 +1,178 @@ +// Boost.Geometry + +// Copyright (c) 2017 Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, 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_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP +#define BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP + + +#include + +#include +#include + +#include +#include + +#include + + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief Compute the arc length of an ellipse. +*/ + +template +class elliptic_arc_length +{ + // Distance computation on meridians using series approximations + // to elliptic integrals. + // https://en.wikipedia.org/wiki/Meridian_arc + // latitudes are assumed to be in radians + template + CT series_expansion(T lat, Spheroid const& spheroid) + { + CT const a = get_radius<0>(spheroid); + CT const f = formula::flattening(spheroid); + CT n = f / (CT(2) - f); + CT M = a/(1+n); + CT C0 = 1; + + if (Order == 0) + { + return M * C0 * lat; + } + + CT C2 = -1.5 * n; + + if (Order == 1) + { + return M * (C0 * lat + C2 * sin(2*lat)); + } + + CT n2 = n * n; + C0 += .25 * n2; + CT C4 = 0.9375 * n2; + + if (Order == 2) + { + return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat)); + } + + CT n3 = n2 * n; + C2 += 0.1875 * n3; + CT C6 = -0.729166667 * n3; + + if (Order == 3) + { + return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat) + + C6 * sin(6*lat)); + } + + CT n4 = n2 * n2; + C4 -= 0.234375 * n4; + CT C8 = 0.615234375 * n4; + + if (Order == 4) + { + return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat) + + C6 * sin(6*lat) + C8 * sin(8*lat)); + } + + CT n5 = n4 * n; + C6 += 0.227864583 * n5; + CT C10 = -0.54140625 * n5; + + // Order 5 otr higher + return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat) + + C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat)); + + } + +public : + + template + < + unsigned int Order = 2, + typename T1, + typename T2, + typename Spheroid + > + CT apply(T1 lat1, T2 lat2, Spheroid const& spheroid) + { + return series_expansion(lat2, spheroid) + - series_expansion(lat1, spheroid); + } + + // Iterative method to elliptic arc length based on + // http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/ + // Geographic-Distance-and-Azimuth-Calculations.htm + // latitudes are assumed to be in radians + template + CT interative_method(T1 lat1, + T2 lat2, + Spheroid const& spheroid) + { + CT result = 0; + CT const zero = 0; + CT const one = 1; + CT const c1 = 2; + CT const c2 = 0.5; + CT const c3 = 4000; + + CT const a = get_radius<0>(spheroid); + CT const f = formula::flattening(spheroid); + + // how many steps to use + + CT lat1_deg = lat1 * geometry::math::r2d(); + CT lat2_deg = lat2 * geometry::math::r2d(); + + int steps = c1 + (c2 + (lat2_deg > lat1_deg) ? CT(lat2_deg - lat1_deg) + : CT(lat1_deg - lat2_deg)); + steps = (steps > c3) ? c3 : steps; + + //std::cout << "Steps=" << steps << std::endl; + + CT snLat1 = sin(lat1); + CT snLat2 = sin(lat2); + CT twoF = 2 * f - f * f; + + // limits of integration + CT x1 = a * cos(lat1) / + sqrt(1 - twoF * snLat1 * snLat1); + CT x2 = a * cos(lat2) / + sqrt(1 - twoF * snLat2 * snLat2); + + CT dx = (x2 - x1) / (steps - one); + CT x, y1, y2, dy, dydx; + CT adx = (dx < zero) ? -dx : dx; // absolute value of dx + + CT a2 = a * a; + CT oneF = 1 - f; + + // now loop through each step adding up all the little + // hypotenuses + for (int i = 0; i < (steps - 1); i++){ + x = x1 + dx * i; + dydx = ((a * oneF * sqrt((one - ((x+dx)*(x+dx))/a2))) - + (a * oneF * sqrt((one - (x*x)/a2)))) / dx; + result += adx * sqrt(one + dydx*dydx); + } + + return result; + } +}; + +}}} // namespace boost::geometry::formula + + +#endif // BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP From 3e97db55ab98de8d6ac79a04e1fddaf2c328038a Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 12 Sep 2017 14:08:50 +0300 Subject: [PATCH 08/27] [test] Tests for geographic distance of pointlike-linear geometries --- .../algorithms/distance/distance_geo_pl_l.cpp | 306 ++++++++++++++++++ 1 file changed, 306 insertions(+) diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 0e611f065..061d9e52c 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -300,13 +300,319 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, */ } +//=========================================================================== + +template +void test_distance_point_linestring(Strategy_pp const& strategy_pp, + Strategy_ps const& strategy_ps) +{ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << std::endl; + std::cout << "point/linestring distance tests" << std::endl; +#endif + typedef test_distance_of_geometries tester; + + tester::apply("p-l-01", + "POINT(0 0)", + "LINESTRING(2 0,2 0)", + pp_distance("POINT(0 0)", "POINT(2 0)", strategy_pp), + strategy_ps); + tester::apply("p-l-02", + "POINT(0 0)", + "LINESTRING(2 0,3 0)", + pp_distance("POINT(0 0)", "POINT(2 0)", strategy_pp), + strategy_ps); + tester::apply("p-l-03", + "POINT(2.5 3)", + "LINESTRING(2 0,3 0)", + pp_distance("POINT(2.5 3)", "POINT(2.5 0)", strategy_pp), + strategy_ps); + tester::apply("p-l-04", + "POINT(2 0)", + "LINESTRING(2 0,3 0)", + 0, + strategy_ps); + tester::apply("p-l-05", + "POINT(3 0)", + "LINESTRING(2 0,3 0)", + 0, + strategy_ps); + tester::apply("p-l-06", + "POINT(2.5 0)", + "LINESTRING(2 0,3 0)", + 0, + strategy_ps); + tester::apply("p-l-07", + "POINT(7.5 10)", + "LINESTRING(1 0,2 0,3 0,4 0,5 0,6 0,7 0,8 0,9 0)", + ps_distance("POINT(7.5 10)", "SEGMENT(7 0,8 0)", strategy_ps), + strategy_ps); + tester::apply("p-l-08", + "POINT(7.5 10)", + "LINESTRING(1 1,2 1,3 1,4 1,5 1,6 1,7 1,20 2,21 2)", + ps_distance("POINT(7.5 10)", "SEGMENT(7 1,20 2)", strategy_ps), + strategy_ps); +} + +void test_distance_point_linestring_strategies() +{ + typedef test_distance_of_geometries tester; + + tester::apply("p-l-03", + "POINT(2.5 3)", + "LINESTRING(2 1,3 1)", + 221147.24332788656, + vincenty_strategy()); + + tester::apply("p-l-03", + "POINT(2.5 3)", + "LINESTRING(2 1,3 1)", + 221147.36682199029, + thomas_strategy()); + + tester::apply("p-l-03", + "POINT(2.5 3)", + "LINESTRING(2 1,3 1)", + 221144.76527049288, + andoyer_strategy()); +} + +//=========================================================================== + +template +void test_distance_point_multilinestring(Strategy_pp const& strategy_pp, + Strategy_ps const& strategy_ps) +{ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << std::endl; + std::cout << "point/multilinestring distance tests" << std::endl; +#endif + typedef test_distance_of_geometries + < + point_type, multi_linestring_type + > tester; + + tester::apply("p-ml-01", + "POINT(0 0)", + "MULTILINESTRING((-5 0,-3 0),(2 0,3 0))", + pp_distance("POINT(0 0)", "POINT(2 0)", strategy_pp), + strategy_ps); + tester::apply("p-ml-02", + "POINT(2.5 3)", + "MULTILINESTRING((-5 0,-3 0),(2 0,3 0))", + pp_distance("POINT(2.5 3)", "POINT(2.5 0)", strategy_pp), + strategy_ps); + tester::apply("p-ml-03", + "POINT(2 0)", + "MULTILINESTRING((-5 0,-3 0),(2 0,3 0))", + 0, + strategy_ps); + tester::apply("p-ml-04", + "POINT(3 0)", + "MULTILINESTRING((-5 0,-3 0),(2 0,3 0))", + 0, + strategy_ps); + tester::apply("p-ml-05", + "POINT(2.5 0)", + "MULTILINESTRING((-5 0,-3 0),(2 0,3 0))", + 0, + strategy_ps); + tester::apply("p-ml-06", + "POINT(7.5 10)", + "MULTILINESTRING((-5 0,-3 0),(2 0,3 0,4 0,5 0,6 0,20 1,21 1))", + ps_distance("POINT(7.5 10)", "SEGMENT(6 0,20 1)", strategy_ps), + strategy_ps); + tester::apply("p-ml-07", + "POINT(-8 10)", + "MULTILINESTRING((-20 10,-19 11,-18 10,-6 0,-5 0,-3 0),(2 0,6 0,20 1,21 1))", + ps_distance("POINT(-8 10)", "SEGMENT(-6 0,-18 10)", strategy_ps), + strategy_ps); +} + +//=========================================================================== + +template +void test_distance_linestring_multipoint(Strategy_pp const& strategy_pp, + Strategy_ps const& strategy_ps) +{ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << std::endl; + std::cout << "linestring/multipoint distance tests" << std::endl; +#endif + typedef test_distance_of_geometries + < + linestring_type, multi_point_type + > tester; + + tester::apply("l-mp-01", + "LINESTRING(2 0,0 2,100 80)", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + ps_distance("POINT(1 1)", "SEGMENT(2 0,0 2)", strategy_ps), + strategy_ps); + + tester::apply("l-mp-02", + "LINESTRING(4 0,0 4,100 80)", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + ps_distance("POINT(1 1)", "SEGMENT(0 4,4 0)", strategy_ps), + strategy_ps); + tester::apply("l-mp-03", + "LINESTRING(1 1,2 2,100 80)", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + 0, + strategy_ps); + tester::apply("l-mp-04", + "LINESTRING(3 3,4 4,100 80)", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + pp_distance("POINT(1 1)", "POINT(3 3)", strategy_pp), + strategy_ps); + tester::apply("l-mp-05", + "LINESTRING(0 0,10 0,10 10,0 10,0 0)", + "MULTIPOINT(1 -1,80 80,5 0,150 90)", + 0, + strategy_ps); +} + +//=========================================================================== +template +void test_distance_multipoint_multilinestring(Strategy_pp const& strategy_pp, + Strategy_ps const& strategy_ps) +{ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << std::endl; + std::cout << "multipoint/multilinestring distance tests" << std::endl; +#endif + typedef test_distance_of_geometries + < + multi_point_type, multi_linestring_type + > tester; + + tester::apply("mp-ml-01", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "MULTILINESTRING((2 0,0 2),(2 2,3 3))", + ps_distance("POINT(1 1)", "SEGMENT(2 0,0 2)", strategy_ps), + strategy_ps); + tester::apply("mp-ml-02", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "MULTILINESTRING((3 0,0 3),(4 4,5 5))", + ps_distance("POINT(1 1)", "SEGMENT(3 0,0 3)", strategy_ps), + strategy_ps); + tester::apply("mp-ml-03", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "MULTILINESTRING((4 4,5 5),(1 1,2 2))", + 0, + strategy_ps); + tester::apply("mp-ml-04", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "MULTILINESTRING((4 4,3 3),(4 4,5 5))", + pp_distance("POINT(1 1)", "POINT(3 3)", strategy_pp), + strategy_ps); +} + +//=========================================================================== + +template +void test_distance_multipoint_segment(Strategy_pp const& strategy_pp, + Strategy_ps const& strategy_ps) +{ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << std::endl; + std::cout << "multipoint/segment distance tests" << std::endl; +#endif + typedef test_distance_of_geometries tester; + + tester::apply("mp-s-01", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "SEGMENT(2 0,0 2)", + ps_distance("POINT(1 1)", "SEGMENT(2 0,0 2)", strategy_ps), + strategy_ps); + tester::apply("mp-s-02", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "SEGMENT(0 -3,1 -10)", + pp_distance("POINT(0 0)", "POINT(0 -3)", strategy_pp), + strategy_ps); + tester::apply("mp-s-03", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "SEGMENT(1 1,2 2)", + 0, + strategy_ps); + tester::apply("mp-s-04", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "SEGMENT(3 3,4 4)", + pp_distance("POINT(1 1)", "POINT(3 3)", strategy_pp), + strategy_ps); + tester::apply("mp-s-05", + "MULTIPOINT(0 0,1 0,0 1,1 1)", + "SEGMENT(0.5 -3,1 -10)", + pp_distance("POINT(1 0)", "POINT(0.5 -3)", strategy_pp), + strategy_ps); +} + +//=========================================================================== + +template +void test_empty_input_pointlike_linear(Strategy const& strategy) +{ +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << std::endl; + std::cout << "testing on empty inputs... " << std::flush; +#endif + bg::model::linestring line_empty; + bg::model::multi_point multipoint_empty; + bg::model::multi_linestring > multiline_empty; + + Point point = from_wkt("POINT(0 0)"); + bg::model::linestring line = + from_wkt >("LINESTRING(0 0,1 1)"); + + // 1st geometry is empty + //test_empty_input(multipoint_empty, line, strategy); + + // 2nd geometry is empty + test_empty_input(point, line_empty, strategy); + test_empty_input(point, multiline_empty, strategy); + + // both geometries are empty + //test_empty_input(multipoint_empty, line_empty, strategy); + //test_empty_input(multipoint_empty, multiline_empty, strategy); + +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << "done!" << std::endl; +#endif +} + //=========================================================================== //=========================================================================== //=========================================================================== BOOST_AUTO_TEST_CASE( test_all_point_segment ) { + //TODO: Operations with multipoints need geographic pt-box strategy + //before activating + test_distance_point_segment(vincenty_pp(), vincenty_strategy()); test_distance_point_segment(thomas_pp(), thomas_strategy()); test_distance_point_segment(andoyer_pp(), andoyer_strategy()); + + test_distance_point_linestring(vincenty_pp(), vincenty_strategy()); + test_distance_point_linestring(thomas_pp(), thomas_strategy()); + test_distance_point_linestring(andoyer_pp(), andoyer_strategy()); + test_distance_point_linestring_strategies(); + + test_distance_point_multilinestring(vincenty_pp(), vincenty_strategy()); + test_distance_point_multilinestring(thomas_pp(), thomas_strategy()); + test_distance_point_multilinestring(andoyer_pp(), andoyer_strategy()); + + // test_distance_linestring_multipoint(vincenty_pp(), vincenty_strategy()); + // test_distance_linestring_multipoint(thomas_pp(), thomas_strategy()); + // test_distance_linestring_multipoint(andoyer_pp(), andoyer_strategy()); + + // test_distance_multipoint_multilinestring(vincenty_pp(), vincenty_strategy()); + // test_distance_multipoint_multilinestring(thomas_pp(), thomas_strategy()); + // test_distance_multipoint_multilinestring(andoyer_pp(), andoyer_strategy()); + + // test_distance_multipoint_segment(vincenty_pp(), vincenty_strategy()); + // test_distance_multipoint_segment(thomas_pp(), thomas_strategy()); + // test_distance_multipoint_segment(andoyer_pp(), andoyer_strategy()); + + test_empty_input_pointlike_linear(vincenty_strategy()); } From 47e2f0b4e96a0a4ee614921de56f46c2c24f6b72 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 12 Sep 2017 16:45:31 +0300 Subject: [PATCH 09/27] [test] Tests for small and large distances from point to segment --- .../formulas/distance_point_segment.hpp | 1 + .../algorithms/distance/distance_geo_pl_l.cpp | 49 ++++++++++--------- 2 files changed, 27 insertions(+), 23 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index afa4eebd7..4ff4dda8f 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -107,6 +107,7 @@ public: } //segment on equator + //TODO: use the meridian distance when it'll be available if (math::equals(lat1, 0) && math::equals(lat2, 0)) { if (lon3 <= lon1) diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 061d9e52c..5ea565fe8 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -250,54 +250,57 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, pp_distance("POINT(2.5 5)", "POINT(2 4)", strategy_pp), strategy_ps); - - //Nearly antipodal points not treated correctly in BG so far -/* // very small distances to segment tester::apply("p-s-07", "POINT(90 1e-3)", "SEGMENT(0.5 0,175.5 0)", - 110.57427582147321, - //1e-3 * d2r * 6335439.3272858784, - //to_comparable(strategy, 1e-3 * d2r * strategy.radius()), + pp_distance("POINT(90 0)", "POINT(90 1e-3)", strategy_pp), strategy_ps); tester::apply("p-s-08", "POINT(90 1e-4)", "SEGMENT(0.5 0,175.5 0)", - 6176.4469462110519, - //1e-4 * d2r * 6335439.3272858784, - //1e-4 * d2r * strategy.radius(), - //to_comparable(strategy, 1e-4 * d2r * strategy.radius()), + pp_distance("POINT(90 0)", "POINT(90 1e-4)", strategy_pp), strategy_ps); tester::apply("p-s-09", "POINT(90 1e-5)", "SEGMENT(0.5 0,175.5 0)", - 4917.5121192812849, - //1e-5 * d2r * strategy.radius(), - //to_comparable(strategy, 1e-5 * d2r * strategy.radius()), + pp_distance("POINT(90 0)", "POINT(90 1e-5)", strategy_pp), strategy_ps); tester::apply("p-s-10", "POINT(90 1e-6)", "SEGMENT(0.5 0,175.5 0)", - 4997.0354014473833, - //1e-6 * d2r * strategy.radius(), - //to_comparable(strategy, 1e-6 * d2r * strategy.radius()), + pp_distance("POINT(90 0)", "POINT(90 1e-6)", strategy_pp), strategy_ps); tester::apply("p-s-11", "POINT(90 1e-7)", "SEGMENT(0.5 0,175.5 0)", - 109854.20732250962, - //1e-7 * d2r * strategy.radius(), - //to_comparable(strategy, 1e-7 * d2r * strategy.radius()), + pp_distance("POINT(90 0)", "POINT(90 1e-7)", strategy_pp), strategy_ps); tester::apply("p-s-12", "POINT(90 1e-8)", "SEGMENT(0.5 0,175.5 0)", - 1372504.9365672076, - //1e-8 * d2r * strategy.radius(), - //to_comparable(strategy, 1e-8 * d2r * strategy.radius()), + pp_distance("POINT(90 0)", "POINT(90 1e-8)", strategy_pp), strategy_ps); -*/ + + // very large distance to segment + tester::apply("p-s-13", + "POINT(90 90)", + "SEGMENT(0.5 0,175.5 0)", + pp_distance("POINT(90 0)", "POINT(90 90)", strategy_pp), + strategy_ps); + tester::apply("p-s-14", + "POINT(90 90)", + "SEGMENT(0.5 -89,175.5 -89)", + pp_distance("POINT(90 -89)", "POINT(90 90)", strategy_pp), + strategy_ps); + // degenerate segment + tester::apply("p-s-14", + "POINT(90 90)", + "SEGMENT(0.5 -90,175.5 -90)", + pp_distance("POINT(90 -90)", "POINT(90 90)", strategy_pp), + strategy_ps); + + } //=========================================================================== From 4e1e18fa39d6ce0549e21e141136d54ff9ce715b Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 13 Sep 2017 11:28:24 +0300 Subject: [PATCH 10/27] [formulas] [distance] Rename inverse distance formula type name in point-segment formula --- .../formulas/distance_point_segment.hpp | 22 +++++++++---------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 4ff4dda8f..bf61b38bc 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -40,7 +40,7 @@ class distance_point_segment{ public: - typedef Inverse inverse_distance_type; + typedef Inverse inverse_distance_quantities_type; typedef Inverse inverse_azimuth_type; typedef Direct direct_distance_type; @@ -64,9 +64,9 @@ public: Spheroid const& spheroid) { result_distance_point_segment result; - result.distance = inverse_distance_type::apply(lon1, lat1, - lon2, lat2, - spheroid).distance; + result.distance = inverse_distance_quantities_type::apply(lon1, lat1, + lon2, lat2, + spheroid).distance; if (EnableClosestPoint) { result.closest_point_lon = lon1; @@ -232,8 +232,8 @@ public: CT a4 = inverse_azimuth_type::apply(res14.lon2, res14.lat2, lon2, lat2, spheroid).azimuth; - res34 = inverse_distance_type::apply(res14.lon2, res14.lat2, - lon3, lat3, spheroid); + res34 = inverse_distance_quantities_type::apply(res14.lon2, res14.lat2, + lon3, lat3, spheroid); g4 = res34.azimuth - a4; // Normalize g4 @@ -306,20 +306,20 @@ public: std::cout << "s34(sph) =" << s34_sph << std::endl; std::cout << "s34(geo) =" - << inverse_distance_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance + << inverse_distance_quantities_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance << ", p4=(" << get<0>(p4) * math::r2d() << "," << get<1>(p4) * math::r2d() << ")" << std::endl; - CT s31 = inverse_distance_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; - CT s32 = inverse_distance_type::apply(lon3, lat3, lon2, lat2, spheroid).distance; + CT s31 = inverse_distance_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; + CT s32 = inverse_distance_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 res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); - CT p4_plus = inverse_distance_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; + CT p4_plus = inverse_distance_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; geometry::formula::result_direct res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a1, spheroid); - CT p4_minus = inverse_distance_type::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; + CT p4_minus = inverse_distance_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() << "," << res4.lat2 * math::r2d() << ")" From 2bc7063e73741db292ca1d0b8d19192dd6fa2420 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 13 Sep 2017 18:06:53 +0300 Subject: [PATCH 11/27] [formulas] [distance] Use mean radius formula in distance point segment formula --- include/boost/geometry/formulas/distance_point_segment.hpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index bf61b38bc..4bd6a9a48 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -12,6 +12,7 @@ #define BOOST_GEOMETRY_FORMULAS_CROSS_TRACK_GEO_HPP #include +#include #ifndef BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS #define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100 @@ -82,8 +83,7 @@ public: CT lon3, CT lat3, //query point p3 Spheroid const& spheroid) { - CT earth_radius = (CT(2) * spheroid.get_radius<1>() - + spheroid.get_radius<2>()) / CT(3); + CT const earth_radius = geometry::formula::mean_radius(spheroid); result_distance_point_segment result; From 02de7cc212e462ae9b0cb7b69713660b2681de40 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 19 Sep 2017 15:13:33 +0300 Subject: [PATCH 12/27] [formulas] [strategies] Fix cross track distance for both spherical and geographic --- .../formulas/distance_point_segment.hpp | 100 ++++++++++++------ .../spherical/distance_cross_track.hpp | 54 +++++++--- .../algorithms/distance/distance_geo_pl_l.cpp | 5 + test/algorithms/distance/distance_se_pl_l.cpp | 9 ++ 4 files changed, 118 insertions(+), 50 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 4bd6a9a48..4fc4e7ad4 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -43,6 +43,7 @@ public: typedef Inverse inverse_distance_quantities_type; typedef Inverse inverse_azimuth_type; + typedef Inverse inverse_azimuth_reverse_type; typedef Direct direct_distance_type; struct result_distance_point_segment @@ -110,6 +111,9 @@ public: //TODO: use the meridian distance when it'll be available if (math::equals(lat1, 0) && math::equals(lat2, 0)) { +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "Equatorial segment" << std::endl; +#endif if (lon3 <= lon1) { return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); @@ -121,51 +125,52 @@ public: return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); } - //segment on meridian - if (math::equals(lon1, lon2)) + //TODO no quantities needed ? + CT d1 = inverse_distance_quantities_type::apply(lon1, lat1, + lon3, lat3, spheroid).distance; + CT d3 = inverse_distance_quantities_type::apply(lon1, lat1, + lon2, lat2, spheroid).distance; + + if (geometry::math::equals(d3, 0.0)) { - if (lat1 > lat2) + // "Degenerate" segment, return either d1 or d2 + result_distance_point_segment result; + result.distance = d1; + if (EnableClosestPoint) { - std::swap(lat1, lat2); + result.closest_point_lon = lon1; + result.closest_point_lat = lat1; } - if (lat3 <= lat1) - { - return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); - } - if (lat3 >= lat2) - { - return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); - } - return non_iterative_case(lon1, lat3, lon3, lat3, spheroid); + return result; } - // Easy cases - if (lon3 <= lon1 && lat3 <= lat1) - { - return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); - } - if (lon3 >= lon2 && lat3 <= lat2) - { - return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); - } + CT d2 = inverse_distance_quantities_type::apply(lon2, lat2, + lon3, lat3, spheroid).distance; // Compute a1 (GEO) - CT a1 = inverse_azimuth_type::apply(lon1, lat1, lon2, lat2, spheroid).azimuth; + geometry::formula::result_inverse res12 = + inverse_azimuth_reverse_type::apply(lon1, lat1, lon2, lat2, spheroid); + CT a1 = res12.azimuth; CT a13 = inverse_azimuth_type::apply(lon1, lat1, lon3, lat3, spheroid).azimuth; - CT a312 = a1 - a13; + CT a312 = a13 - a1; - if (a312 > half_pi) + if (geometry::math::equals(a312, 0.0)) { - // 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); +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "point on segment" << std::endl; +#endif + result_distance_point_segment result; + result.distance = 0.0; + if (EnableClosestPoint) + { + result.closest_point_lon = lon3; + result.closest_point_lat = lat3; + } + return result; } - CT a2 = pi + inverse_azimuth_type::apply(lon2, lat2, lon1, lat1, spheroid).azimuth; - CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; - - CT a321 = a2 - a23; + CT projection1 = cos( a312 ) * d1 / d3; #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG std::cout << "segment=(" << lon1 * math::r2d(); @@ -177,12 +182,41 @@ public: std::cout << ")\na1=" << a1 * math::r2d() << std::endl; std::cout << "a13=" << a13 * math::r2d() << std::endl; std::cout << "a312=" << a312 * math::r2d() << std::endl; + std::cout << "cos(a312)=" << cos(a312) << std::endl; +#endif + //if (a312 > half_pi) + if (projection1 < 0.0) + { +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "projection closer to p1" << std::endl; +#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 a2 = pi + inverse_azimuth_type::apply(lon2, lat2, lon1, lat1, spheroid).azimuth; +// CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; + + CT a2 = res12.reverse_azimuth - pi; + CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; + + CT a321 = a23 - a2; + +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG std::cout << "a2=" << a2 * math::r2d() << std::endl; std::cout << "a23=" << a13 * math::r2d() << std::endl; std::cout << "a321=" << a321 * math::r2d() << std::endl; + std::cout << "cos(a321)=" << cos(a321) << std::endl; #endif - if (a321 < half_pi) + CT projection2 = cos( a321 ) * d2 / d3; + + //if (a321 < half_pi) + if (projection2 < 0.0) { +#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG + std::cout << "projection closer to p2" << std::endl; +#endif // projection of p3 on geodesic spanned by segment (p1,p2) fall // outside of segment on the side of p2 return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); diff --git a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp index 7daafa4a1..8f7df9b2a 100644 --- a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -373,19 +373,6 @@ public : typedef typename return_type::type return_type; -#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK - std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " - << crs_AD * geometry::math::r2d() << std::endl; - std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " - << crs_AB * geometry::math::r2d() << std::endl; - std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " - << crs_BD * geometry::math::r2d << std::endl; - std::cout << "Projection AD-AB " << projection1 << " : " - << d_crs1 * geometry::math::r2d() << std::endl; - std::cout << "Projection BD-BA " << projection2 << " : " - << d_crs2 * geometry::math::r2d() << std::endl; -#endif - // http://williams.best.vwh.net/avform.htm#XTE return_type d1 = m_strategy.apply(sp1, p); return_type d3 = m_strategy.apply(sp1, sp2); @@ -398,10 +385,25 @@ public : return_type d2 = m_strategy.apply(sp2, p); - return_type crs_AD = geometry::detail::course(sp1, p); - return_type crs_AB = geometry::detail::course(sp1, sp2); - return_type crs_BA = crs_AB - geometry::math::pi(); - return_type crs_BD = geometry::detail::course(sp2, p); + return_type lon1 = geometry::get_as_radian<0>(sp1); + return_type lat1 = geometry::get_as_radian<1>(sp1); + return_type lon2 = geometry::get_as_radian<0>(sp2); + return_type lat2 = geometry::get_as_radian<1>(sp2); + return_type lon = geometry::get_as_radian<0>(p); + return_type lat = geometry::get_as_radian<1>(p); + + return_type crs_AD = geometry::formula::spherical_azimuth + (lon1, lat1, lon, lat).azimuth; + + geometry::formula::result_spherical result = + geometry::formula::spherical_azimuth + (lon1, lat1, lon2, lat2); + return_type crs_AB = result.azimuth; + return_type crs_BA = result.reverse_azimuth - geometry::math::pi(); + + return_type crs_BD = geometry::formula::spherical_azimuth + (lon2, lat2, lon, lat).azimuth; + return_type d_crs1 = crs_AD - crs_AB; return_type d_crs2 = crs_BD - crs_BA; @@ -409,6 +411,24 @@ public : return_type projection1 = cos( d_crs1 ) * d1 / d3; return_type projection2 = cos( d_crs2 ) * d2 / d3; +#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK + std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " " + << crs_AD * geometry::math::r2d() << std::endl; + std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " " + << crs_AB * geometry::math::r2d() << std::endl; + std::cout << "Course " << dsv(sp2) << " to " << dsv(sp1) << " " + << crs_BA * geometry::math::r2d() << std::endl; + std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " " + << crs_BD * geometry::math::r2d() << std::endl; + std::cout << "Projection AD-AB " << projection1 << " : " + << d_crs1 * geometry::math::r2d() << std::endl; + std::cout << "Projection BD-BA " << projection2 << " : " + << d_crs2 * geometry::math::r2d() << std::endl; + std::cout << " d1: " << (d1 ) + << " d2: " << (d2 ) + << std::endl; +#endif + if (projection1 > 0.0 && projection2 > 0.0) { #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 5ea565fe8..5efecc171 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -249,6 +249,11 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, "SEGMENT(2 2,2 4)", pp_distance("POINT(2.5 5)", "POINT(2 4)", strategy_pp), strategy_ps); + tester::apply("p-s-mer4", + "POINT(1 80)", + "SEGMENT(0 0,0 90)", + pp_distance("POINT(1 80)", "POINT(0 80.00149225834545)", andoyer_pp()), + andoyer_strategy()); // very small distances to segment tester::apply("p-s-07", diff --git a/test/algorithms/distance/distance_se_pl_l.cpp b/test/algorithms/distance/distance_se_pl_l.cpp index f298bae25..5b3dfd5a9 100644 --- a/test/algorithms/distance/distance_se_pl_l.cpp +++ b/test/algorithms/distance/distance_se_pl_l.cpp @@ -10,6 +10,7 @@ // http://www.boost.org/users/license.html #include +#define BOOST_GEOMETRY_TEST_DEBUG #ifndef BOOST_TEST_MODULE #define BOOST_TEST_MODULE test_distance_spherical_equatorial_pl_l @@ -160,6 +161,14 @@ void test_distance_point_segment(Strategy const& strategy) "POINT(3.5 3)", strategy), strategy); + tester::apply("p-s-07", + "POINT(0 0)", + "SEGMENT(0 10,10 10)", + ps_distance("POINT(0 0)", "SEGMENT(10 10,0 10)", strategy), + pp_comparable_distance("POINT(0 0)", + "POINT(0 10)", + strategy), + strategy); // very small distances to segment tester::apply("p-s-07", "POINT(90 1e-3)", From 3b05e5ca4604aeab7e716a3e52ee95b39a5acdd6 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 19 Sep 2017 17:25:47 +0300 Subject: [PATCH 13/27] [strategies] [formulas] Replace course algorithm by azimuth formula --- .../geographic/distance_cross_track.hpp | 2 -- .../spherical/distance_cross_track.hpp | 7 +++-- .../spherical/side_by_cross_track.hpp | 26 +++++++++++++++---- 3 files changed, 24 insertions(+), 11 deletions(-) diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp index 295858f6e..4a220d1c9 100644 --- a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -23,8 +23,6 @@ #include #include -#include - #include #include #include diff --git a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp index 8f7df9b2a..82a2fb3d6 100644 --- a/include/boost/geometry/strategies/spherical/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp @@ -2,9 +2,10 @@ // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2017. +// Modifications copyright (c) 2014-2017, Oracle and/or its affiliates. +// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, @@ -26,8 +27,6 @@ #include #include -#include - #include #include #include diff --git a/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp b/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp index 3f7be0555..1f14ce9aa 100644 --- a/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp +++ b/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp @@ -2,9 +2,10 @@ // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands. -// This file was modified by Oracle on 2014. -// Modifications copyright (c) 2014, Oracle and/or its affiliates. +// This file was modified by Oracle on 2014-2017. +// Modifications copyright (c) 2014-2017, 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 // Use, modification and distribution is subject to the Boost Software License, @@ -18,7 +19,7 @@ #include #include -#include +#include #include #include @@ -59,8 +60,23 @@ public : >::type calc_t; calc_t d1 = 0.001; // m_strategy.apply(sp1, p); - calc_t crs_AD = geometry::detail::course(p1, p); - calc_t crs_AB = geometry::detail::course(p1, p2); + + calc_t lon1 = geometry::get_as_radian<0>(p1); + calc_t lat1 = geometry::get_as_radian<1>(p1); + calc_t lon2 = geometry::get_as_radian<0>(p2); + calc_t lat2 = geometry::get_as_radian<1>(p2); + calc_t lon = geometry::get_as_radian<0>(p); + calc_t lat = geometry::get_as_radian<1>(p); + + calc_t crs_AD = geometry::formula::spherical_azimuth + (lon1, lat1, lon, lat).azimuth; + + //calc_t crs_AD = geometry::detail::course(p1, p); + + calc_t crs_AB = geometry::formula::spherical_azimuth + (lon1, lat1, lon2, lat2).azimuth; + + //calc_t crs_AB = geometry::detail::course(p1, p2); calc_t XTD = asin(sin(d1) * sin(crs_AD - crs_AB)); return math::equals(XTD, 0) ? 0 : XTD < 0 ? 1 : -1; From 19eda49034fc4c594c6b920c8e350ee3205f34cc Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 19 Sep 2017 17:28:08 +0300 Subject: [PATCH 14/27] [strategies] [formulas] Some code cleaning --- include/boost/geometry/formulas/distance_point_segment.hpp | 5 ----- .../geometry/strategies/spherical/side_by_cross_track.hpp | 3 --- 2 files changed, 8 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 4fc4e7ad4..790823f4a 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -184,7 +184,6 @@ public: std::cout << "a312=" << a312 * math::r2d() << std::endl; std::cout << "cos(a312)=" << cos(a312) << std::endl; #endif - //if (a312 > half_pi) if (projection1 < 0.0) { #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG @@ -195,9 +194,6 @@ public: return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } -// CT a2 = pi + inverse_azimuth_type::apply(lon2, lat2, lon1, lat1, spheroid).azimuth; -// CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; - CT a2 = res12.reverse_azimuth - pi; CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; @@ -211,7 +207,6 @@ public: #endif CT projection2 = cos( a321 ) * d2 / d3; - //if (a321 < half_pi) if (projection2 < 0.0) { #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG diff --git a/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp b/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp index 1f14ce9aa..2b195e96f 100644 --- a/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp +++ b/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp @@ -71,12 +71,9 @@ public : calc_t crs_AD = geometry::formula::spherical_azimuth (lon1, lat1, lon, lat).azimuth; - //calc_t crs_AD = geometry::detail::course(p1, p); - calc_t crs_AB = geometry::formula::spherical_azimuth (lon1, lat1, lon2, lat2).azimuth; - //calc_t crs_AB = geometry::detail::course(p1, p2); calc_t XTD = asin(sin(d1) * sin(crs_AD - crs_AB)); return math::equals(XTD, 0) ? 0 : XTD < 0 ? 1 : -1; From d9f3641795fc370b5508974c7f74db1b0ec32d0d Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 19 Sep 2017 18:01:43 +0300 Subject: [PATCH 15/27] [formulas] [strategies] [tests] Code simplification for special cases and some unit test cases --- .../formulas/distance_point_segment.hpp | 91 +++++++++---------- .../algorithms/distance/distance_geo_pl_l.cpp | 19 ++++ 2 files changed, 64 insertions(+), 46 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 790823f4a..1462fd8dd 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -42,6 +42,7 @@ 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; @@ -59,24 +60,34 @@ public: CT closest_point_lat; }; + result_distance_point_segment + static inline non_iterative_case(CT lon, CT lat, CT distance) + { + result_distance_point_segment result; + result.distance = distance; + + if (EnableClosestPoint) + { + result.closest_point_lon = lon; + result.closest_point_lat = lat; + } + return result; + } + template result_distance_point_segment static inline non_iterative_case(CT lon1, CT lat1, //p1 CT lon2, CT lat2, //p2 Spheroid const& spheroid) { - result_distance_point_segment result; - result.distance = inverse_distance_quantities_type::apply(lon1, lat1, - lon2, lat2, - spheroid).distance; - if (EnableClosestPoint) - { - result.closest_point_lon = lon1; - result.closest_point_lat = lat1; - } - return result; + CT distance = inverse_distance_quantities_type::apply(lon1, lat1, + lon2, lat2, + spheroid).distance; + return non_iterative_case(lon1, lat1, distance); } + + template result_distance_point_segment static inline apply(CT lon1, CT lat1, //p1 @@ -92,6 +103,7 @@ public: CT const f = flattening(spheroid); CT const pi = math::pi(); CT const half_pi = pi / CT(2); + CT const c0 = CT(0); // Convert to radians lon1 = math::as_radian(lon1); @@ -125,49 +137,33 @@ public: return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); } - //TODO no quantities needed ? - CT d1 = inverse_distance_quantities_type::apply(lon1, lat1, - lon3, lat3, spheroid).distance; - CT d3 = inverse_distance_quantities_type::apply(lon1, lat1, - lon2, lat2, spheroid).distance; + CT d1 = inverse_distance_type::apply(lon1, lat1, + lon3, lat3, spheroid).distance; + CT d3 = inverse_distance_type::apply(lon1, lat1, + lon2, lat2, spheroid).distance; - if (geometry::math::equals(d3, 0.0)) + if (geometry::math::equals(d3, c0)) { - // "Degenerate" segment, return either d1 or d2 - result_distance_point_segment result; - result.distance = d1; - if (EnableClosestPoint) - { - result.closest_point_lon = lon1; - result.closest_point_lat = lat1; - } - return result; + return non_iterative_case(lon1, lat2, d1); } - CT d2 = inverse_distance_quantities_type::apply(lon2, lat2, - lon3, lat3, spheroid).distance; + CT d2 = inverse_distance_type::apply(lon2, lat2, + lon3, lat3, spheroid).distance; - // Compute a1 (GEO) + // Compute a12 (GEO) geometry::formula::result_inverse res12 = inverse_azimuth_reverse_type::apply(lon1, lat1, lon2, lat2, spheroid); - CT a1 = res12.azimuth; + CT a12 = res12.azimuth; CT a13 = inverse_azimuth_type::apply(lon1, lat1, lon3, lat3, spheroid).azimuth; - CT a312 = a13 - a1; + CT a312 = a13 - a12; - if (geometry::math::equals(a312, 0.0)) + if (geometry::math::equals(a312, c0)) { #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG std::cout << "point on segment" << std::endl; #endif - result_distance_point_segment result; - result.distance = 0.0; - if (EnableClosestPoint) - { - result.closest_point_lon = lon3; - result.closest_point_lat = lat3; - } - return result; + return non_iterative_case(lon3, lat3, c0); } CT projection1 = cos( a312 ) * d1 / d3; @@ -179,7 +175,7 @@ public: std::cout << "," << lat2 * math::r2d(); std::cout << ")\np=(" << lon3 * math::r2d(); std::cout << "," << lat3 * math::r2d(); - std::cout << ")\na1=" << a1 * math::r2d() << std::endl; + std::cout << ")\na1=" << a12 * math::r2d() << std::endl; std::cout << "a13=" << a13 * math::r2d() << std::endl; std::cout << "a312=" << a312 * math::r2d() << std::endl; std::cout << "cos(a312)=" << cos(a312) << std::endl; @@ -194,13 +190,13 @@ public: return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); } - CT a2 = res12.reverse_azimuth - pi; + CT a21 = res12.reverse_azimuth - pi; CT a23 = inverse_azimuth_type::apply(lon2, lat2, lon3, lat3, spheroid).azimuth; - CT a321 = a23 - a2; + CT a321 = a23 - a21; #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "a2=" << a2 * math::r2d() << std::endl; + std::cout << "a21=" << a21 * math::r2d() << std::endl; std::cout << "a23=" << a13 * math::r2d() << std::endl; std::cout << "a321=" << a321 * math::r2d() << std::endl; std::cout << "cos(a321)=" << cos(a321) << std::endl; @@ -218,7 +214,10 @@ public: } // Guess s14 (SPHERICAL) - typedef geometry::model::point > point; + typedef geometry::model::point< + CT, 2, + geometry::cs::spherical_equatorial + > point; CT bet1 = atan((1 - f) * tan(lon1)); CT bet2 = atan((1 - f) * tan(lon2)); @@ -254,7 +253,7 @@ public: prev_distance = res34.distance; // Solve the direct problem to find p4 (GEO) - res14 = direct_distance_type::apply(lon1, lat1, s14, a1, spheroid); + 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) @@ -347,7 +346,7 @@ public: geometry::formula::result_direct res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); CT p4_plus = inverse_distance_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; - geometry::formula::result_direct res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a1, spheroid); + geometry::formula::result_direct res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a12, spheroid); CT p4_minus = inverse_distance_quantities_type::apply(res1.lon2, res1.lat2, lon3, lat3, spheroid).distance; std::cout << "s31=" << s31 << "\ns32=" << s32 diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 5efecc171..f687ed293 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -109,6 +109,11 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, "SEGMENT(2 0,3 0)", 0, strategy_ps); + tester::apply("p-s-05", + "POINT(2.5 0)", + "SEGMENT(2 0,3 0)", + 0, + strategy_ps); tester::apply("p-s-06", "POINT(3.5 3)", "SEGMENT(2 0,3 0)", @@ -255,6 +260,20 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, pp_distance("POINT(1 80)", "POINT(0 80.00149225834545)", andoyer_pp()), andoyer_strategy()); + //degenerate segment + tester::apply("p-s-deg", + "POINT(1 80)", + "SEGMENT(0 0,0 0)", + pp_distance("POINT(0 0)", "POINT(1 80)", strategy_pp), + strategy_ps); + + //point on segment + tester::apply("p-s-deg", + "POINT(0 80)", + "SEGMENT(0 0,0 90)", + 0, + strategy_ps); + // very small distances to segment tester::apply("p-s-07", "POINT(90 1e-3)", From 866b297482138a8fe7976b27603ced3db05e7900 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Thu, 21 Sep 2017 13:22:51 +0100 Subject: [PATCH 16/27] [formulas] [tests] Distance point-segment; normalization and some changes in the method and tests --- .../formulas/distance_point_segment.hpp | 50 ++++++++++--------- .../algorithms/distance/distance_geo_pl_l.cpp | 4 +- 2 files changed, 28 insertions(+), 26 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 1462fd8dd..1700a8842 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -86,7 +86,24 @@ public: return non_iterative_case(lon1, lat1, distance); } - + template + T static inline normalize(T g4) + { + CT const pi = math::pi(); + if (g4 < 0 && g4 < -pi)//close to -270 + { + return g4 + 1.5 * pi; + } + else if (g4 > 0 && g4 > pi)//close to 270 + { + return - g4 + 1.5 * pi; + } + else if (g4 < 0 && g4 > -pi)//close to -90 + { + return -g4 - pi/2; + } + return g4 - pi/2; + } template result_distance_point_segment @@ -121,7 +138,7 @@ public: //segment on equator //TODO: use the meridian distance when it'll be available - if (math::equals(lat1, 0) && math::equals(lat2, 0)) + if (math::equals(lat1, c0) && math::equals(lat2, c0)) { #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG std::cout << "Equatorial segment" << std::endl; @@ -214,10 +231,11 @@ public: } // Guess s14 (SPHERICAL) - typedef geometry::model::point< - CT, 2, - geometry::cs::spherical_equatorial - > point; + typedef geometry::model::point + < + CT, 2, + geometry::cs::spherical_equatorial + > point; CT bet1 = atan((1 - f) * tan(lon1)); CT bet2 = atan((1 - f) * tan(lon2)); @@ -264,27 +282,11 @@ public: lon3, lat3, spheroid); g4 = res34.azimuth - a4; - // Normalize g4 - if (g4 < 0 && g4 < -pi)//close to -270 - { - delta_g4 = g4 + 1.5 * pi; - } - else if (g4 > 0 && g4 > pi)//close to 270 - { - delta_g4 = - g4 + 1.5 * pi; - } - else if (g4 < 0 && g4 > -pi)//close to -90 - { - delta_g4 = -g4 - pi/2; - } - else //close to 90 - { - delta_g4 = g4 - pi/2; - } + delta_g4 = normalize(g4); CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit CT m34 = res34.reduced_length; - CT der = M43 / m34; + CT der = (M43 / m34) * sin(g4); s14 = s14 - delta_g4 / der; #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index f687ed293..6da51ca55 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -247,8 +247,8 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, tester::apply("p-s-mer2", "POINT(2.5 3)", "SEGMENT(2 2,2 4)", - pp_distance("POINT(2.5 3)", "POINT(2 3)", strategy_pp), - strategy_ps); + pp_distance("POINT(2.5 3)", "POINT(2 3.000114792872075)", andoyer_pp()), + andoyer_strategy()); tester::apply("p-s-mer3", "POINT(2.5 5)", "SEGMENT(2 2,2 4)", From eca31647a81282e40e2d5c65a48f317587b3341b Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Thu, 28 Sep 2017 16:45:46 +0100 Subject: [PATCH 17/27] [distance] [test] Meridian distance computation actived in distance strategy --- .../geometry/formulas/andoyer_inverse.hpp | 58 +++------ .../geometry/formulas/elliptic_arc_length.hpp | 32 ++--- .../strategies/geographic/distance.hpp | 53 +++++++- .../strategies/geographic/parameters.hpp | 1 + test/strategies/andoyer.cpp | 10 +- test/strategies/distance.cpp | 117 ++++++++++++++++++ 6 files changed, 197 insertions(+), 74 deletions(-) create mode 100644 test/strategies/distance.cpp diff --git a/include/boost/geometry/formulas/andoyer_inverse.hpp b/include/boost/geometry/formulas/andoyer_inverse.hpp index ef8dca9ab..5c195dbdb 100644 --- a/include/boost/geometry/formulas/andoyer_inverse.hpp +++ b/include/boost/geometry/formulas/andoyer_inverse.hpp @@ -2,7 +2,6 @@ // Copyright (c) 2015-2017 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 // Use, modification and distribution is subject to the Boost Software License, @@ -74,44 +73,28 @@ public: CT const c0 = CT(0); CT const c1 = CT(1); + CT const pi = math::pi(); CT const f = formula::flattening(spheroid); - CT d; //the spherical distance - CT dlon; - CT cos_d; - CT cos_dlon; - CT sin_dlon; - CT cos_lat1; - CT cos_lat2; + CT const dlon = lon2 - lon1; + CT const sin_dlon = sin(dlon); + CT const cos_dlon = cos(dlon); CT const sin_lat1 = sin(lat1); + CT const cos_lat1 = cos(lat1); CT const sin_lat2 = sin(lat2); + CT const cos_lat2 = cos(lat2); - if (math::equals(lon1, lon2) ) - { - d = lat1 - lat2; - cos_d = cos(d); - } - else - { - dlon = lon2 - lon1; - sin_dlon = sin(dlon); - cos_dlon = cos(dlon); - cos_lat1 = cos(lat1); - cos_lat2 = cos(lat2); - - // H,G,T = infinity if cos_d = 1 or cos_d = -1 - // lat1 == +-90 && lat2 == +-90 - // lat1 == lat2 && lon1 == lon2 - cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; - // on some platforms cos_d may be outside valid range - if (cos_d < -c1) - cos_d = -c1; - else if (cos_d > c1) - cos_d = c1; - - d = acos(cos_d); // [0, pi] - } + // H,G,T = infinity if cos_d = 1 or cos_d = -1 + // lat1 == +-90 && lat2 == +-90 + // lat1 == lat2 && lon1 == lon2 + CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon; + // on some platforms cos_d may be outside valid range + if (cos_d < -c1) + cos_d = -c1; + else if (cos_d > c1) + cos_d = c1; + CT const d = acos(cos_d); // [0, pi] CT const sin_d = sin(d); // [-1, 1] if ( BOOST_GEOMETRY_CONDITION(EnableDistance) ) @@ -140,8 +123,6 @@ public: if ( BOOST_GEOMETRY_CONDITION(CalcAzimuths) ) { - CT const pi = math::pi(); - // sin_d = 0 <=> antipodal points if (math::equals(sin_d, c0)) { @@ -150,11 +131,6 @@ public: // azimuth = -inf result.azimuth = lat1 <= lat2 ? c0 : pi; } - // meridian segment - else if (math::equals(lon1, lon2)) - { - result.azimuth = c0; - } else { CT const c2 = CT(2); @@ -238,7 +214,7 @@ private: static inline void normalize_azimuth(CT & azimuth, CT const& A, CT const& dA) { CT const c0 = 0; - + if (A >= c0) // A indicates Eastern hemisphere { if (dA >= c0) // A altered towards 0 diff --git a/include/boost/geometry/formulas/elliptic_arc_length.hpp b/include/boost/geometry/formulas/elliptic_arc_length.hpp index 3d26aa14d..057085eae 100644 --- a/include/boost/geometry/formulas/elliptic_arc_length.hpp +++ b/include/boost/geometry/formulas/elliptic_arc_length.hpp @@ -30,15 +30,18 @@ namespace boost { namespace geometry { namespace formula \brief Compute the arc length of an ellipse. */ -template +template class elliptic_arc_length { + +public : + // Distance computation on meridians using series approximations - // to elliptic integrals. + // to elliptic integrals. Formula to compute distance from lattitude 0 to lat // https://en.wikipedia.org/wiki/Meridian_arc - // latitudes are assumed to be in radians - template - CT series_expansion(T lat, Spheroid const& spheroid) + // latitudes are assumed to be in radians and in [-pi/2,pi/2] + template + static CT apply(T lat, Spheroid const& spheroid) { CT const a = get_radius<0>(spheroid); CT const f = formula::flattening(spheroid); @@ -91,31 +94,16 @@ class elliptic_arc_length C6 += 0.227864583 * n5; CT C10 = -0.54140625 * n5; - // Order 5 otr higher + // Order 5 or higher return M * (C0 * lat + C2 * sin(2*lat) + C4 * sin(4*lat) + C6 * sin(6*lat) + C8 * sin(8*lat) + C10 * sin(10*lat)); } -public : - - template - < - unsigned int Order = 2, - typename T1, - typename T2, - typename Spheroid - > - CT apply(T1 lat1, T2 lat2, Spheroid const& spheroid) - { - return series_expansion(lat2, spheroid) - - series_expansion(lat1, spheroid); - } - // Iterative method to elliptic arc length based on // http://www.codeguru.com/cpp/cpp/algorithms/article.php/c5115/ // Geographic-Distance-and-Azimuth-Calculations.htm - // latitudes are assumed to be in radians + // latitudes are assumed to be in radians and in [-pi/2,pi/2] template CT interative_method(T1 lat1, T2 lat2, diff --git a/include/boost/geometry/strategies/geographic/distance.hpp b/include/boost/geometry/strategies/geographic/distance.hpp index d3656f449..bd0e0d5fc 100644 --- a/include/boost/geometry/strategies/geographic/distance.hpp +++ b/include/boost/geometry/strategies/geographic/distance.hpp @@ -5,6 +5,7 @@ // This file was modified by Oracle on 2014-2017. // Modifications copyright (c) 2014-2017 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 // Use, modification and distribution is subject to the Boost Software License, @@ -28,6 +29,7 @@ #include #include +#include #include #include @@ -74,13 +76,52 @@ public : inline typename calculation_type::type apply(Point1 const& point1, Point2 const& point2) const { + typedef typename calculation_type::type CT; + + CT lon1 = get_as_radian<0>(point1); + CT lat1 = get_as_radian<1>(point1); + 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 - < - typename calculation_type::type, - true, false, false, false, false - >::apply(get_as_radian<0>(point1), get_as_radian<1>(point1), - get_as_radian<0>(point2), get_as_radian<1>(point2), - m_spheroid).distance; + < + CT, true, false, false, false, false + >::apply(lon1, lat1, lon2, lat2, m_spheroid).distance; } inline Spheroid const& model() const diff --git a/include/boost/geometry/strategies/geographic/parameters.hpp b/include/boost/geometry/strategies/geographic/parameters.hpp index 5638db50f..c70b0cabe 100644 --- a/include/boost/geometry/strategies/geographic/parameters.hpp +++ b/include/boost/geometry/strategies/geographic/parameters.hpp @@ -14,6 +14,7 @@ #include #include #include +#include #include #include diff --git a/test/strategies/andoyer.cpp b/test/strategies/andoyer.cpp index 01b31c73d..0b4fff12c 100644 --- a/test/strategies/andoyer.cpp +++ b/test/strategies/andoyer.cpp @@ -214,10 +214,10 @@ void test_all() // antipodal // ok? in those cases shorter path would pass through a pole // but 90 or -90 would be consistent with distance? - test_distazi(0, 0, 180, 0, 20037.5, 0.0); - test_distazi(0, 0, -180, 0, 20037.5, 0.0); - test_distazi(-90, 0, 90, 0, 20037.5, 0.0); - test_distazi(90, 0, -90, 0, 20037.5, 0.0); + //test_distazi(0, 0, 180, 0, 20037.5, 0.0); + //test_distazi(0, 0, -180, 0, 20037.5, 0.0); + //test_distazi(-90, 0, 90, 0, 20037.5, 0.0); + //test_distazi(90, 0, -90, 0, 20037.5, 0.0); // 0, 45, 90 ... for (int i = 0 ; i < 360 ; i += 45) @@ -264,7 +264,7 @@ void test_all() test_distazi_symm(normlized_deg(l-44.99), -44.99, normlized_deg(l+135), 45, 20008.1, 0.0); test_distazi_symm(normlized_deg(l-44.999), -44.999, normlized_deg(l+135), 45, 20009.4, 0.0); // antipodal - test_distazi_symm(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20020.7, 0.0); + //test_distazi_symm(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20020.7, 0.0); } /* SQL Server gives: diff --git a/test/strategies/distance.cpp b/test/strategies/distance.cpp new file mode 100644 index 000000000..dfe47b9ec --- /dev/null +++ b/test/strategies/distance.cpp @@ -0,0 +1,117 @@ +// Boost.Geometry + +// Copyright (c) 2017 Oracle and/or its affiliates. + +// Contributed and/or modified by Vissarion Fysikopoulos, 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) + +#include + +#include + +#include +#include +#include +#include + +#include + +#ifdef HAVE_TTMATH +# include +#endif + +typedef bg::srs::spheroid stype; + +typedef bg::strategy::andoyer andoyer_formula; +typedef bg::strategy::thomas thomas_formula; +typedef bg::strategy::vincenty vincenty_formula; + +template +bool non_precise_ct() +{ + typedef typename bg::coordinate_type

::type ct; + return boost::is_integral::value || boost::is_float::value; +} + +template +void test_distance(double lon1, double lat1, double lon2, double lat2) +{ + typedef typename bg::promote_floating_point + < + typename bg::select_calculation_type::type + >::type calc_t; + + calc_t tolerance = non_precise_ct() || non_precise_ct() ? + 5.0 : 0.001; + + P1 p1; + P2 p2; + + bg::assign_values(p1, lon1, lat1); + bg::assign_values(p2, lon2, lat2); + + // Test strategy that implements meridian distance against formula + // that implements general distance + // That may change in the future but in any case these calls must return + // the same result + + calc_t dist_formula = FormulaPolicy::template inverse + < + double, true, false, false, false, false + >::apply(lon1 * bg::math::d2r(), + lat1 * bg::math::d2r(), + lon2 * bg::math::d2r(), + lat2 * bg::math::d2r(), + stype()).distance; + + bg::strategy::distance::geographic strategy; + calc_t dist_strategy = strategy.apply(p1, p2); + BOOST_CHECK_CLOSE(dist_formula, dist_strategy, tolerance); +} + +template +void test_distance_reverse(double lon1, double lat1, + double lon2, double lat2) +{ + test_distance(lon1, lat1, lon2, lat2); + test_distance(lon2, lat2, lon1, lat1); +} + +template +void test_meridian() +{ + test_distance_reverse(0., 70., 0., 80.); + test_distance_reverse(0, 70, 0., -80.); + test_distance_reverse(0., -70., 0., 80.); + test_distance_reverse(0., -70., 0., -80.); + + test_distance_reverse(0., 70., 180., 80.); + test_distance_reverse(0., 70., 180., -80.); + test_distance_reverse(0., -70., 180., 80.); + test_distance_reverse(0., -70., 180., -80.); + + test_distance_reverse(350., 70., 170., 80.); + test_distance_reverse(350., 70., 170., -80.); + test_distance_reverse(350., -70., 170., 80.); + test_distance_reverse(350., -70., 170., -80.); +} + +template +void test_all() +{ + test_meridian(); + test_meridian(); + test_meridian(); +} + +int test_main(int, char* []) +{ + test_all > >(); + test_all > >(); + test_all > >(); + + return 0; +} From e019d4b728a009a61a4fb2c7af330ae39f609604 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Tue, 31 Oct 2017 16:21:34 +0200 Subject: [PATCH 18/27] [tests] Use meridian distance formula results in andoyer tests --- test/strategies/andoyer.cpp | 32 +++++++++++++++++++++++++------- 1 file changed, 25 insertions(+), 7 deletions(-) diff --git a/test/strategies/andoyer.cpp b/test/strategies/andoyer.cpp index 0b4fff12c..74a7323ec 100644 --- a/test/strategies/andoyer.cpp +++ b/test/strategies/andoyer.cpp @@ -5,9 +5,10 @@ // Copyright (c) 2008-2016 Bruno Lalande, Paris, France. // Copyright (c) 2009-2016 Mateusz Loskot, London, UK. -// This file was modified by Oracle on 2014, 2015, 2016. +// This file was modified by Oracle on 2014-2017. // Modifications copyright (c) 2014-2016 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 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library @@ -94,7 +95,24 @@ void test_distance(double lon1, double lat1, double lon2, double lat2, double ex return_type d_strategy = andoyer.apply(p1, p2); return_type d_function = bg::distance(p1, p2, andoyer); - return_type d_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1), to_rad(lon2), to_rad(lat2), stype()).distance; + + double diff = bg::math::longitude_distance_signed(lon1, lon2); + return_type d_formula; + + // if the points lay on a meridian, distance strategy calls the special formula + // for meridian distance that returns different result than andoyer formula + // for nearly antipodal points + if (bg::math::equals(diff, 0.0) + || bg::math::equals(bg::math::abs(diff), 180.0)) + { + d_formula = d_strategy; + } + else + { + d_formula = andoyer_inverse_type::apply(to_rad(lon1), to_rad(lat1), + to_rad(lon2), to_rad(lat2), + stype()).distance; + } BOOST_CHECK_CLOSE(d_strategy / 1000.0, expected_km, 0.001); BOOST_CHECK_CLOSE(d_function / 1000.0, expected_km, 0.001); @@ -214,10 +232,10 @@ void test_all() // antipodal // ok? in those cases shorter path would pass through a pole // but 90 or -90 would be consistent with distance? - //test_distazi(0, 0, 180, 0, 20037.5, 0.0); - //test_distazi(0, 0, -180, 0, 20037.5, 0.0); - //test_distazi(-90, 0, 90, 0, 20037.5, 0.0); - //test_distazi(90, 0, -90, 0, 20037.5, 0.0); + test_distazi(0, 0, 180, 0, 20003.9, 0.0); + test_distazi(0, 0, -180, 0, 20003.9, 0.0); + test_distazi(-90, 0, 90, 0, 20003.9, 0.0); + test_distazi(90, 0, -90, 0, 20003.9, 0.0); // 0, 45, 90 ... for (int i = 0 ; i < 360 ; i += 45) @@ -264,7 +282,7 @@ void test_all() test_distazi_symm(normlized_deg(l-44.99), -44.99, normlized_deg(l+135), 45, 20008.1, 0.0); test_distazi_symm(normlized_deg(l-44.999), -44.999, normlized_deg(l+135), 45, 20009.4, 0.0); // antipodal - //test_distazi_symm(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20020.7, 0.0); + test_distazi_symm(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20003.92, 0.0); } /* SQL Server gives: From b17ad43f7fbccd5ee985d936b9c9d94cbc3bb252 Mon Sep 17 00:00:00 2001 From: Adam Wulkiewicz Date: Tue, 31 Oct 2017 16:28:41 +0100 Subject: [PATCH 19/27] [strategy] Add missing include in geographic distance strategy. --- include/boost/geometry/strategies/geographic/distance.hpp | 1 + 1 file changed, 1 insertion(+) diff --git a/include/boost/geometry/strategies/geographic/distance.hpp b/include/boost/geometry/strategies/geographic/distance.hpp index bd0e0d5fc..6b195b65e 100644 --- a/include/boost/geometry/strategies/geographic/distance.hpp +++ b/include/boost/geometry/strategies/geographic/distance.hpp @@ -23,6 +23,7 @@ #include #include +#include #include #include From f458d8d28e7273b59375ccccee947428c320fcdb Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 17:42:48 +0200 Subject: [PATCH 20/27] [formulas] [strategies] Distance point-segment use the meridian formula. Use new static version of distance strategy. --- .../formulas/distance_point_segment.hpp | 45 +++++++++---- .../geometry/formulas/elliptic_arc_length.hpp | 50 ++++++++++++++- .../strategies/geographic/distance.hpp | 63 +++++++------------ .../geographic/distance_cross_track.hpp | 3 +- .../algorithms/distance/distance_geo_pl_l.cpp | 8 ++- 5 files changed, 112 insertions(+), 57 deletions(-) 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); - } //=========================================================================== From 5a7ccaff355e016ac2a9a75260687aeea9c3764f Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 19:03:52 +0200 Subject: [PATCH 21/27] [tests] Fix expected results in perimeter tests and add few cases --- test/algorithms/distance/distance_geo_pl_l.cpp | 3 --- test/algorithms/perimeter/perimeter_geo.cpp | 14 ++++++++++++-- .../perimeter/perimeter_polygon_cases.hpp | 4 +++- test/algorithms/test_perimeter.hpp | 4 ++-- 4 files changed, 17 insertions(+), 8 deletions(-) diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index d3bfd0851..16feacb70 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -8,9 +8,6 @@ // 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 diff --git a/test/algorithms/perimeter/perimeter_geo.cpp b/test/algorithms/perimeter/perimeter_geo.cpp index 9b4452dcd..403f87b0b 100644 --- a/test/algorithms/perimeter/perimeter_geo.cpp +++ b/test/algorithms/perimeter/perimeter_geo.cpp @@ -1,7 +1,7 @@ // Boost.Geometry (aka GGL, Generic Geometry Library) // Unit Test -// Copyright (c) 2016 Oracle and/or its affiliates. +// Copyright (c) 2016-2017 Oracle and/or its affiliates. // Contributed and/or modified by Vissarion Fisikopoulos, on behalf of Oracle // Use, modification and distribution is subject to the Boost Software License, @@ -39,9 +39,19 @@ void test_default() //this should use andoyer strategy 1116814.237 + 1116152.605); } + //since the geodesic distance is the shortest path it should go over the pole + //in this case; thus the correct perimeter is the meridian length (below) + //and not 40075160 which is the legth of the equator + test_geometry >(poly_data_geo[3], + 40007834.7139); + + //force to use equator path + test_geometry >(poly_data_geo[4], + 40075016.6856); + // Multipolygon test_geometry > > - (multipoly_data[0], 60078934.0425); + (multipoly_data[0], 60011752.0709); // Geometries with length zero test_geometry

("POINT(0 0)", 0); diff --git a/test/algorithms/perimeter/perimeter_polygon_cases.hpp b/test/algorithms/perimeter/perimeter_polygon_cases.hpp index 68c72a2b1..837829f5f 100644 --- a/test/algorithms/perimeter/perimeter_polygon_cases.hpp +++ b/test/algorithms/perimeter/perimeter_polygon_cases.hpp @@ -17,7 +17,9 @@ static std::string poly_data_geo[] = { "POLYGON((0 90,1 80,1 70))", "POLYGON((0 90,1 80,1 80,1 80,1 70,1 70))", "POLYGON((0 90,1 80,1 79,1 78,1 77,1 76,1 75,1 74,\ - 1 73,1 72,1 71,1 70))" + 1 73,1 72,1 71,1 70))",\ + "POLYGON((0 0,180 0,0 0))", + "POLYGON((0 0,90 0,180 0,90 0,0 0))" }; static std::string poly_data_sph[] = { diff --git a/test/algorithms/test_perimeter.hpp b/test/algorithms/test_perimeter.hpp index b0e703796..32a5defb6 100644 --- a/test/algorithms/test_perimeter.hpp +++ b/test/algorithms/test_perimeter.hpp @@ -29,7 +29,7 @@ void test_perimeter(Geometry const& geometry, long double expected_perimeter) std::ostringstream out; out << typeid(typename bg::coordinate_type::type).name() << std::endl - << typeid(typename bg::default_perimeter_result::type).name() + //<< typeid(typename bg::default_perimeter_result::type).name() << std::endl << "perimeter : " << bg::perimeter(geometry) << std::endl; @@ -50,7 +50,7 @@ void test_perimeter(Geometry const& geometry, long double expected_perimeter, St std::ostringstream out; out << typeid(typename bg::coordinate_type::type).name() << std::endl - << typeid(typename bg::default_perimeter_result::type).name() + //<< typeid(typename bg::default_perimeter_result::type).name() << std::endl << "perimeter : " << bg::perimeter(geometry, strategy) << std::endl; From 1c5d10fb7742ceb9ad28313eb4cae0d1877352ad Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 19:44:28 +0200 Subject: [PATCH 22/27] [formulas] Thomas formula return zero for distance between poles. Fixed by calling meridian formula for this case. --- .../boost/geometry/formulas/elliptic_arc_length.hpp | 12 +++++++----- test/algorithms/distance/distance_geo_pl_l.cpp | 5 +++++ 2 files changed, 12 insertions(+), 5 deletions(-) diff --git a/include/boost/geometry/formulas/elliptic_arc_length.hpp b/include/boost/geometry/formulas/elliptic_arc_length.hpp index e9d7ea7db..06e2b1ad7 100644 --- a/include/boost/geometry/formulas/elliptic_arc_length.hpp +++ b/include/boost/geometry/formulas/elliptic_arc_length.hpp @@ -56,13 +56,15 @@ public : CT half_pi = pi/CT(2); CT diff = math::longitude_distance_signed(lon1, lon2); - if (math::equals(diff, c0)) + if (lat1 > lat2) + { + std::swap(lat1, lat2); + } + + if ( math::equals(diff, c0) || + (math::equals(lat2, half_pi) && math::equals(lat1, -half_pi)) ) { // single meridian not crossing pole - if (lat1 > lat2) - { - std::swap(lat1, lat2); - } res.distance = apply(lat2, spheroid) - apply(lat1, spheroid); res.meridian = true; } diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index 16feacb70..a71ab0aab 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -323,6 +323,11 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, "SEGMENT(0.5 -90,175.5 -90)", pp_distance("POINT(0.5 -90)", "POINT(90 90)", strategy_pp), strategy_ps); + tester::apply("p-s-16", + "POINT(90 90)", + "SEGMENT(0.5 -90,175.5 -90)", + pp_distance("POINT(90 -90)", "POINT(90 90)", strategy_pp), + strategy_ps); } From e099a819e7622f81601d3bb2afb6a0c524d9f6b0 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 19:55:41 +0200 Subject: [PATCH 23/27] [formulas] Clean-up code/comments --- .../formulas/distance_point_segment.hpp | 18 +----------------- 1 file changed, 1 insertion(+), 17 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index c3eafc3a7..8aa93eebc 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -38,8 +38,6 @@ template typename CT, typename Units, typename FormulaPolicy, - //template class Inverse, - //template class Direct, bool EnableClosestPoint = false > class distance_point_segment{ @@ -85,14 +83,6 @@ public: CT lon2, CT lat2, //p2 Spheroid const& spheroid) { - //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); @@ -165,12 +155,6 @@ 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); @@ -182,7 +166,7 @@ public: { #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG std::cout << "Degenerate segment" << std::endl; - std::cout << "d1=" << d1 << std::endl; + std::cout << "distance between points=" << d1 << std::endl; #endif return non_iterative_case(lon1, lat2, d1); } From b1fe2ad934516bd435a135c92ea4989ff9f48e63 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 20:02:49 +0200 Subject: [PATCH 24/27] [formulas] Use new distance function with meridian formula optimization everywhere in point-segment distance formula --- .../boost/geometry/formulas/distance_point_segment.hpp | 9 ++++----- .../strategies/geographic/distance_cross_track.hpp | 2 +- 2 files changed, 5 insertions(+), 6 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 8aa93eebc..04e1994c7 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -45,7 +45,6 @@ class distance_point_segment{ public: 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; @@ -83,7 +82,8 @@ public: CT lon2, CT lat2, //p2 Spheroid const& spheroid) { - CT distance = geometry::strategy::distance::geographic::apply(lon1, lat1, lon2, lat2, spheroid); + CT distance = geometry::strategy::distance::geographic + ::apply(lon1, lat1, lon2, lat2, spheroid); return non_iterative_case(lon1, lat1, distance); } @@ -139,7 +139,6 @@ public: } //segment on equator - //TODO: use the meridian distance when it'll be available if (math::equals(lat1, c0) && math::equals(lat2, c0)) { #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG @@ -171,8 +170,8 @@ public: return non_iterative_case(lon1, lat2, d1); } - CT d2 = inverse_distance_type::apply(lon2, lat2, - lon3, lat3, spheroid).distance; + CT d2 = geometry::strategy::distance::geographic + ::apply(lon2, lat2, lon3, lat3, spheroid); // Compute a12 (GEO) geometry::formula::result_inverse res12 = diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp index 872dad358..7c3e58b19 100644 --- a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -211,7 +211,7 @@ struct get_comparable::type comparable_type; public : static inline comparable_type - apply(geographic_cross_track const& strategy) + apply(geographic_cross_track const& ) { return comparable_type(); } From 8d29c629d01955a752f7d7157e83e97ecb074cb5 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 21:01:35 +0200 Subject: [PATCH 25/27] [test] Special case handling for Pt-Sgmt distance when segment contains the pole --- .../formulas/distance_point_segment.hpp | 8 +++-- .../algorithms/distance/distance_geo_pl_l.cpp | 35 ++++++++++++++++++- 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index 04e1994c7..cd82233bf 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -139,7 +139,11 @@ public: } //segment on equator - if (math::equals(lat1, c0) && math::equals(lat2, c0)) + //Note: antipodal points on equator does not define segment on equator + //but pass by the pole + CT diff = math::longitude_distance_signed(lon1, lon2); + if (math::equals(lat1, c0) && math::equals(lat2, c0) + && !math::equals(math::abs(diff), pi)) { #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG std::cout << "Equatorial segment" << std::endl; @@ -220,7 +224,7 @@ public: #ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG std::cout << "a21=" << a21 * math::r2d() << std::endl; - std::cout << "a23=" << a13 * math::r2d() << std::endl; + std::cout << "a23=" << a23 * math::r2d() << std::endl; std::cout << "a321=" << a321 * math::r2d() << std::endl; std::cout << "cos(a321)=" << cos(a321) << std::endl; #endif diff --git a/test/algorithms/distance/distance_geo_pl_l.cpp b/test/algorithms/distance/distance_geo_pl_l.cpp index a71ab0aab..fd5c7518a 100644 --- a/test/algorithms/distance/distance_geo_pl_l.cpp +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -7,7 +7,8 @@ // 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 @@ -328,7 +329,39 @@ void test_distance_point_segment(Strategy_pp const& strategy_pp, "SEGMENT(0.5 -90,175.5 -90)", pp_distance("POINT(90 -90)", "POINT(90 90)", strategy_pp), strategy_ps); + // equatorial segment + tester::apply("p-s-17", + "POINT(90 90)", + "SEGMENT(0 0,175.5 0)", + pp_distance("POINT(90 90)", "POINT(0 0)", strategy_pp), + strategy_ps); + // segment pass by pole + tester::apply("p-s-18", + "POINT(90 90)", + "SEGMENT(0 0,180 0)", + 0, + strategy_ps); +/* The following fails for thomas strategy + tester::apply("p-s-19", + "POINT(1 80)", + "SEGMENT(0 0,0 89)", + pp_distance("POINT(1 80)", "POINT(0 80.00149297334342)", strategy_pp), + strategy_ps); +*/ +// Tests failing, convergence issues of the method when dealing with poles +/* + tester::apply("p-s-20", + "POINT(1 80)", + "SEGMENT(0 0,0 90)", + pp_distance("POINT(1 80)", "POINT(0 80.00149297334342)", strategy_pp), + strategy_ps); + tester::apply("p-s-21", + "POINT(90 89)", + "SEGMENT(0 0,180 0)", + pp_distance("POINT(90 89)", "POINT(90 90)", strategy_pp), + strategy_ps); +*/ } //=========================================================================== From 730d91c3b78b42bf99a4956609cc9542df621822 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 22:05:44 +0200 Subject: [PATCH 26/27] [formulas] Missing include files and namespaces --- include/boost/geometry/formulas/distance_point_segment.hpp | 4 +++- include/boost/geometry/formulas/elliptic_arc_length.hpp | 4 ++-- 2 files changed, 5 insertions(+), 3 deletions(-) diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp index cd82233bf..fd003a778 100644 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ b/include/boost/geometry/formulas/distance_point_segment.hpp @@ -14,6 +14,8 @@ #include #include +#include + #ifndef BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS #define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100 #endif @@ -141,7 +143,7 @@ public: //segment on equator //Note: antipodal points on equator does not define segment on equator //but pass by the pole - CT diff = math::longitude_distance_signed(lon1, lon2); + CT diff = geometry::math::longitude_distance_signed(lon1, lon2); if (math::equals(lat1, c0) && math::equals(lat2, c0) && !math::equals(math::abs(diff), pi)) { diff --git a/include/boost/geometry/formulas/elliptic_arc_length.hpp b/include/boost/geometry/formulas/elliptic_arc_length.hpp index 06e2b1ad7..75881f3d0 100644 --- a/include/boost/geometry/formulas/elliptic_arc_length.hpp +++ b/include/boost/geometry/formulas/elliptic_arc_length.hpp @@ -11,7 +11,6 @@ #ifndef BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP #define BOOST_GEOMETRY_FORMULAS_ELLIPTIC_ARC_LENGTH_HPP - #include #include @@ -19,6 +18,7 @@ #include #include +#include #include @@ -54,7 +54,7 @@ public : CT c0 = 0; CT pi = math::pi(); CT half_pi = pi/CT(2); - CT diff = math::longitude_distance_signed(lon1, lon2); + CT diff = geometry::math::longitude_distance_signed(lon1, lon2); if (lat1 > lat2) { From 8d412963c8c27dec63b020fbe2a196221a2ecdf1 Mon Sep 17 00:00:00 2001 From: Vissarion Fysikopoulos Date: Wed, 1 Nov 2017 23:41:36 +0200 Subject: [PATCH 27/27] [formulas] [strategies] Move point-segment distance formula to strategies --- .../formulas/distance_point_segment.hpp | 387 ------------------ .../geographic/distance_cross_track.hpp | 381 ++++++++++++++++- 2 files changed, 365 insertions(+), 403 deletions(-) delete mode 100644 include/boost/geometry/formulas/distance_point_segment.hpp diff --git a/include/boost/geometry/formulas/distance_point_segment.hpp b/include/boost/geometry/formulas/distance_point_segment.hpp deleted file mode 100644 index fd003a778..000000000 --- a/include/boost/geometry/formulas/distance_point_segment.hpp +++ /dev/null @@ -1,387 +0,0 @@ -// Boost.Geometry - -// Copyright (c) 2016-2017 Oracle and/or its affiliates. - -// Contributed and/or modified by Vissarion Fysikopoulos, 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_FORMULAS_CROSS_TRACK_GEO_HPP -#define BOOST_GEOMETRY_FORMULAS_CROSS_TRACK_GEO_HPP - -#include -#include - -#include - -#ifndef BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS -#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 - approximates the distance by an iterative Newton method. -\see C.F.F.Karney - Geodesics on an ellipsoid of revolution, - https://arxiv.org/abs/1102.1215 -*/ - -namespace boost { namespace geometry { namespace formula -{ - -template -< - typename CT, - typename Units, - typename FormulaPolicy, - bool EnableClosestPoint = false -> -class distance_point_segment{ - -public: - - typedef typename FormulaPolicy::template inverse inverse_distance_quantities_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 - { - result_distance_point_segment() - : distance(0) - , closest_point_lon(0) - , closest_point_lat(0) - {} - - CT distance; - CT closest_point_lon; - CT closest_point_lat; - }; - - result_distance_point_segment - static inline non_iterative_case(CT lon, CT lat, CT distance) - { - result_distance_point_segment result; - result.distance = distance; - - if (EnableClosestPoint) - { - result.closest_point_lon = lon; - result.closest_point_lat = lat; - } - return result; - } - - template - result_distance_point_segment - static inline non_iterative_case(CT lon1, CT lat1, //p1 - CT lon2, CT lat2, //p2 - Spheroid const& spheroid) - { - CT distance = geometry::strategy::distance::geographic - ::apply(lon1, lat1, lon2, lat2, spheroid); - - return non_iterative_case(lon1, lat1, distance); - } - - template - T static inline normalize(T g4) - { - CT const pi = math::pi(); - if (g4 < 0 && g4 < -pi)//close to -270 - { - return g4 + 1.5 * pi; - } - else if (g4 > 0 && g4 > pi)//close to 270 - { - return - g4 + 1.5 * pi; - } - else if (g4 < 0 && g4 > -pi)//close to -90 - { - return -g4 - pi/2; - } - return g4 - pi/2; - } - - template - result_distance_point_segment - static inline apply(CT lon1, CT lat1, //p1 - CT lon2, CT lat2, //p2 - CT lon3, CT lat3, //query point p3 - Spheroid const& spheroid) - { - CT const earth_radius = geometry::formula::mean_radius(spheroid); - - result_distance_point_segment result; - - // Constants - CT const f = flattening(spheroid); - CT const pi = math::pi(); - CT const half_pi = pi / CT(2); - CT const c0 = CT(0); - - // Convert to radians - lon1 = math::as_radian(lon1); - lat1 = math::as_radian(lat1); - lon2 = math::as_radian(lon2); - lat2 = math::as_radian(lat2); - lon3 = math::as_radian(lon3); - lat3 = math::as_radian(lat3); - - if (lon1 > lon2) - { - std::swap(lon1, lon2); - std::swap(lat1, lat2); - } - - //segment on equator - //Note: antipodal points on equator does not define segment on equator - //but pass by the pole - CT diff = geometry::math::longitude_distance_signed(lon1, lon2); - if (math::equals(lat1, c0) && math::equals(lat2, c0) - && !math::equals(math::abs(diff), pi)) - { -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "Equatorial segment" << std::endl; -#endif - if (lon3 <= lon1) - { - return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); - } - if (lon3 >= lon2) - { - return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); - } - return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); - } - - 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 << "distance between points=" << d1 << std::endl; -#endif - return non_iterative_case(lon1, lat2, d1); - } - - CT d2 = geometry::strategy::distance::geographic - ::apply(lon2, lat2, lon3, lat3, spheroid); - - // Compute a12 (GEO) - geometry::formula::result_inverse 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; - - if (geometry::math::equals(a312, c0)) - { -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "point on segment" << std::endl; -#endif - return non_iterative_case(lon3, lat3, c0); - } - - CT projection1 = cos( a312 ) * d1 / d3; - -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "segment=(" << lon1 * math::r2d(); - std::cout << "," << lat1 * math::r2d(); - std::cout << "),(" << lon2 * math::r2d(); - std::cout << "," << lat2 * math::r2d(); - std::cout << ")\np=(" << lon3 * math::r2d(); - std::cout << "," << lat3 * math::r2d(); - std::cout << ")\na1=" << a12 * math::r2d() << std::endl; - std::cout << "a13=" << a13 * math::r2d() << std::endl; - std::cout << "a312=" << a312 * math::r2d() << std::endl; - std::cout << "cos(a312)=" << cos(a312) << std::endl; -#endif - if (projection1 < 0.0) - { -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "projection closer to p1" << std::endl; -#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; - -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "a21=" << a21 * math::r2d() << std::endl; - std::cout << "a23=" << a23 * math::r2d() << std::endl; - std::cout << "a321=" << a321 * math::r2d() << std::endl; - std::cout << "cos(a321)=" << cos(a321) << std::endl; -#endif - CT projection2 = cos( a321 ) * d2 / d3; - - if (projection2 < 0.0) - { -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "projection closer to p2" << std::endl; -#endif - // projection of p3 on geodesic spanned by segment (p1,p2) fall - // outside of segment on the side of p2 - return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); - } - - // Guess s14 (SPHERICAL) - typedef geometry::model::point - < - CT, 2, - geometry::cs::spherical_equatorial - > point; - - CT bet1 = atan((1 - f) * tan(lon1)); - CT bet2 = atan((1 - f) * tan(lon2)); - CT bet3 = atan((1 - f) * tan(lon3)); - point p1 = point(bet1, lat1); - point p2 = point(bet2, lat2); - point p3 = point(bet3, lat3); - - geometry::strategy::distance::cross_track cross_track(earth_radius); - CT s34 = cross_track.apply(p3, p1, p2); - - geometry::strategy::distance::haversine str(earth_radius); - CT s13 = str.apply(p1, p3); - CT s14 = acos( cos(s13/earth_radius) / cos(s34/earth_radius) ) * earth_radius; - -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - 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 = 0; - geometry::formula::result_direct res14; - geometry::formula::result_inverse res34; - - int counter = 0; // robustness - CT g4; - CT delta_g4; - - 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_azimuth_type::apply(res14.lon2, res14.lat2, - lon2, lat2, spheroid).azimuth; - res34 = inverse_distance_quantities_type::apply(res14.lon2, res14.lat2, - lon3, lat3, spheroid); - g4 = res34.azimuth - a4; - - delta_g4 = normalize(g4); - - CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit - CT m34 = res34.reduced_length; - CT der = (M43 / m34) * sin(g4); - s14 = s14 - delta_g4 / der; - -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - std::cout << "p4=" << res14.lon2 * math::r2d() << - "," << res14.lat2 * math::r2d() << std::endl; - std::cout << "delta_g4=" << delta_g4 << std::endl; - std::cout << "g4=" << g4 * math::r2d() << 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; -#endif - result.distance = prev_distance; - -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - if (g4 == half_pi) - { - std::cout << "Stop msg: g4 == half_pi" << std::endl; - } - if (res34.distance >= prev_distance && prev_distance != 0) - { - 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 == 19) - { - std::cout << "Stop msg: counter" << std::endl; - } -#endif - - } while (g4 != half_pi - && (prev_distance > res34.distance || prev_distance == 0) - && delta_g4 != 0 - && ++counter < BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS ) ; - -#ifdef BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG - 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_quantities_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance - << ", p4=(" << get<0>(p4) * math::r2d() << "," - << get<1>(p4) * math::r2d() << ")" - << std::endl; - - CT s31 = inverse_distance_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; - CT s32 = inverse_distance_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 res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); - CT p4_plus = inverse_distance_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; - - geometry::formula::result_direct res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a12, spheroid); - CT p4_minus = inverse_distance_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() << "," << res4.lat2 * math::r2d() << ")" - << "\np4_minus=" << p4_minus << ", p4=(" << res1.lon2 * math::r2d() << "," << res1.lat2 * math::r2d() << ")" - << 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 - - return result; - } - -}; - -}}} // namespace boost::geometry::formula -#endif // BOOST_GEOMETRY_FORMULAS_CROSS_TRACK_GEO_HPP diff --git a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp index 7c3e58b19..1f8fbc41e 100644 --- a/include/boost/geometry/strategies/geographic/distance_cross_track.hpp +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -30,16 +30,26 @@ #include #include -#include #include #include #include +#include -#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK +#include +#include + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK # include #endif +#ifndef BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS +#define BOOST_GEOMETRY_DETAIL_POINT_SEGMENT_DISTANCE_MAX_STEPS 100 +#endif + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK +#include +#endif namespace boost { namespace geometry { @@ -49,19 +59,24 @@ namespace strategy { namespace distance /*! \brief Strategy functor for distance point to segment calculation on ellipsoid + Algorithm uses direct and inverse geodesic problems as subroutines. + The algorithm approximates the distance by an iterative Newton method. \ingroup strategies \details Class which calculates the distance of a point to a segment, for points on the ellipsoid -\see https://arxiv.org/abs/1102.1215 +\see C.F.F.Karney - Geodesics on an ellipsoid of revolution, + https://arxiv.org/abs/1102.1215 \tparam FormulaPolicy underlying point-point distance strategy \tparam Spheroid is the spheroidal model used \tparam CalculationType \tparam_calculation +\tparam EnableClosestPoint computes the closest point on segment if true */ template < typename FormulaPolicy = strategy::andoyer, typename Spheroid = srs::spheroid, - typename CalculationType = void + typename CalculationType = void, + bool EnableClosestPoint = false > class geographic_cross_track { @@ -100,22 +115,356 @@ public : { typedef typename coordinate_system::type::units units_type; - typedef typename return_type::type CT; - - return (geometry::formula::distance_point_segment - < - CT, - units_type, - FormulaPolicy - >::apply(get<0>(sp1), get<1>(sp1), - get<0>(sp2), get<1>(sp2), - get<0>(p), get<1>(p), - m_spheroid)).distance; + return (apply(get<0>(sp1), get<1>(sp1), + get<0>(sp2), get<1>(sp2), + get<0>(p), get<1>(p), + m_spheroid)).distance; } private : - Spheroid m_spheroid; + template + struct result_distance_point_segment + { + result_distance_point_segment() + : distance(0) + , closest_point_lon(0) + , closest_point_lat(0) + {} + + CT distance; + CT closest_point_lon; + CT closest_point_lat; + }; + + template + result_distance_point_segment + static inline non_iterative_case(CT lon, CT lat, CT distance) + { + result_distance_point_segment result; + result.distance = distance; + + if (EnableClosestPoint) + { + result.closest_point_lon = lon; + result.closest_point_lat = lat; + } + return result; + } + + template + result_distance_point_segment + static inline non_iterative_case(CT lon1, CT lat1, //p1 + CT lon2, CT lat2, //p2 + Spheroid const& spheroid) + { + CT distance = geometry::strategy::distance::geographic + ::apply(lon1, lat1, lon2, lat2, spheroid); + + return non_iterative_case(lon1, lat1, distance); + } + + template + CT static inline normalize(CT g4) + { + CT const pi = math::pi(); + if (g4 < 0 && g4 < -pi)//close to -270 + { + return g4 + 1.5 * pi; + } + else if (g4 > 0 && g4 > pi)//close to 270 + { + return - g4 + 1.5 * pi; + } + else if (g4 < 0 && g4 > -pi)//close to -90 + { + return -g4 - pi/2; + } + return g4 - pi/2; + } + + template + result_distance_point_segment + static inline apply(CT lon1, CT lat1, //p1 + CT lon2, CT lat2, //p2 + CT lon3, CT lat3, //query point p3 + Spheroid const& spheroid) + { + typedef typename FormulaPolicy::template inverse + inverse_distance_quantities_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; + + CT const earth_radius = geometry::formula::mean_radius(spheroid); + + result_distance_point_segment result; + + // Constants + CT const f = geometry::formula::flattening(spheroid); + CT const pi = math::pi(); + CT const half_pi = pi / CT(2); + CT const c0 = CT(0); + + // Convert to radians + lon1 = math::as_radian(lon1); + lat1 = math::as_radian(lat1); + lon2 = math::as_radian(lon2); + lat2 = math::as_radian(lat2); + lon3 = math::as_radian(lon3); + lat3 = math::as_radian(lat3); + + if (lon1 > lon2) + { + std::swap(lon1, lon2); + std::swap(lat1, lat2); + } + + //segment on equator + //Note: antipodal points on equator does not define segment on equator + //but pass by the pole + CT diff = geometry::math::longitude_distance_signed(lon1, lon2); + if (math::equals(lat1, c0) && math::equals(lat2, c0) + && !math::equals(math::abs(diff), pi)) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "Equatorial segment" << std::endl; +#endif + if (lon3 <= lon1) + { + return non_iterative_case(lon1, lat1, lon3, lat3, spheroid); + } + if (lon3 >= lon2) + { + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); + } + return non_iterative_case(lon3, lat1, lon3, lat3, spheroid); + } + + 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_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "Degenerate segment" << std::endl; + std::cout << "distance between points=" << d1 << std::endl; +#endif + return non_iterative_case(lon1, lat2, d1); + } + + CT d2 = geometry::strategy::distance::geographic + ::apply(lon2, lat2, lon3, lat3, spheroid); + + // Compute a12 (GEO) + geometry::formula::result_inverse 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; + + if (geometry::math::equals(a312, c0)) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "point on segment" << std::endl; +#endif + return non_iterative_case(lon3, lat3, c0); + } + + CT projection1 = cos( a312 ) * d1 / d3; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "segment=(" << lon1 * math::r2d(); + std::cout << "," << lat1 * math::r2d(); + std::cout << "),(" << lon2 * math::r2d(); + std::cout << "," << lat2 * math::r2d(); + std::cout << ")\np=(" << lon3 * math::r2d(); + std::cout << "," << lat3 * math::r2d(); + std::cout << ")\na1=" << a12 * math::r2d() << std::endl; + std::cout << "a13=" << a13 * math::r2d() << std::endl; + std::cout << "a312=" << a312 * math::r2d() << std::endl; + std::cout << "cos(a312)=" << cos(a312) << std::endl; +#endif + if (projection1 < 0.0) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "projection closer to p1" << std::endl; +#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; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "a21=" << a21 * math::r2d() << std::endl; + std::cout << "a23=" << a23 * math::r2d() << std::endl; + std::cout << "a321=" << a321 * math::r2d() << std::endl; + std::cout << "cos(a321)=" << cos(a321) << std::endl; +#endif + CT projection2 = cos( a321 ) * d2 / d3; + + if (projection2 < 0.0) + { +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "projection closer to p2" << std::endl; +#endif + // projection of p3 on geodesic spanned by segment (p1,p2) fall + // outside of segment on the side of p2 + return non_iterative_case(lon2, lat2, lon3, lat3, spheroid); + } + + // Guess s14 (SPHERICAL) + typedef geometry::model::point + < + CT, 2, + geometry::cs::spherical_equatorial + > point; + + CT bet1 = atan((1 - f) * tan(lon1)); + CT bet2 = atan((1 - f) * tan(lon2)); + CT bet3 = atan((1 - f) * tan(lon3)); + point p1 = point(bet1, lat1); + point p2 = point(bet2, lat2); + point p3 = point(bet3, lat3); + + geometry::strategy::distance::cross_track cross_track(earth_radius); + CT s34 = cross_track.apply(p3, p1, p2); + + geometry::strategy::distance::haversine str(earth_radius); + CT s13 = str.apply(p1, p3); + CT s14 = acos( cos(s13/earth_radius) / cos(s34/earth_radius) ) * 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 = 0; + geometry::formula::result_direct res14; + geometry::formula::result_inverse res34; + + int counter = 0; // robustness + CT g4; + CT delta_g4; + + 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_azimuth_type::apply(res14.lon2, res14.lat2, + lon2, lat2, spheroid).azimuth; + res34 = inverse_distance_quantities_type::apply(res14.lon2, res14.lat2, + lon3, lat3, spheroid); + g4 = res34.azimuth - a4; + + delta_g4 = normalize(g4); + + CT M43 = res34.geodesic_scale; // cos(s14/earth_radius) is the spherical limit + CT m34 = res34.reduced_length; + CT der = (M43 / m34) * sin(g4); + s14 = s14 - delta_g4 / der; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + std::cout << "p4=" << res14.lon2 * math::r2d() << + "," << res14.lat2 * math::r2d() << std::endl; + std::cout << "delta_g4=" << delta_g4 << std::endl; + std::cout << "g4=" << g4 * math::r2d() << 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; +#endif + result.distance = prev_distance; + +#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK + if (g4 == half_pi) + { + std::cout << "Stop msg: g4 == half_pi" << std::endl; + } + if (res34.distance >= prev_distance && prev_distance != 0) + { + 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 == 19) + { + std::cout << "Stop msg: counter" << std::endl; + } +#endif + + } while (g4 != half_pi + && (prev_distance > res34.distance || prev_distance == 0) + && 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_quantities_type::apply(get<0>(p4), get<1>(p4), lon3, lat3, spheroid).distance + << ", p4=(" << get<0>(p4) * math::r2d() << "," + << get<1>(p4) * math::r2d() << ")" + << std::endl; + + CT s31 = inverse_distance_quantities_type::apply(lon3, lat3, lon1, lat1, spheroid).distance; + CT s32 = inverse_distance_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 res4 = direct_distance_type::apply(get<0>(p4), get<1>(p4), .04, a4, spheroid); + CT p4_plus = inverse_distance_quantities_type::apply(res4.lon2, res4.lat2, lon3, lat3, spheroid).distance; + + geometry::formula::result_direct res1 = direct_distance_type::apply(lon1, lat1, s14-.04, a12, spheroid); + CT p4_minus = inverse_distance_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() << "," << res4.lat2 * math::r2d() << ")" + << "\np4_minus=" << p4_minus << ", p4=(" << res1.lon2 * math::r2d() << "," << res1.lat2 * math::r2d() << ")" + << 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 + + return result; + } + + Spheroid m_spheroid; };