Merge branch 'develop' into bg-prepare

include/boost/geometry/algorithms/detail/distance/point_to_geometry.hpp

@@ -391,7 +391,7 @@ struct distance
 template <typename Point, typename Polygon, typename Strategy>
 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
         <

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

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 <boost/math/constants/constants.hpp>
+
+#include <boost/geometry/core/radius.hpp>
+#include <boost/geometry/core/srs.hpp>
+
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
+
+#include <boost/geometry/formulas/flattening.hpp>
+
+namespace boost { namespace geometry { namespace formula
+{
+
+/*!
+\brief Compute the arc length of an ellipse.
+*/
+
+template <typename CT, unsigned int Order = 1>
+class elliptic_arc_length
+{
+
+public :
+
+    struct result
+    {
+        result()
+            : distance(0)
+            , meridian(false)
+        {}
+
+        CT distance;
+        bool meridian;
+    };
+
+    template <typename T, typename Spheroid>
+    static result apply(T lon1, T lat1, T lon2, T lat2, Spheroid const& spheroid)
+    {
+        result res;
+
+        CT c0 = 0;
+        CT pi = math::pi<CT>();
+        CT half_pi = pi/CT(2);
+        CT diff = geometry::math::longitude_distance_signed<geometry::radian>(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 <typename T, typename Spheroid>
+    static CT apply(T lat, Spheroid const& spheroid)
+    {
+        CT const a = get_radius<0>(spheroid);
+        CT const f = formula::flattening<CT>(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 <typename T1, typename T2, typename Spheroid>
+    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<CT>(spheroid);
+
+        // how many steps to use
+
+        CT lat1_deg = lat1 * geometry::math::r2d<CT>();
+        CT lat2_deg = lat2 * geometry::math::r2d<CT>();
+
+        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

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 <boost/geometry/core/srs.hpp>
 
 #include <boost/geometry/formulas/andoyer_inverse.hpp>
+#include <boost/geometry/formulas/elliptic_arc_length.hpp>
 #include <boost/geometry/formulas/flattening.hpp>
 
 #include <boost/geometry/strategies/distance.hpp>
 #include <boost/geometry/strategies/geographic/parameters.hpp>
 
 #include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
 #include <boost/geometry/util/promote_floating_point.hpp>
 #include <boost/geometry/util/select_calculation_type.hpp>
 
@@ -70,17 +73,41 @@ public :
         : m_spheroid(spheroid)
     {}
 
+    template <typename CT>
+    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<FormulaPolicy>::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 <typename Point1, typename Point2>
     inline typename calculation_type<Point1, Point2>::type
     apply(Point1 const& point1, Point2 const& point2) const
     {
-        return FormulaPolicy::template inverse
-            <
-                typename calculation_type<Point1, Point2>::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<Point1, Point2>::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

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 <algorithm>
+
+#include <boost/config.hpp>
+#include <boost/concept_check.hpp>
+#include <boost/mpl/if.hpp>
+#include <boost/type_traits/is_void.hpp>
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/radian_access.hpp>
+#include <boost/geometry/core/tags.hpp>
+
+#include <boost/geometry/strategies/distance.hpp>
+#include <boost/geometry/strategies/concepts/distance_concept.hpp>
+#include <boost/geometry/strategies/spherical/distance_haversine.hpp>
+#include <boost/geometry/strategies/geographic/azimuth.hpp>
+#include <boost/geometry/strategies/geographic/parameters.hpp>
+
+#include <boost/geometry/formulas/vincenty_direct.hpp>
+
+#include <boost/geometry/util/math.hpp>
+#include <boost/geometry/util/promote_floating_point.hpp>
+#include <boost/geometry/util/select_calculation_type.hpp>
+#include <boost/geometry/util/normalize_spheroidal_coordinates.hpp>
+
+#include <boost/geometry/formulas/result_direct.hpp>
+#include <boost/geometry/formulas/mean_radius.hpp>
+
+#ifdef BOOST_GEOMETRY_DEBUG_GEOGRAPHIC_CROSS_TRACK
+#  include <boost/geometry/io/dsv/write.hpp>
+#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 <iostream>
+#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<double>,
+    typename CalculationType = void,
+    bool EnableClosestPoint = false
+>
+class geographic_cross_track
+{
+public :
+    template <typename Point, typename PointOfSegment>
+    struct return_type
+        : promote_floating_point
+          <
+              typename select_calculation_type
+                  <
+                      Point,
+                      PointOfSegment,
+                      CalculationType
+                  >::type
+          >
+    {};
+
+    struct distance_strategy
+    {
+        typedef geographic<FormulaPolicy, Spheroid, CalculationType> 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 <typename Point, typename PointOfSegment>
+    inline typename return_type<Point, PointOfSegment>::type
+    apply(Point const& p, PointOfSegment const& sp1, PointOfSegment const& sp2) const
+    {
+        typedef typename coordinate_system<Point>::type::units units_type;
+
+        return (apply<units_type>(get<0>(sp1), get<1>(sp1),
+                                  get<0>(sp2), get<1>(sp2),
+                                  get<0>(p), get<1>(p),
+                                  m_spheroid)).distance;
+    }
+
+private :
+
+    template <typename CT>
+    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 <typename CT>
+    result_distance_point_segment<CT>
+    static inline non_iterative_case(CT lon, CT lat, CT distance)
+    {
+        result_distance_point_segment<CT> result;
+        result.distance = distance;
+
+        if (EnableClosestPoint)
+        {
+            result.closest_point_lon = lon;
+            result.closest_point_lat = lat;
+        }
+        return result;
+    }
+
+    template <typename CT>
+    result_distance_point_segment<CT>
+    static inline non_iterative_case(CT lon1, CT lat1, //p1
+                                     CT lon2, CT lat2, //p2
+                                     Spheroid const& spheroid)
+    {
+        CT distance = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT>
+                              ::apply(lon1, lat1, lon2, lat2, spheroid);
+
+        return non_iterative_case(lon1, lat1, distance);
+    }
+
+    template <typename CT>
+    CT static inline normalize(CT g4)
+    {
+        CT const pi = math::pi<CT>();
+        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 <typename Units, typename CT>
+    result_distance_point_segment<CT>
+    static inline apply(CT lon1, CT lat1, //p1
+                        CT lon2, CT lat2, //p2
+                        CT lon3, CT lat3, //query point p3
+                        Spheroid const& spheroid)
+    {
+        typedef typename FormulaPolicy::template inverse<CT, true, false, false, true, true>
+                inverse_distance_quantities_type;
+        typedef typename FormulaPolicy::template inverse<CT, false, true, false, false, false>
+                inverse_azimuth_type;
+        typedef typename FormulaPolicy::template inverse<CT, false, true, true, false, false>
+                inverse_azimuth_reverse_type;
+        typedef typename FormulaPolicy::template direct<CT, true, false, false, false>
+                direct_distance_type;
+
+        CT const earth_radius = geometry::formula::mean_radius<CT>(spheroid);
+
+        result_distance_point_segment<CT> result;
+
+        // Constants
+        CT const f = geometry::formula::flattening<CT>(spheroid);
+        CT const pi = math::pi<CT>();
+        CT const half_pi = pi / CT(2);
+        CT const c0 = CT(0);
+
+        // Convert to radians
+        lon1 = math::as_radian<Units>(lon1);
+        lat1 = math::as_radian<Units>(lat1);
+        lon2 = math::as_radian<Units>(lon2);
+        lat2 = math::as_radian<Units>(lat2);
+        lon3 = math::as_radian<Units>(lon3);
+        lat3 = math::as_radian<Units>(lat3);
+
+        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<geometry::radian>(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<FormulaPolicy, Spheroid, CT>
+                ::apply(lon1, lat1, lon3, lat3, spheroid);
+
+        CT d3 = geometry::strategy::distance::geographic<FormulaPolicy, Spheroid, CT>
+                ::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<FormulaPolicy, Spheroid, CT>
+                ::apply(lon2, lat2, lon3, lat3, spheroid);
+
+        // Compute a12 (GEO)
+        geometry::formula::result_inverse<CT> res12 =
+                inverse_azimuth_reverse_type::apply(lon1, lat1, lon2, lat2, spheroid);
+        CT a12 = res12.azimuth;
+        CT a13 = inverse_azimuth_type::apply(lon1, lat1, lon3, lat3, spheroid).azimuth;
+
+        CT a312 = a13 - a12;
+
+        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<CT>();
+        std::cout << "," << lat1 * math::r2d<CT>();
+        std::cout << "),(" << lon2 * math::r2d<CT>();
+        std::cout << "," << lat2 * math::r2d<CT>();
+        std::cout << ")\np=(" << lon3 * math::r2d<CT>();
+        std::cout << "," << lat3 * math::r2d<CT>();
+        std::cout << ")\na1=" << a12 * math::r2d<CT>() << std::endl;
+        std::cout << "a13=" << a13 * math::r2d<CT>() << std::endl;
+        std::cout << "a312=" << a312 * math::r2d<CT>() << 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<CT>() << std::endl;
+        std::cout << "a23=" << a23 * math::r2d<CT>() << std::endl;
+        std::cout << "a321=" << a321 * math::r2d<CT>() << 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<geometry::radian>
+                > 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<CT> cross_track(earth_radius);
+        CT s34 = cross_track.apply(p3, p1, p2);
+
+        geometry::strategy::distance::haversine<CT> 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<CT> res14;
+        geometry::formula::result_inverse<CT> 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<CT>() <<
+                         "," << res14.lat2 * math::r2d<CT>() << std::endl;
+            std::cout << "delta_g4=" << delta_g4  << std::endl;
+            std::cout << "g4=" << g4 * math::r2d<CT>() << std::endl;
+            std::cout << "der=" << der  << std::endl;
+            std::cout << "M43=" << M43 << std::endl;
+            std::cout << "spherical limit=" << cos(s14/earth_radius) << std::endl;
+            std::cout << "m34=" << m34 << std::endl;
+            std::cout << "new_s14=" << s14 << std::endl;
+            std::cout << std::setprecision(16) << "dist     =" << res34.distance << std::endl;
+            std::cout << "---------end of step " << counter << std::endl<< std::endl;
+#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<double>() << ","
+                              << get<1>(p4) * math::r2d<double>() << ")"
+                  << 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<CT> 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<CT> 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<double>() << "," << res4.lat2 * math::r2d<double>() << ")"
+                  << "\np4_minus=" << p4_minus << ", p4=(" << res1.lon2 * math::r2d<double>() << "," << res1.lat2 * math::r2d<double>() << ")"
+                  << std::endl;
+
+        if (res34.distance <= p4_plus && res34.distance <= p4_minus)
+        {
+            std::cout << "Closest point computed" << std::endl;
+        }
+        else
+        {
+            std::cout << "There is a closer point nearby" << std::endl;
+        }
+#endif
+
+        return result;
+    }
+
+    Spheroid m_spheroid;
+};
+
+
+
+#ifndef DOXYGEN_NO_STRATEGY_SPECIALIZATIONS
+namespace services
+{
+
+//tags
+template <typename FormulaPolicy>
+struct tag<geographic_cross_track<FormulaPolicy> >
+{
+    typedef strategy_tag_distance_point_segment type;
+};
+
+template
+<
+        typename FormulaPolicy,
+        typename Spheroid
+>
+struct tag<geographic_cross_track<FormulaPolicy, Spheroid> >
+{
+    typedef strategy_tag_distance_point_segment type;
+};
+
+template
+<
+        typename FormulaPolicy,
+        typename Spheroid,
+        typename CalculationType
+>
+struct tag<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> >
+{
+    typedef strategy_tag_distance_point_segment type;
+};
+
+
+//return types
+template <typename FormulaPolicy, typename P, typename PS>
+struct return_type<geographic_cross_track<FormulaPolicy>, P, PS>
+    : geographic_cross_track<FormulaPolicy>::template return_type<P, PS>
+{};
+
+template
+<
+        typename FormulaPolicy,
+        typename Spheroid,
+        typename P,
+        typename PS
+>
+struct return_type<geographic_cross_track<FormulaPolicy, Spheroid>, P, PS>
+    : geographic_cross_track<FormulaPolicy, Spheroid>::template return_type<P, PS>
+{};
+
+template
+<
+        typename FormulaPolicy,
+        typename Spheroid,
+        typename CalculationType,
+        typename P,
+        typename PS
+>
+struct return_type<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType>, P, PS>
+    : geographic_cross_track<FormulaPolicy, Spheroid, CalculationType>::template return_type<P, PS>
+{};
+
+//comparable types
+template
+<
+        typename FormulaPolicy,
+        typename Spheroid,
+        typename CalculationType
+>
+struct comparable_type<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> >
+{
+    typedef geographic_cross_track
+        <
+            FormulaPolicy, Spheroid, CalculationType
+        >  type;
+};
+
+template
+<
+        typename FormulaPolicy,
+        typename Spheroid,
+        typename CalculationType
+>
+struct get_comparable<geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> >
+{
+    typedef typename comparable_type
+        <
+            geographic_cross_track<FormulaPolicy, Spheroid, CalculationType>
+        >::type comparable_type;
+public :
+    static inline comparable_type
+    apply(geographic_cross_track<FormulaPolicy, Spheroid, CalculationType> const& )
+    {
+        return comparable_type();
+    }
+};
+
+
+template
+<
+    typename FormulaPolicy,
+    typename P,
+    typename PS
+>
+struct result_from_distance<geographic_cross_track<FormulaPolicy>, P, PS>
+{
+private :
+    typedef typename geographic_cross_track
+        <
+            FormulaPolicy
+        >::template return_type<P, PS>::type return_type;
+public :
+    template <typename T>
+    static inline return_type
+    apply(geographic_cross_track<FormulaPolicy> const& , T const& distance)
+    {
+        return distance;
+    }
+};
+
+
+template <typename Point, typename PointOfSegment>
+struct default_strategy
+    <
+        point_tag, segment_tag, Point, PointOfSegment,
+        geographic_tag, geographic_tag
+    >
+{
+    typedef geographic_cross_track<> type;
+};
+
+
+template <typename PointOfSegment, typename Point>
+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

include/boost/geometry/strategies/geographic/parameters.hpp

@@ -12,7 +12,9 @@
 
 
 #include <boost/geometry/formulas/andoyer_inverse.hpp>
+#include <boost/geometry/formulas/thomas_direct.hpp>
 #include <boost/geometry/formulas/thomas_inverse.hpp>
+#include <boost/geometry/formulas/vincenty_direct.hpp>
 #include <boost/geometry/formulas/vincenty_inverse.hpp>
 
 #include <boost/mpl/assert.hpp>
@@ -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,

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 <boost/geometry/core/radian_access.hpp>
 #include <boost/geometry/core/tags.hpp>
 
-#include <boost/geometry/algorithms/detail/course.hpp>
-
 #include <boost/geometry/strategies/distance.hpp>
 #include <boost/geometry/strategies/concepts/distance_concept.hpp>
 #include <boost/geometry/strategies/spherical/distance_haversine.hpp>
@@ -373,19 +372,6 @@ public :
 
         typedef typename return_type<Point, PointOfSegment>::type return_type;
 
-#ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK
-        std::cout << "Course " << dsv(sp1) << " to " << dsv(p) << " "
-                  << crs_AD * geometry::math::r2d<return_type>() << std::endl;
-        std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " "
-                  << crs_AB * geometry::math::r2d<return_type>() << 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<return_type>() << std::endl;
-        std::cout << "Projection BD-BA " << projection2 << " : "
-                  << d_crs2 * geometry::math::r2d<return_type>() << 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<return_type>(sp1, p);
-        return_type crs_AB = geometry::detail::course<return_type>(sp1, sp2);
-        return_type crs_BA = crs_AB - geometry::math::pi<return_type>();
-        return_type crs_BD = geometry::detail::course<return_type>(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<return_type, false>
+                             (lon1, lat1, lon, lat).azimuth;
+
+        geometry::formula::result_spherical<return_type> result =
+                geometry::formula::spherical_azimuth<return_type, true>
+                    (lon1, lat1, lon2, lat2);
+        return_type crs_AB = result.azimuth;
+        return_type crs_BA = result.reverse_azimuth - geometry::math::pi<return_type>();
+
+        return_type crs_BD = geometry::formula::spherical_azimuth<return_type, false>
+                             (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<return_type>() << std::endl;
+        std::cout << "Course " << dsv(sp1) << " to " << dsv(sp2) << " "
+                  << crs_AB * geometry::math::r2d<return_type>() << std::endl;
+        std::cout << "Course " << dsv(sp2) << " to " << dsv(sp1) << " "
+                  << crs_BA * geometry::math::r2d<return_type>() << std::endl;
+        std::cout << "Course " << dsv(sp2) << " to " << dsv(p) << " "
+                  << crs_BD * geometry::math::r2d<return_type>() << std::endl;
+        std::cout << "Projection AD-AB " << projection1 << " : "
+                  << d_crs1 * geometry::math::r2d<return_type>() << std::endl;
+        std::cout << "Projection BD-BA " << projection2 << " : "
+                  << d_crs2 * geometry::math::r2d<return_type>() << std::endl;
+        std::cout << " d1: " << (d1 )
+                  << " d2: " << (d2 )
+                  << std::endl;
+#endif
+
         if (projection1 > 0.0 && projection2 > 0.0)
         {
 #ifdef BOOST_GEOMETRY_DEBUG_CROSS_TRACK

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 <boost/geometry/core/access.hpp>
 #include <boost/geometry/core/radian_access.hpp>
 
-#include <boost/geometry/algorithms/detail/course.hpp>
+#include <boost/geometry/formulas/spherical.hpp>
 
 #include <boost/geometry/util/math.hpp>
 #include <boost/geometry/util/promote_floating_point.hpp>
@@ -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<calc_t>(p1, p);
-        calc_t crs_AB = geometry::detail::course<calc_t>(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<calc_t, false>
+                             (lon1, lat1, lon, lat).azimuth;
+
+        calc_t crs_AB = geometry::formula::spherical_azimuth<calc_t, false>
+                             (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;

include/boost/geometry/strategies/strategies.hpp

@@ -85,6 +85,7 @@
 #include <boost/geometry/strategies/geographic/disjoint_segment_box.hpp>
 #include <boost/geometry/strategies/geographic/distance.hpp>
 #include <boost/geometry/strategies/geographic/distance_andoyer.hpp>
+#include <boost/geometry/strategies/geographic/distance_cross_track.hpp>
 #include <boost/geometry/strategies/geographic/distance_thomas.hpp>
 #include <boost/geometry/strategies/geographic/distance_vincenty.hpp>
 #include <boost/geometry/strategies/geographic/envelope_segment.hpp>

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 ]
     ;

test/algorithms/distance/distance_brute_force.hpp

@@ -170,7 +170,7 @@ struct distance_brute_force<Geometry1, Geometry2, Strategy, Tag1, Tag2, true>
 };
 
 
-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,

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 <iostream>
+
+#ifndef BOOST_TEST_MODULE
+#define BOOST_TEST_MODULE test_distance_geographic_pl_l
+#endif
+
+#include <boost/geometry/core/srs.hpp>
+#include <boost/test/included/unit_test.hpp>
+
+#include "test_distance_geo_common.hpp"
+
+typedef bg::cs::geographic<bg::degree> cs_type;
+typedef bg::model::point<double, 2, cs_type> point_type;
+typedef bg::model::segment<point_type> segment_type;
+typedef bg::model::multi_point<point_type> multi_point_type;
+typedef bg::model::segment<point_type> segment_type;
+typedef bg::model::linestring<point_type> linestring_type;
+typedef bg::model::multi_linestring<linestring_type> multi_linestring_type;
+
+namespace services = bg::strategy::distance::services;
+typedef bg::default_distance_result<point_type>::type return_type;
+
+typedef bg::srs::spheroid<double> stype;
+
+// Strategies for point-point distance
+
+typedef bg::strategy::distance::andoyer<stype> andoyer_pp;
+typedef bg::strategy::distance::thomas<stype> thomas_pp;
+typedef bg::strategy::distance::vincenty<stype> vincenty_pp;
+
+// Strategies for point-segment distance
+
+typedef bg::strategy::distance::geographic_cross_track<bg::strategy::andoyer, stype, double>
+        andoyer_strategy;
+
+typedef bg::strategy::distance::geographic_cross_track<bg::strategy::thomas, stype, double>
+        thomas_strategy;
+
+typedef bg::strategy::distance::geographic_cross_track<bg::strategy::vincenty, stype, double>
+        vincenty_strategy;
+
+//===========================================================================
+
+template <typename Strategy>
+inline bg::default_distance_result<point_type>::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 <typename Strategy>
+inline bg::default_distance_result<point_type>::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 <typename Strategy_pp, typename Strategy_ps>
+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<point_type, segment_type> 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 <typename Strategy_pp, typename Strategy_ps>
+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<point_type, linestring_type> 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<point_type, linestring_type> 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 <typename Strategy_pp, typename Strategy_ps>
+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 <typename Strategy_pp, typename Strategy_ps>
+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 <typename Strategy_pp, typename Strategy_ps>
+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 <typename Strategy_pp, typename Strategy_ps>
+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<multi_point_type, segment_type> 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 <typename Point, typename Strategy>
+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<Point> line_empty;
+    bg::model::multi_point<Point> multipoint_empty;
+    bg::model::multi_linestring<bg::model::linestring<Point> > multiline_empty;
+
+    Point point = from_wkt<Point>("POINT(0 0)");
+    bg::model::linestring<Point> line =
+        from_wkt<bg::model::linestring<Point> >("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<point_type>(vincenty_strategy());
+}

test/algorithms/distance/distance_se_pl_l.cpp

@@ -10,6 +10,7 @@
 // http://www.boost.org/users/license.html
 
 #include <iostream>
+#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)",

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 <iostream>
+#include <string>
+
+#include <boost/mpl/assert.hpp>
+#include <boost/type_traits/is_integral.hpp>
+#include <boost/type_traits/is_same.hpp>
+
+#include <boost/geometry/geometries/point.hpp>
+#include <boost/geometry/geometries/point_xy.hpp>
+#include <boost/geometry/geometries/segment.hpp>
+#include <boost/geometry/geometries/linestring.hpp>
+#include <boost/geometry/geometries/polygon.hpp>
+#include <boost/geometry/geometries/ring.hpp>
+#include <boost/geometry/geometries/box.hpp>
+#include <boost/geometry/geometries/multi_point.hpp>
+#include <boost/geometry/geometries/multi_linestring.hpp>
+#include <boost/geometry/geometries/multi_polygon.hpp>
+
+#include <boost/geometry/io/wkt/write.hpp>
+#include <boost/geometry/io/dsv/write.hpp>
+
+#include <boost/geometry/algorithms/num_interior_rings.hpp>
+#include <boost/geometry/algorithms/distance.hpp>
+
+#include <boost/geometry/strategies/strategies.hpp>
+
+#include <from_wkt.hpp>
+#include <string_from_type.hpp>
+
+#include "distance_brute_force.hpp"
+
+namespace bg = ::boost::geometry;
+
+//========================================================================
+
+
+template <typename T>
+struct check_equal
+{
+    template <typename Value, typename = void>
+    struct equal_to
+    {
+        static inline void apply(Value const& x, Value const& y)
+        {
+            BOOST_CHECK(x == y);
+        }
+    };
+
+    template <typename Dummy>
+    struct equal_to<double, Dummy>
+    {
+        static inline void apply(double x, double y)
+        {
+            BOOST_CHECK_CLOSE(x, y, 0.001);
+        }
+    };
+
+    template <typename Geometry1, typename Geometry2>
+    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<T>::apply(expected, detected);
+    }
+};
+
+//========================================================================
+
+template
+<
+    typename Geometry1, typename Geometry2,
+    int id1 = bg::geometry_id<Geometry1>::value,
+    int id2 = bg::geometry_id<Geometry2>::value
+>
+struct test_distance_of_geometries
+    : public test_distance_of_geometries<Geometry1, Geometry2, 0, 0>
+{};
+
+
+template <typename Geometry1, typename Geometry2>
+struct test_distance_of_geometries<Geometry1, Geometry2, 0, 0>
+{
+    template <typename DistanceType, typename Strategy>
+    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<Geometry1>(wkt1);
+        Geometry2 geometry2 = from_wkt<Geometry2>(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<typename bg::coordinate_type<Geometry1>::type>::name()
+                  << string_from_type<typename bg::coordinate_type<Geometry2>::type>::name()
+                  << " -> "
+                  << string_from_type<default_distance_result>::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 <typename Geometry1, typename Geometry2, typename Strategy>
+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

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<bg::model::polygon<P> >(poly_data_geo[3],
+                                          40007834.7139);
+
+    //force to use equator path
+    test_geometry<bg::model::polygon<P> >(poly_data_geo[4],
+                                          40075016.6856);
+
     // Multipolygon
     test_geometry<bg::model::multi_polygon<bg::model::polygon<P> > >
-                                            (multipoly_data[0], 60078934.0425);
+                                            (multipoly_data[0], 60011752.0709);
 
     // Geometries with length zero
     test_geometry<P>("POINT(0 0)", 0);

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[] = {

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<Geometry>::type).name()
         << std::endl
-        << typeid(typename bg::default_perimeter_result<Geometry>::type).name()
+        //<< typeid(typename bg::default_perimeter_result<Geometry>::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<Geometry>::type).name()
         << std::endl
-        << typeid(typename bg::default_perimeter_result<Geometry>::type).name()
+        //<< typeid(typename bg::default_perimeter_result<Geometry>::type).name()
         << std::endl
         << "perimeter : " << bg::perimeter(geometry, strategy)
         << std::endl;

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<bg::degree>(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<P1, P2>(0, 0,  180, 0, 20037.5, 0.0);
-    test_distazi<P1, P2>(0, 0, -180, 0, 20037.5, 0.0);
-    test_distazi<P1, P2>(-90, 0, 90, 0, 20037.5, 0.0);
-    test_distazi<P1, P2>(90, 0, -90, 0, 20037.5, 0.0);
+    test_distazi<P1, P2>(0, 0,  180, 0, 20003.9, 0.0);
+    test_distazi<P1, P2>(0, 0, -180, 0, 20003.9, 0.0);
+    test_distazi<P1, P2>(-90, 0, 90, 0, 20003.9, 0.0);
+    test_distazi<P1, P2>(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<P1, P2>(normlized_deg(l-44.99), -44.99, normlized_deg(l+135), 45, 20008.1, 0.0);
         test_distazi_symm<P1, P2>(normlized_deg(l-44.999), -44.999, normlized_deg(l+135), 45, 20009.4, 0.0);
         // antipodal
-        test_distazi_symm<P1, P2>(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20020.7, 0.0);
+        test_distazi_symm<P1, P2>(normlized_deg(l-45), -45, normlized_deg(l+135), 45, 20003.92, 0.0);
     }
 
     /* SQL Server gives:

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 <geometry_test_common.hpp>
+
+#include <boost/concept_check.hpp>
+
+#include <boost/geometry/core/srs.hpp>
+#include <boost/geometry/algorithms/assign.hpp>
+#include <boost/geometry/geometries/point.hpp>
+#include <boost/geometry/strategies/strategies.hpp>
+
+#include <test_common/test_point.hpp>
+
+#ifdef HAVE_TTMATH
+#  include <boost/geometry/extensions/contrib/ttmath_stub.hpp>
+#endif
+
+typedef bg::srs::spheroid<double> stype;
+
+typedef bg::strategy::andoyer andoyer_formula;
+typedef bg::strategy::thomas thomas_formula;
+typedef bg::strategy::vincenty vincenty_formula;
+
+template <typename P>
+bool non_precise_ct()
+{
+    typedef typename bg::coordinate_type<P>::type ct;
+    return boost::is_integral<ct>::value || boost::is_float<ct>::value;
+}
+
+template <typename P1, typename P2, typename FormulaPolicy>
+void test_distance(double lon1, double lat1, double lon2, double lat2)
+{
+    typedef typename bg::promote_floating_point
+        <
+            typename bg::select_calculation_type<P1, P2, void>::type
+        >::type calc_t;
+
+    calc_t tolerance = non_precise_ct<P1>() || non_precise_ct<P2>() ?
+                       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<double>(),
+                     lat1 * bg::math::d2r<double>(),
+                     lon2 * bg::math::d2r<double>(),
+                     lat2 * bg::math::d2r<double>(),
+                     stype()).distance;
+
+    bg::strategy::distance::geographic<FormulaPolicy, stype> strategy;
+    calc_t dist_strategy = strategy.apply(p1, p2);
+    BOOST_CHECK_CLOSE(dist_formula, dist_strategy, tolerance);    
+}
+
+template <typename P1, typename P2, typename FormulaPolicy>
+void test_distance_reverse(double lon1, double lat1,
+                           double lon2, double lat2)
+{
+    test_distance<P1, P2, FormulaPolicy>(lon1, lat1, lon2, lat2);
+    test_distance<P1, P2, FormulaPolicy>(lon2, lat2, lon1, lat1);
+}
+
+template <typename P1, typename P2, typename FormulaPolicy>
+void test_meridian()
+{
+    test_distance_reverse<P1, P2, FormulaPolicy>(0., 70., 0., 80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(0, 70, 0., -80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(0., -70., 0., 80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(0., -70., 0., -80.);
+
+    test_distance_reverse<P1, P2, FormulaPolicy>(0., 70., 180., 80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(0., 70., 180., -80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(0., -70., 180., 80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(0., -70., 180., -80.);
+
+    test_distance_reverse<P1, P2, FormulaPolicy>(350., 70., 170., 80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(350., 70., 170., -80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(350., -70., 170., 80.);
+    test_distance_reverse<P1, P2, FormulaPolicy>(350., -70., 170., -80.);
+}
+
+template <typename P>
+void test_all()
+{
+    test_meridian<P, P, andoyer_formula>();
+    test_meridian<P, P, thomas_formula>();
+    test_meridian<P, P, vincenty_formula>();
+}
+
+int test_main(int, char* [])
+{
+    test_all<bg::model::point<double, 2, bg::cs::geographic<bg::degree> > >();
+    test_all<bg::model::point<float, 2, bg::cs::geographic<bg::degree> > >();
+    test_all<bg::model::point<int, 2, bg::cs::geographic<bg::degree> > >();
+
+    return 0;
+}