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/andoyer_inverse.hpp b/include/boost/geometry/formulas/andoyer_inverse.hpp index 902fd7d8f..5c195dbdb 100644 --- a/include/boost/geometry/formulas/andoyer_inverse.hpp +++ b/include/boost/geometry/formulas/andoyer_inverse.hpp @@ -214,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 new file mode 100644 index 000000000..75881f3d0 --- /dev/null +++ b/include/boost/geometry/formulas/elliptic_arc_length.hpp @@ -0,0 +1,214 @@ +// 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 + +#include + +namespace boost { namespace geometry { namespace formula +{ + +/*! +\brief Compute the arc length of an ellipse. +*/ + +template +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 = geometry::math::longitude_distance_signed(lon1, lon2); + + 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 + 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 + // 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); + 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 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)); + + } + + // 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 and in [-pi/2,pi/2] + 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 diff --git a/include/boost/geometry/strategies/geographic/distance.hpp b/include/boost/geometry/strategies/geographic/distance.hpp index d3656f449..98ee46a36 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, @@ -22,12 +23,14 @@ #include #include +#include #include #include #include #include +#include #include #include @@ -70,17 +73,41 @@ 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 { - 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; + 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); + + 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 new file mode 100644 index 000000000..1f8fbc41e --- /dev/null +++ b/include/boost/geometry/strategies/geographic/distance_cross_track.hpp @@ -0,0 +1,624 @@ +// 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 +#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 +{ + +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 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, + bool EnableClosestPoint = false +> +class geographic_cross_track +{ +public : + template + struct return_type + : promote_floating_point + < + typename select_calculation_type + < + Point, + PointOfSegment, + CalculationType + >::type + > + {}; + + 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 geographic_cross_track(Spheroid const& spheroid = Spheroid()) + : m_spheroid(spheroid) + {} + + 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; + + return (apply(get<0>(sp1), get<1>(sp1), + get<0>(sp2), get<1>(sp2), + get<0>(p), get<1>(p), + m_spheroid)).distance; + } + +private : + + 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; +}; + + + +#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> + : geographic_cross_track::template return_type +{}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename P, + typename PS +> +struct return_type, P, PS> + : geographic_cross_track::template return_type +{}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType, + typename P, + typename PS +> +struct return_type, P, PS> + : geographic_cross_track::template return_type +{}; + +//comparable types +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct comparable_type > +{ + typedef geographic_cross_track + < + FormulaPolicy, Spheroid, CalculationType + > type; +}; + +template +< + typename FormulaPolicy, + typename Spheroid, + typename CalculationType +> +struct get_comparable > +{ + typedef typename comparable_type + < + geographic_cross_track + >::type comparable_type; +public : + static inline comparable_type + apply(geographic_cross_track const& ) + { + return comparable_type(); + } +}; + + +template +< + typename FormulaPolicy, + typename P, + typename PS +> +struct result_from_distance, P, PS> +{ +private : + typedef typename geographic_cross_track + < + FormulaPolicy + >::template return_type::type return_type; +public : + template + static inline return_type + apply(geographic_cross_track const& , T const& distance) + { + return distance; + } +}; + + +template +struct default_strategy + < + point_tag, segment_tag, Point, PointOfSegment, + geographic_tag, geographic_tag + > +{ + typedef geographic_cross_track<> 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..4896e42e8 100644 --- a/include/boost/geometry/strategies/geographic/parameters.hpp +++ b/include/boost/geometry/strategies/geographic/parameters.hpp @@ -12,7 +12,9 @@ #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/spherical/distance_cross_track.hpp b/include/boost/geometry/strategies/spherical/distance_cross_track.hpp index 7daafa4a1..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 @@ -373,19 +372,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 +384,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 +410,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/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp b/include/boost/geometry/strategies/spherical/side_by_cross_track.hpp index 3f7be0555..2b195e96f 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,20 @@ 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_AB = geometry::formula::spherical_azimuth + (lon1, lat1, lon2, lat2).azimuth; + calc_t XTD = asin(sin(d1) * sin(crs_AD - crs_AB)); return math::equals(XTD, 0) ? 0 : XTD < 0 ? 1 : -1; diff --git a/include/boost/geometry/strategies/strategies.hpp b/include/boost/geometry/strategies/strategies.hpp index ef7793cae..e7f3604ab 100644 --- a/include/boost/geometry/strategies/strategies.hpp +++ b/include/boost/geometry/strategies/strategies.hpp @@ -85,6 +85,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..fd5c7518a --- /dev/null +++ b/test/algorithms/distance/distance_geo_pl_l.cpp @@ -0,0 +1,682 @@ +// 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 +#define BOOST_GEOMETRY_TEST_DEBUG +#define BOOST_GEOMETRY_DISTANCE_POINT_SEGMENT_DEBUG +#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::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::geographic_cross_track + andoyer_strategy; + +typedef bg::strategy::distance::geographic_cross_track + thomas_strategy; + +typedef bg::strategy::distance::geographic_cross_track + 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_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 +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-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)", + 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.000114792872075)", andoyer_pp()), + andoyer_strategy()); + 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); + 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()); + + //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)", + "SEGMENT(0.5 0,175.5 0)", + 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)", + 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)", + 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)", + 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)", + 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)", + 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-15", + "POINT(90 90)", + "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); + // 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); +*/ +} + +//=========================================================================== + +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()); +} 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)", 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..f1a06e0a2 --- /dev/null +++ b/test/algorithms/distance/test_distance_geo_common.hpp @@ -0,0 +1,277 @@ +// 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 "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); + } +}; + +//======================================================================== + +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 + 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, + strategy, test_reversed); + } + + + template + < + typename DistanceType, + typename Strategy + > + static inline + void apply(std::string const& case_id, + Geometry1 const& geometry1, + Geometry2 const& geometry2, + DistanceType const& expected_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); + + // 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); + +#ifdef BOOST_GEOMETRY_TEST_DEBUG + std::cout << string_from_type::type>::name() + << string_from_type::type>::name() + << " -> " + << string_from_type::name() + << std::endl; + std::cout << "expected distance = " + << expected_distance << " ; " + << std::endl; + std::cout << "distance = " + << dist << " ; " + << 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); + +#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 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; diff --git a/test/strategies/andoyer.cpp b/test/strategies/andoyer.cpp index 01b31c73d..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: 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; +}