Imported Upstream version 1.63.0upstream/1.63.0

Change-Id: Iac85556a04b7e58d63ba636dedb0986e3555714a
Signed-off-by: DongHun Kwak <dh0128.kwak@samsung.com>

15 files changed, 470 insertions, 991 deletions

diff --git a/boost/geometry/algorithms/detail/andoyer_inverse.hpp b/boost/geometry/algorithms/detail/andoyer_inverse.hpp
deleted file mode 100644
index 66ad6446e9..0000000000
--- a/boost/geometry/algorithms/detail/andoyer_inverse.hpp
+++ /dev/null
@@ -1,205 +0,0 @@
-// Boost.Geometry
-
-// Copyright (c) 2015-2016 Oracle and/or its affiliates.
-
-// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
-
-// Use, modification and distribution is subject to the Boost Software License,
-// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
-// http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_ANDOYER_INVERSE_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_ANDOYER_INVERSE_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/algorithms/detail/flattening.hpp>
-#include <boost/geometry/algorithms/detail/result_inverse.hpp>
-
-
-namespace boost { namespace geometry { namespace detail
-{
-
-/*!
-\brief The solution of the inverse problem of geodesics on latlong coordinates,
-       Forsyth-Andoyer-Lambert type approximation with first order terms.
-\author See
-    - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
-      http://www.dtic.mil/docs/citations/AD0627893
-    - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
-      http://www.dtic.mil/docs/citations/AD703541
-*/
-
-template <typename CT, bool EnableDistance, bool EnableAzimuth>
-struct andoyer_inverse
-{
-    typedef result_inverse<CT> result_type;
-
-    template <typename T1, typename T2, typename Spheroid>
-    static inline result_type apply(T1 const& lon1,
-                                    T1 const& lat1,
-                                    T2 const& lon2,
-                                    T2 const& lat2,
-                                    Spheroid const& spheroid)
-    {
-        CT const c0 = CT(0);
-        CT const c1 = CT(1);
-        CT const pi = math::pi<CT>();
-
-        result_type result;
-
-        // coordinates in radians
-
-        if ( math::equals(lon1, lon2) && math::equals(lat1, lat2) )
-        {
-            result.set(c0, c0);
-            return result;
-        }
-
-        CT const dlon = lon2 - lon1;
-        CT const sin_dlon = sin(dlon);
-        CT const cos_dlon = cos(dlon);
-        CT const sin_lat1 = sin(lat1);
-        CT const cos_lat1 = cos(lat1);
-        CT const sin_lat2 = sin(lat2);
-        CT const cos_lat2 = cos(lat2);
-
-        // H,G,T = infinity if cos_d = 1 or cos_d = -1
-        // lat1 == +-90 && lat2 == +-90
-        // lat1 == lat2 && lon1 == lon2
-        CT cos_d = sin_lat1*sin_lat2 + cos_lat1*cos_lat2*cos_dlon;
-        // on some platforms cos_d may be outside valid range
-        if (cos_d < -c1)
-            cos_d = -c1;
-        else if (cos_d > c1)
-            cos_d = c1;
-
-        CT const d = acos(cos_d); // [0, pi]
-        CT const sin_d = sin(d);  // [-1, 1]
-        
-        CT const f = detail::flattening<CT>(spheroid);
-
-        if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
-        {
-            CT const K = math::sqr(sin_lat1-sin_lat2);
-            CT const L = math::sqr(sin_lat1+sin_lat2);
-            CT const three_sin_d = CT(3) * sin_d;
-
-            CT const one_minus_cos_d = c1 - cos_d;
-            CT const one_plus_cos_d = c1 + cos_d;
-            // cos_d = 1 or cos_d = -1 means that the points are antipodal
-
-            CT const H = math::equals(one_minus_cos_d, c0) ?
-                            c0 :
-                            (d + three_sin_d) / one_minus_cos_d;
-            CT const G = math::equals(one_plus_cos_d, c0) ?
-                            c0 :
-                            (d - three_sin_d) / one_plus_cos_d;
-
-            CT const dd = -(f/CT(4))*(H*K+G*L);
-
-            CT const a = get_radius<0>(spheroid);
-
-            result.distance = a * (d + dd);
-        }
-        else
-        {
-            result.distance = c0;
-        }
-
-        if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) )
-        {
-            // sin_d = 0 <=> antipodal points
-            if (math::equals(sin_d, c0))
-            {
-                // T = inf
-                // dA = inf
-                // azimuth = -inf
-                if (lat1 <= lat2)
-                    result.azimuth = c0;
-                else
-                    result.azimuth = pi;
-            }
-            else
-            {
-                CT const c2 = CT(2);
-
-                CT A = c0;
-                CT U = c0;
-                if ( ! math::equals(cos_lat2, c0) )
-                {
-                    CT const tan_lat2 = sin_lat2/cos_lat2;
-                    CT const M = cos_lat1*tan_lat2-sin_lat1*cos_dlon;
-                    A = atan2(sin_dlon, M);
-                    CT const sin_2A = sin(c2*A);
-                    U = (f/ c2)*math::sqr(cos_lat1)*sin_2A;
-                }
-
-                CT V = c0;
-                if ( ! math::equals(cos_lat1, c0) )
-                {
-                    CT const tan_lat1 = sin_lat1/cos_lat1;
-                    CT const N = cos_lat2*tan_lat1-sin_lat2*cos_dlon;
-                    CT const B = atan2(sin_dlon, N);
-                    CT const sin_2B = sin(c2*B);
-                    V = (f/ c2)*math::sqr(cos_lat2)*sin_2B;
-                }
-
-                CT const T = d / sin_d;
-                CT const dA = V*T-U;
-
-                result.azimuth = A - dA;
-
-                // even with sin_d == 0 checked above if the second point
-                // is somewhere in the antipodal area T may still be great
-                // therefore dA may be great and the resulting azimuth
-                // may be some more or less arbitrary angle
-                if (A >= c0) // A indicates Eastern hemisphere
-                {
-                    if (dA >= c0) // A altered towards 0
-                    {
-                        if ((result.azimuth) < c0)
-                            result.azimuth = c0;
-                    }
-                    else // dA < 0, A altered towards pi
-                    {
-                        if (result.azimuth > pi)
-                            result.azimuth = pi;
-                    }
-                }
-                else // A indicates Western hemisphere
-                {
-                    if (dA <= c0) // A altered towards 0
-                    {
-                        if (result.azimuth > c0)
-                            result.azimuth = c0;
-                    }
-                    else // dA > 0, A altered towards -pi
-                    {
-                        CT const minus_pi = -pi;
-                        if ((result.azimuth) < minus_pi)
-                            result.azimuth = minus_pi;
-                    }
-                }
-            }
-        }
-        else
-        {
-            result.azimuth = c0;
-        }
-
-        return result;
-    }
-};
-
-}}} // namespace boost::geometry::detail
-
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_ANDOYER_INVERSE_HPP
diff --git a/boost/geometry/algorithms/detail/azimuth.hpp b/boost/geometry/algorithms/detail/azimuth.hpp
index 8948bf6afd..7e0d1691ed 100644
--- a/boost/geometry/algorithms/detail/azimuth.hpp
+++ b/boost/geometry/algorithms/detail/azimuth.hpp
@@ -2,8 +2,8 @@
 
 // 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, 2016.
+// Modifications copyright (c) 2014-2016, Oracle and/or its affiliates.
 
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
@@ -22,7 +22,8 @@
 #include <boost/geometry/util/math.hpp>
 
 #include <boost/geometry/algorithms/not_implemented.hpp>
-#include <boost/geometry/algorithms/detail/vincenty_inverse.hpp>
+
+#include <boost/geometry/formulas/vincenty_inverse.hpp>
 
 namespace boost { namespace geometry
 {
@@ -49,7 +50,7 @@ struct azimuth<ReturnType, geographic_tag>
     template <typename P1, typename P2, typename Spheroid>
     static inline ReturnType apply(P1 const& p1, P2 const& p2, Spheroid const& spheroid)
     {
-        return geometry::detail::vincenty_inverse<ReturnType, false, true>().apply
+        return geometry::formula::vincenty_inverse<ReturnType, false, true>().apply
                     ( get_as_radian<0>(p1), get_as_radian<1>(p1),
                       get_as_radian<0>(p2), get_as_radian<1>(p2),
                       spheroid ).azimuth;
diff --git a/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp b/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp
new file mode 100644
index 0000000000..df62a1f2f6
--- /dev/null
+++ b/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp
@@ -0,0 +1,100 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2016 Barend Gehrels, Amsterdam, the Netherlands.
+
+// 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_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP
+
+#include <set>
+
+#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
+
+namespace boost { namespace geometry
+{
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail { namespace overlay { namespace sort_by_side
+{
+
+struct ring_with_direction
+{
+    ring_identifier ring_id;
+    direction_type direction;
+    bool only_turn_on_ring;
+
+
+    inline bool operator<(ring_with_direction const& other) const
+    {
+        return this->ring_id != other.ring_id
+                ? this->ring_id < other.ring_id
+                : this->direction < other.direction;
+    }
+
+    ring_with_direction()
+        : direction(dir_unknown)
+        , only_turn_on_ring(false)
+    {}
+};
+
+struct rank_with_rings
+{
+    std::size_t rank;
+    std::set<ring_with_direction> rings;
+
+    rank_with_rings()
+        : rank(0)
+    {
+    }
+
+    inline bool all_to() const
+    {
+        for (std::set<ring_with_direction>::const_iterator it = rings.begin();
+             it != rings.end(); ++it)
+        {
+            if (it->direction == sort_by_side::dir_from)
+            {
+                return false;
+            }
+        }
+        return true;
+    }
+};
+
+template <typename Sbs>
+inline void aggregate_operations(Sbs const& sbs, std::vector<rank_with_rings>& aggregation)
+{
+    aggregation.clear();
+    for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
+    {
+        typename Sbs::rp const& ranked_point = sbs.m_ranked_points[i];
+
+        if (aggregation.empty() || aggregation.back().rank != ranked_point.rank)
+        {
+            rank_with_rings current;
+            current.rank = ranked_point.rank;
+            aggregation.push_back(current);
+        }
+
+        ring_with_direction rwd;
+        segment_identifier const& sid = ranked_point.seg_id;
+        rwd.ring_id = ring_identifier(sid.source_index, sid.multi_index, sid.ring_index);
+        rwd.direction = ranked_point.direction;
+        rwd.only_turn_on_ring = ranked_point.only_turn_on_ring;
+
+
+        aggregation.back().rings.insert(rwd);
+    }
+}
+
+
+}}} // namespace detail::overlay::sort_by_side
+#endif //DOXYGEN_NO_DETAIL
+
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP
diff --git a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
index cc55414870..2643415343 100644
--- a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
@@ -38,6 +38,7 @@ struct enrichment_info
         , count_left(0)
         , count_right(0)
         , zone(-1)
+        , only_turn_on_ring(false)
     {}
 
     // vertex to which is free travel after this IP,
@@ -57,6 +58,7 @@ struct enrichment_info
     std::size_t count_left;
     std::size_t count_right;
     signed_size_type zone; // open zone, in cluster
+    bool only_turn_on_ring; // True if it is the only turn on a ring (for clusters)
 };
 
 
diff --git a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
index 110e1be2f1..400ed3b881 100644
--- a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
+++ b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
@@ -14,6 +14,7 @@
 #include <map>
 #include <vector>
 
+#include <boost/core/ignore_unused.hpp>
 #include <boost/range.hpp>
 #include <boost/geometry/core/point_order.hpp>
 #include <boost/geometry/algorithms/detail/overlay/cluster_info.hpp>
@@ -198,7 +199,9 @@ inline signed_size_type add_turn_to_cluster(Turn const& turn,
 
     // Both operations.seg_id/fraction were already part of any cluster, and
     // these clusters are not the same. Merge of two clusters is necessary
+#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS)
     std::cout << " TODO: merge " << cid0 << " and " << cid1 << std::endl;
+#endif
     return cid0;
 }
 
@@ -314,11 +317,13 @@ inline void assign_cluster_to_turns(Turns& turns,
             typename ClusterPerSegment::const_iterator it = cluster_per_segment.find(seg_frac);
             if (it != cluster_per_segment.end())
             {
+#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS)
                 if (turn.cluster_id != -1
                         && turn.cluster_id != it->second)
                 {
                     std::cout << " CONFLICT " << std::endl;
                 }
+#endif
                 turn.cluster_id = it->second;
                 clusters[turn.cluster_id].turn_indices.insert(turn_index);
             }
@@ -392,6 +397,8 @@ inline bool is_ie_turn(segment_identifier const& ext_seg_0,
     bool const same_multi1 = ! Reverse1
                              && ext_seg_1.multi_index == other_seg_1.multi_index;
 
+    boost::ignore_unused(same_multi1);
+
     return same_multi0
             && same_multi1
             && ! is_interior<Reverse0>(ext_seg_0)
@@ -680,6 +687,7 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns,
         sbs.apply(turn_point);
 
         sbs.find_open();
+        sbs.assign_zones(for_operation);
 
         // Unset the startable flag for all 'closed' zones
         for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
@@ -706,7 +714,7 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns,
             }
         }
 
-        cinfo.open_count = sbs.open_count(turns);
+        cinfo.open_count = sbs.open_count(for_operation);
     }
 }
 
diff --git a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
index 91b0f8ae24..bbba623eee 100644
--- a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
+++ b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
@@ -37,10 +37,12 @@ struct ranked_point
         , count_left(0)
         , count_right(0)
         , operation(operation_none)
+        , only_turn_on_ring(false)
     {}
 
-    ranked_point(const Point& p, signed_size_type ti, signed_size_type oi,
-                 direction_type d, operation_type op, segment_identifier sid)
+    template <typename Op>
+    ranked_point(const Point& p, signed_size_type ti, int oi,
+                 direction_type d, Op op)
         : point(p)
         , rank(0)
         , zone(-1)
@@ -49,20 +51,22 @@ struct ranked_point
         , direction(d)
         , count_left(0)
         , count_right(0)
-        , operation(op)
-        , seg_id(sid)
+        , operation(op.operation)
+        , seg_id(op.seg_id)
+        , only_turn_on_ring(op.enriched.only_turn_on_ring)
     {}
 
     Point point;
     std::size_t rank;
     signed_size_type zone; // index of closed zone, in uu turn there would be 2 zones
     signed_size_type turn_index;
-    signed_size_type operation_index;
+    int operation_index; // 0,1
     direction_type direction;
     std::size_t count_left;
     std::size_t count_right;
     operation_type operation;
     segment_identifier seg_id;
+    bool only_turn_on_ring;
 };
 
 struct less_by_turn_index
@@ -181,11 +185,36 @@ private :
     Point m_p1, m_p2;
 };
 
+// Sorts vectors in counter clockwise order (by default)
 template <bool Reverse1, bool Reverse2, typename Point, typename Compare>
 struct side_sorter
 {
     typedef ranked_point<Point> rp;
 
+private :
+    struct include_union
+    {
+        inline bool operator()(rp const& ranked_point) const
+        {
+            // New candidate if there are no polygons on left side,
+            // but there are on right side
+            return ranked_point.count_left == 0
+                && ranked_point.count_right > 0;
+        }
+    };
+
+    struct include_intersection
+    {
+        inline bool operator()(rp const& ranked_point) const
+        {
+            // New candidate if there are two polygons on right side,
+            // and less on the left side
+            return ranked_point.count_left < 2
+                && ranked_point.count_right >= 2;
+        }
+    };
+
+public :
     inline void set_origin(Point const& origin)
     {
         m_origin = origin;
@@ -202,8 +231,8 @@ struct side_sorter
                 op.seg_id, point1, point2, point3);
         Point const& point_to = op.fraction.is_one() ? point3 : point2;
 
-        m_ranked_points.push_back(rp(point1, turn_index, op_index, dir_from, op.operation, op.seg_id));
-        m_ranked_points.push_back(rp(point_to, turn_index, op_index, dir_to, op.operation, op.seg_id));
+        m_ranked_points.push_back(rp(point1, turn_index, op_index, dir_from, op));
+        m_ranked_points.push_back(rp(point_to, turn_index, op_index, dir_to, op));
 
         if (is_origin)
         {
@@ -292,8 +321,6 @@ struct side_sorter
                     &segment_identifier::source_index
                 >(handled);
         }
-
-        assign_zones();
     }
 
     void reverse()
@@ -303,7 +330,7 @@ struct side_sorter
             return;
         }
 
-        int const last = 1 + m_ranked_points.back().rank;
+        std::size_t const last = 1 + m_ranked_points.back().rank;
 
         // Move iterator after rank==0
         bool has_first = false;
@@ -334,13 +361,18 @@ struct side_sorter
         }
     }
 
+//private :
+
+    typedef std::vector<rp> container_type;
+    container_type m_ranked_points;
+    Point m_origin;
+
+private :
+
     //! Check how many open spaces there are
-    template <typename Turns>
-    std::size_t open_count(Turns const& turns) const
+    template <typename Include>
+    inline std::size_t open_count(Include const& include_functor) const
     {
-        typedef typename boost::range_value<Turns>::type turn_type;
-        typedef typename turn_type::turn_operation_type turn_operation_type;
-
         std::size_t result = 0;
         std::size_t last_rank = 0;
         for (std::size_t i = 0; i < m_ranked_points.size(); i++)
@@ -348,31 +380,16 @@ struct side_sorter
             rp const& ranked_point = m_ranked_points[i];
 
             if (ranked_point.rank > last_rank
-                && ranked_point.direction == sort_by_side::dir_to)
+                && ranked_point.direction == sort_by_side::dir_to
+                && include_functor(ranked_point))
             {
-                // TODO: take count-left / count_right from rank itself
-                turn_type const& ranked_turn = turns[ranked_point.turn_index];
-                turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.operation_index];
-                if (ranked_op.enriched.count_left == 0
-                     && ranked_op.enriched.count_right > 0)
-                {
-                    result++;
-                    last_rank = ranked_point.rank;
-                }
+                result++;
+                last_rank = ranked_point.rank;
             }
         }
         return result;
     }
 
-//protected :
-
-    typedef std::vector<rp> container_type;
-    container_type m_ranked_points;
-    Point m_origin;
-
-private :
-
-
     std::size_t move(std::size_t index) const
     {
         std::size_t const result = index + 1;
@@ -458,7 +475,8 @@ private :
     }
 
     //! Find closed zones and assign it
-    void assign_zones()
+    template <typename Include>
+    std::size_t assign_zones(Include const& include_functor)
     {
         // Find a starting point (the first rank after an outgoing rank
         // with no polygons on the left side)
@@ -473,8 +491,7 @@ private :
                 max_rank = ranked_point.rank;
             }
             if (ranked_point.direction == sort_by_side::dir_to
-                 && ranked_point.count_left == 0
-                 && ranked_point.count_right > 0)
+                && include_functor(ranked_point))
             {
                 start_rank = ranked_point.rank + 1;
             }
@@ -510,8 +527,7 @@ private :
                 }
 
                 if (ranked_point.direction == sort_by_side::dir_to
-                     && ranked_point.count_left == 0
-                     && ranked_point.count_right > 0)
+                    && include_functor(ranked_point))
                 {
                     rank_at_next_zone = ranked_point.rank + 1;
                     if (rank_at_next_zone > max_rank)
@@ -525,8 +541,25 @@ private :
 
             ranked_point.zone = zone_id;
         }
+        return zone_id;
     }
 
+public :
+    inline std::size_t open_count(operation_type for_operation) const
+    {
+        return for_operation == operation_union
+            ? open_count(include_union())
+            : open_count(include_intersection())
+            ;
+    }
+
+    inline std::size_t assign_zones(operation_type for_operation)
+    {
+        return for_operation == operation_union
+            ? assign_zones(include_union())
+            : assign_zones(include_intersection())
+            ;
+    }
 
 };
 
diff --git a/boost/geometry/algorithms/detail/overlay/traversal.hpp b/boost/geometry/algorithms/detail/overlay/traversal.hpp
index 1156644bc3..5adc0fcf69 100644
--- a/boost/geometry/algorithms/detail/overlay/traversal.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traversal.hpp
@@ -13,6 +13,7 @@
 
 #include <boost/range.hpp>
 
+#include <boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp>
 #include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
 #include <boost/geometry/core/access.hpp>
@@ -152,8 +153,8 @@ struct traversal
         }
     }
 
-    inline bool is_visited(turn_type const& turn, turn_operation_type const& op,
-                           signed_size_type turn_index, int op_index) const
+    inline bool is_visited(turn_type const& , turn_operation_type const& op,
+                         signed_size_type , int) const
     {
         return op.visited.visited();
     }
@@ -162,39 +163,30 @@ struct traversal
                               segment_identifier const& seg_id1,
                               segment_identifier const& seg_id2) const
     {
-        if (target_operation == operation_intersection)
+        // For uu/ii, only switch sources if indicated
+        turn_type const& turn = m_turns[turn_index];
+
+        if (OverlayType == overlay_buffer)
         {
-            // For intersections always switch sources
-            return seg_id1.source_index != seg_id2.source_index;
+            // Buffer does not use source_index (always 0)
+            return turn.switch_source
+                    ? seg_id1.multi_index != seg_id2.multi_index
+                    : seg_id1.multi_index == seg_id2.multi_index;
         }
-        else if (target_operation == operation_union)
-        {
-            // For uu, only switch sources if indicated
-            turn_type const& turn = m_turns[turn_index];
-
-            if (OverlayType == overlay_buffer)
-            {
-                // Buffer does not use source_index (always 0)
-                return turn.switch_source
-                        ? seg_id1.multi_index != seg_id2.multi_index
-                        : seg_id1.multi_index == seg_id2.multi_index;
-            }
 
 #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
-            if (turn.switch_source == 1)
-            {
-                std::cout << "Switch source at " << turn_index << std::endl;
-            }
-            else
-            {
-                std::cout << "DON'T SWITCH SOURCES at " << turn_index << std::endl;
-            }
-#endif
-            return turn.switch_source
-                    ? seg_id1.source_index != seg_id2.source_index
-                    : seg_id1.source_index == seg_id2.source_index;
+        if (turn.switch_source)
+        {
+            std::cout << "Switch source at " << turn_index << std::endl;
         }
-        return false;
+        else
+        {
+            std::cout << "DON'T SWITCH SOURCES at " << turn_index << std::endl;
+        }
+#endif
+        return turn.switch_source
+                ? seg_id1.source_index != seg_id2.source_index
+                : seg_id1.source_index == seg_id2.source_index;
     }
 
     inline
@@ -273,10 +265,6 @@ struct traversal
                 segment_identifier const& seg_id,
                 int& selected_op_index) const
     {
-        // For "ii", take the other one (alternate)
-        //           UNLESS the other one is already visited
-        // For "uu", take the same one (see above);
-
         bool result = false;
 
         for (int i = 0; i < 2; i++)
@@ -347,13 +335,10 @@ struct traversal
         return true;
     }
 
-    inline bool select_from_cluster(signed_size_type& turn_index,
+    inline bool select_from_cluster_union(signed_size_type& turn_index,
         int& op_index, signed_size_type start_turn_index,
         sbs_type const& sbs, bool is_touching) const
     {
-        bool const is_union = target_operation == operation_union;
-        bool const is_intersection = target_operation == operation_intersection;
-
         std::size_t selected_rank = 0;
         std::size_t min_rank = 0;
         bool result = false;
@@ -386,11 +371,8 @@ struct traversal
                 && (ranked_point.rank > min_rank
                     || ranked_turn.both(operation_continue)))
             {
-                if ((is_union
-                     && ranked_op.enriched.count_left == 0
+                if (ranked_op.enriched.count_left == 0
                      && ranked_op.enriched.count_right > 0)
-                || (is_intersection
-                     && ranked_op.enriched.count_right == 2))
                 {
                     if (result && ranked_point.turn_index != start_turn_index)
                     {
@@ -401,16 +383,6 @@ struct traversal
                     turn_index = ranked_point.turn_index;
                     op_index = ranked_point.operation_index;
 
-                    if (is_intersection
-                        && ranked_turn.both(operation_intersection)
-                        && ranked_op.visited.finalized())
-                    {
-                        // Override:
-                        // For a ii turn, even though one operation might be selected,
-                        // it should take the other one if the first one is used in a completed ring
-                        op_index = 1 - ranked_point.operation_index;
-                    }
-
                     result = true;
                     selected_rank = ranked_point.rank;
                 }
@@ -423,10 +395,94 @@ struct traversal
         return result;
     }
 
+    inline bool analyze_cluster_intersection(signed_size_type& turn_index,
+                int& op_index,
+                sbs_type const& sbs) const
+    {
+        std::vector<sort_by_side::rank_with_rings> aggregation;
+        sort_by_side::aggregate_operations(sbs, aggregation);
+
+        std::size_t selected_rank = 0;
+
+        for (std::size_t i = 0; i < aggregation.size(); i++)
+        {
+            sort_by_side::rank_with_rings const& rwr = aggregation[i];
+
+            if (i > 1
+                && i - 1 == selected_rank
+                && rwr.rings.size() == 1)
+            {
+                sort_by_side::ring_with_direction const& rwd = *rwr.rings.begin();
+
+                if (rwd.only_turn_on_ring)
+                {
+                    // Find if this arriving ring was leaving previously
+                    sort_by_side::ring_with_direction leaving = rwd;
+                    leaving.direction = sort_by_side::dir_to;
+
+                    sort_by_side::rank_with_rings const& previous = aggregation[i - 1];
+
+                    if (previous.rings.size() == 1
+                        && previous.rings.count(leaving) == 1)
+                    {
+                        // It arrives back - if this is one of the selected, unselect it
+                        selected_rank = 0;
+                    }
+                }
+            }
+
+            if (rwr.all_to())
+            {
+                if (selected_rank == 0)
+                {
+                    // Take the first (= right) where segments leave,
+                    // having the polygon on the right side
+                    selected_rank = rwr.rank;
+                }
+            }
+        }
+
+        if (selected_rank > 0)
+        {
+            std::size_t selected_index = sbs.m_ranked_points.size();
+            for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
+            {
+                typename sbs_type::rp const& ranked_point = sbs.m_ranked_points[i];
+
+                if (ranked_point.rank == selected_rank)
+                {
+                    turn_type const& ranked_turn = m_turns[ranked_point.turn_index];
+                    turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.operation_index];
+
+                    if (ranked_op.visited.finalized())
+                    {
+                        // This direction is already traveled before, the same
+                        // cannot be traveled again
+                        return false;
+                    }
+
+                    // Take the last turn from this rank
+                    selected_index = i;
+                }
+            }
+
+            if (selected_index < sbs.m_ranked_points.size())
+            {
+                typename sbs_type::rp const& ranked_point = sbs.m_ranked_points[selected_index];
+                turn_index = ranked_point.turn_index;
+                op_index = ranked_point.operation_index;
+                return true;
+            }
+        }
+
+        return false;
+    }
+
     inline bool select_turn_from_cluster(signed_size_type& turn_index,
             int& op_index, bool& is_touching,
             signed_size_type start_turn_index,
-            segment_identifier const& previous_seg_id) const
+            segment_identifier const& previous_seg_id,
+            bool is_start) const
     {
         bool const is_union = target_operation == operation_union;
 
@@ -488,30 +544,76 @@ struct traversal
 
         sbs.apply(turn.point);
 
-#if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
-        is_touching = is_union && cinfo.open_count > 1;
-        if (is_touching)
+        bool result = false;
+
+        if (is_union)
         {
-            if (cinfo.switch_source)
+    #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
+            is_touching = cinfo.open_count > 1;
+            if (is_touching)
             {
-                is_touching = false;
-                std::cout << "CLUSTER: SWITCH SOURCES at " << turn_index << std::endl;
+                if (cinfo.switch_source)
+                {
+                    is_touching = false;
+                    std::cout << "CLUSTER: SWITCH SOURCES at " << turn_index << std::endl;
+                }
+                else
+                {
+                    std::cout << "CLUSTER: CONTINUE at " << turn_index << std::endl;
+                }
             }
-            else
+    #else
+            is_touching = cinfo.open_count > 1 && ! cinfo.switch_source;
+    #endif
+
+            if (is_touching)
             {
-                std::cout << "CLUSTER: CONTINUE at " << turn_index << std::endl;
+                sbs.reverse();
             }
+
+            result = select_from_cluster_union(turn_index, op_index, start_turn_index, sbs,
+                                       is_touching);
         }
-#else
-        is_touching = is_union && cinfo.open_count > 1 && ! cinfo.switch_source;
-#endif
-        if (is_touching)
+        else
         {
-            sbs.reverse();
+            if (is_start
+                && turn.both(operation_intersection)
+                && turn.operations[op_index].enriched.only_turn_on_ring)
+            {
+                // For an ii (usually interior ring), only turn on ring,
+                // reverse to take first exit
+                sbs.reverse();
+            }
+
+            result = analyze_cluster_intersection(turn_index, op_index, sbs);
         }
+        return result;
+    }
 
-        return select_from_cluster(turn_index, op_index, start_turn_index, sbs,
-                                   is_touching);
+    inline bool analyze_ii_intersection(signed_size_type& turn_index, int& op_index,
+                    turn_type const& current_turn,
+                    segment_identifier const& previous_seg_id)
+    {
+        sbs_type sbs;
+
+        // Add this turn to the sort-by-side sorter
+        bool has_origin = false;
+        for (int i = 0; i < 2; i++)
+        {
+            turn_operation_type const& op = current_turn.operations[i];
+            bool const is_origin = op.seg_id.source_index
+                                   == previous_seg_id.source_index;
+            has_origin = has_origin || is_origin;
+            sbs.add(op, turn_index, i, m_geometry1, m_geometry2, is_origin);
+        }
+
+        if (! has_origin)
+        {
+            return false;
+        }
+
+        sbs.apply(current_turn.point);
+        return analyze_cluster_intersection(turn_index, op_index, sbs);
     }
 
     inline void change_index_for_self_turn(signed_size_type& to_vertex_index,
@@ -607,7 +709,7 @@ struct traversal
         return true;
     }
 
-    bool select_turn(signed_size_type start_turn_index,
+    bool select_turn(signed_size_type start_turn_index, int start_op_index,
                      signed_size_type& turn_index,
                      int& op_index,
                      bool& is_touching,
@@ -616,10 +718,37 @@ struct traversal
                      segment_identifier const& previous_seg_id,
                      bool is_start)
     {
-        if (m_turns[turn_index].cluster_id >= 0)
+        turn_type const& current_turn = m_turns[turn_index];
+
+        if (target_operation == operation_intersection)
+        {
+            bool const back_at_start_cluster
+                    = current_turn.cluster_id >= 0
+                    && m_turns[start_turn_index].cluster_id == current_turn.cluster_id;
+
+            if (turn_index == start_turn_index || back_at_start_cluster)
+            {
+                // Intersection can always be finished if returning
+                turn_index = start_turn_index;
+                op_index = start_op_index;
+                return true;
+            }
+
+            if (current_turn.cluster_id < 0
+                && current_turn.both(operation_intersection))
+            {
+                if (analyze_ii_intersection(turn_index, op_index, current_turn,
+                            previous_seg_id))
+                {
+                    return true;
+                }
+            }
+        }
+
+        if (current_turn.cluster_id >= 0)
         {
             if (! select_turn_from_cluster(turn_index, op_index, is_touching,
-                    start_turn_index, previous_seg_id))
+                    start_turn_index, previous_seg_id, is_start))
             {
                 return false;
             }
@@ -631,8 +760,6 @@ struct traversal
         }
         else
         {
-            turn_type const& current_turn = m_turns[turn_index];
-
             op_index = starting_operation_index(current_turn);
             if (op_index == -1)
             {
diff --git a/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp b/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp
index 104bd6b8e7..e0dfee19a8 100644
--- a/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp
@@ -64,9 +64,7 @@ struct traversal_ring_creator
         , m_clusters(clusters)
         , m_robust_policy(robust_policy)
         , m_visitor(visitor)
-        , m_has_uu(false)
     {
-
     }
 
     template <typename Ring>
@@ -123,7 +121,8 @@ struct traversal_ring_creator
         }
 
         bool is_touching = false;
-        if (! m_trav.select_turn(start_turn_index, turn_index, op_index,
+        if (! m_trav.select_turn(start_turn_index, start_op_index,
+                turn_index, op_index,
                 is_touching,
                 previous_op_index, previous_turn_index, previous_seg_id,
                 is_start))
@@ -189,6 +188,18 @@ struct traversal_ring_creator
             return traverse_error_none;
         }
 
+        if (start_turn.cluster_id >= 0)
+        {
+            turn_type const& turn = m_turns[current_turn_index];
+            if (turn.cluster_id == start_turn.cluster_id)
+            {
+                turn_operation_type& op = m_turns[start_turn_index].operations[current_op_index];
+                op.visited.set_finished();
+                m_visitor.visit_traverse(m_turns, m_turns[current_turn_index], start_op, "Early finish (cluster)");
+                return traverse_error_none;
+            }
+        }
+
         std::size_t const max_iterations = 2 + 2 * m_turns.size();
         for (std::size_t i = 0; i <= max_iterations; i++)
         {
@@ -276,49 +287,21 @@ struct traversal_ring_creator
 
     template <typename Rings>
     void iterate(Rings& rings, std::size_t& finalized_ring_size,
-                 typename Backtrack::state_type& state,
-                 int pass)
+                 typename Backtrack::state_type& state)
     {
-        if (pass == 1)
-        {
-            if (target_operation == operation_intersection)
-            {
-                // Second pass currently only used for uu
-                return;
-            }
-            if (! m_has_uu)
-            {
-                // There is no uu found in first pass
-                return;
-            }
-        }
-
-        // Iterate through all unvisited points
         for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index)
         {
-            turn_type const& start_turn = m_turns[turn_index];
+            turn_type const& turn = m_turns[turn_index];
 
-            if (start_turn.discarded || start_turn.blocked())
+            if (turn.discarded || turn.blocked())
             {
                 // Skip discarded and blocked turns
                 continue;
             }
-            if (target_operation == operation_union)
-            {
-                if (start_turn.both(operation_union))
-                {
-                    // Start with a uu-turn only in the second pass
-                    m_has_uu = true;
-                    if (pass == 0)
-                    {
-                        continue;
-                    }
-                }
-            }
 
             for (int op_index = 0; op_index < 2; op_index++)
             {
-                traverse_with_operation(start_turn, turn_index, op_index,
+                traverse_with_operation(turn, turn_index, op_index,
                         rings, finalized_ring_size, state);
             }
         }
@@ -333,10 +316,6 @@ private:
     Clusters const& m_clusters;
     RobustPolicy const& m_robust_policy;
     Visitor& m_visitor;
-
-    // Next member is only used for operation union
-    bool m_has_uu;
-
 };
 
 }} // namespace detail::overlay
diff --git a/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp b/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp
index 9381c66e0e..183131c74b 100644
--- a/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp
@@ -71,8 +71,8 @@ struct traversal_switch_detector
     {
         if (turn.cluster_id == -1)
         {
-            // If it is a uu-turn (non clustered), it is never same zone
-            return ! turn.both(operation_union);
+            // If it is a uu/ii-turn (non clustered), it is never same zone
+            return ! (turn.both(operation_union) || turn.both(operation_intersection));
         }
 
         // It is a cluster, check zones of both operations
@@ -110,11 +110,11 @@ struct traversal_switch_detector
         for (set_iterator sit = ring_turn_indices.begin();
              sit != ring_turn_indices.end(); ++sit)
         {
-            int const turn_index = *sit;
+            signed_size_type const turn_index = *sit;
             turn_type const& turn = m_turns[turn_index];
             if (! connects_same_zone(turn))
             {
-                // This is a non clustered uu-turn, or a cluster connecting different 'zones'
+                // This is a non clustered uu/ii-turn, or a cluster connecting different 'zones'
                 continue;
             }
 
@@ -132,6 +132,56 @@ struct traversal_switch_detector
         }
     }
 
+    void check_turns_per_ring(ring_identifier const& ring_id,
+            std::set<signed_size_type> const& ring_turn_indices)
+    {
+        bool only_turn_on_ring = true;
+        if (ring_turn_indices.size() > 1)
+        {
+            // More turns on this ring. Only leave only_turn_on_ring true
+            // if they are all of the same cluster
+            int cluster_id = -1;
+            for (set_iterator sit = ring_turn_indices.begin();
+                 sit != ring_turn_indices.end(); ++sit)
+            {
+                turn_type const& turn = m_turns[*sit];
+                if (turn.cluster_id == -1)
+                {
+                    // Unclustered turn - and there are 2 or more turns
+                    // so the ring has different turns
+                    only_turn_on_ring = false;
+                    break;
+                }
+
+                // Clustered turn, check if it is the first or same as previous
+                if (cluster_id == -1)
+                {
+                    cluster_id = turn.cluster_id;
+                }
+                else if (turn.cluster_id != cluster_id)
+                {
+                    only_turn_on_ring = false;
+                    break;
+                }
+            }
+        }
+
+        // Assign result to matching operation (a turn is always on two rings)
+        for (set_iterator sit = ring_turn_indices.begin();
+             sit != ring_turn_indices.end(); ++sit)
+        {
+            turn_type& turn = m_turns[*sit];
+            for (int i = 0; i < 2; i++)
+            {
+                turn_operation_type& op = turn.operations[i];
+                if (ring_id_by_seg_id(op.seg_id) == ring_id)
+                {
+                    op.enriched.only_turn_on_ring = only_turn_on_ring;
+                }
+            }
+        }
+    }
+
     void propagate_region(ring_identifier const& ring_id, int region_id)
     {
         std::map<ring_identifier, std::set<signed_size_type> >::const_iterator it = m_turns_per_ring.find(ring_id);
@@ -168,6 +218,7 @@ struct traversal_switch_detector
              = m_turns_per_ring.begin(); it != m_turns_per_ring.end(); ++it)
         {
             create_region(it->first, it->second);
+            check_turns_per_ring(it->first, it->second);
         }
 
         // Now that all regions are filled, assign switch_source property
@@ -204,7 +255,7 @@ struct traversal_switch_detector
             cinfo.switch_source = regions.size() == 1;
         }
 
-        // Iterate through all uu turns (non-clustered)
+        // Iterate through all uu/ii turns (non-clustered)
         for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index)
         {
             turn_type& turn = m_turns[turn_index];
@@ -212,9 +263,9 @@ struct traversal_switch_detector
             if (turn.discarded
                     || turn.blocked()
                     || turn.cluster_id >= 0
-                    || ! turn.both(operation_union))
+                    || ! (turn.both(operation_union) || turn.both(operation_intersection)))
             {
-                // Skip discarded, blocked, non-uu and clustered turns
+                // Skip discarded, blocked, non-uu/ii and clustered turns
                 continue;
             }
 
@@ -242,9 +293,10 @@ struct traversal_switch_detector
         {
             turn_type const& turn = m_turns[turn_index];
 
-            if (turn.both(operation_union) && turn.cluster_id < 0)
+            if ((turn.both(operation_union) || turn.both(operation_intersection))
+                 && turn.cluster_id < 0)
             {
-                std::cout << "UU SWITCH RESULT "
+                std::cout << "UU/II SWITCH RESULT "
                              << turn_index << " -> "
                           << turn.switch_source << std::endl;
             }
diff --git a/boost/geometry/algorithms/detail/overlay/traverse.hpp b/boost/geometry/algorithms/detail/overlay/traverse.hpp
index 2d2933ebdd..f01e50eb03 100644
--- a/boost/geometry/algorithms/detail/overlay/traverse.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traverse.hpp
@@ -97,10 +97,7 @@ public :
 
         typename Backtrack::state_type state;
 
-        for (int pass = 0; pass < 2; pass++)
-        {
-            trav.iterate(rings, finalized_ring_size, state, pass);
-        }
+        trav.iterate(rings, finalized_ring_size, state);
     }
 };
 
diff --git a/boost/geometry/algorithms/detail/overlay/turn_info.hpp b/boost/geometry/algorithms/detail/overlay/turn_info.hpp
index 73f266a04a..e09af126c6 100644
--- a/boost/geometry/algorithms/detail/overlay/turn_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/turn_info.hpp
@@ -13,6 +13,7 @@
 #include <boost/array.hpp>
 
 #include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/algorithms/detail/signed_size_type.hpp>
 #include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp>
 #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp>
 
@@ -88,7 +89,7 @@ struct turn_info
 
     Point point;
     method_type method;
-    int cluster_id;
+    signed_size_type cluster_id; // For multiple turns on same location, >= 0. Else -1
     bool discarded;
     bool colocated;
     bool switch_source; // For u/u turns which can either switch or not
diff --git a/boost/geometry/algorithms/detail/result_inverse.hpp b/boost/geometry/algorithms/detail/result_inverse.hpp
deleted file mode 100644
index 01a1997e49..0000000000
--- a/boost/geometry/algorithms/detail/result_inverse.hpp
+++ /dev/null
@@ -1,44 +0,0 @@
-// Boost.Geometry
-
-// Copyright (c) 2015 Oracle and/or its affiliates.
-
-// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
-
-// Use, modification and distribution is subject to the Boost Software License,
-// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
-// http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_RESULT_INVERSE_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_RESULT_INVERSE_HPP
-
-
-#include <boost/math/constants/constants.hpp>
-
-#include <boost/geometry/core/radius.hpp>
-#include <boost/geometry/core/srs.hpp>
-
-#include <boost/geometry/util/math.hpp>
-
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
-
-namespace boost { namespace geometry { namespace detail
-{
-
-template <typename T>
-struct result_inverse
-{
-    void set(T const& d, T const& a)
-    {
-        distance = d;
-        azimuth = a;
-    }
-
-    T distance;
-    T azimuth;
-};
-
-}}} // namespace boost::geometry::detail
-
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_RESULT_INVERSE_HPP
diff --git a/boost/geometry/algorithms/detail/thomas_inverse.hpp b/boost/geometry/algorithms/detail/thomas_inverse.hpp
deleted file mode 100644
index 96b237e054..0000000000
--- a/boost/geometry/algorithms/detail/thomas_inverse.hpp
+++ /dev/null
@@ -1,191 +0,0 @@
-// Boost.Geometry
-
-// Copyright (c) 2015 Oracle and/or its affiliates.
-
-// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
-
-// Use, modification and distribution is subject to the Boost Software License,
-// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
-// http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_THOMAS_INVERSE_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_THOMAS_INVERSE_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/algorithms/detail/flattening.hpp>
-#include <boost/geometry/algorithms/detail/result_inverse.hpp>
-
-namespace boost { namespace geometry { namespace detail
-{
-
-/*!
-\brief The solution of the inverse problem of geodesics on latlong coordinates,
-       Forsyth-Andoyer-Lambert type approximation with second order terms.
-\author See
-    - Technical Report: PAUL D. THOMAS, MATHEMATICAL MODELS FOR NAVIGATION SYSTEMS, 1965
-      http://www.dtic.mil/docs/citations/AD0627893
-    - Technical Report: PAUL D. THOMAS, SPHEROIDAL GEODESICS, REFERENCE SYSTEMS, AND LOCAL GEOMETRY, 1970
-      http://www.dtic.mil/docs/citations/AD703541
-*/
-template <typename CT, bool EnableDistance, bool EnableAzimuth>
-struct thomas_inverse
-{
-    typedef result_inverse<CT> result_type;
-
-    template <typename T1, typename T2, typename Spheroid>
-    static inline result_type apply(T1 const& lon1,
-                                    T1 const& lat1,
-                                    T2 const& lon2,
-                                    T2 const& lat2,
-                                    Spheroid const& spheroid)
-    {
-        result_type result;
-
-        // coordinates in radians
-
-        if ( math::equals(lon1, lon2)
-          && math::equals(lat1, lat2) )
-        {
-            result.set(CT(0), CT(0));
-            return result;
-        }
-
-        CT const f = detail::flattening<CT>(spheroid);
-        CT const one_minus_f = CT(1) - f;
-
-//        CT const tan_theta1 = one_minus_f * tan(lat1);
-//        CT const tan_theta2 = one_minus_f * tan(lat2);
-//        CT const theta1 = atan(tan_theta1);
-//        CT const theta2 = atan(tan_theta2);
-
-        CT const pi_half = math::pi<CT>() / CT(2);
-        CT const theta1 = math::equals(lat1, pi_half) ? lat1 :
-                          math::equals(lat1, -pi_half) ? lat1 :
-                          atan(one_minus_f * tan(lat1));
-        CT const theta2 = math::equals(lat2, pi_half) ? lat2 :
-                          math::equals(lat2, -pi_half) ? lat2 :
-                          atan(one_minus_f * tan(lat2));
-
-        CT const theta_m = (theta1 + theta2) / CT(2);
-        CT const d_theta_m = (theta2 - theta1) / CT(2);
-        CT const d_lambda = lon2 - lon1;
-        CT const d_lambda_m = d_lambda / CT(2);
-
-        CT const sin_theta_m = sin(theta_m);
-        CT const cos_theta_m = cos(theta_m);
-        CT const sin_d_theta_m = sin(d_theta_m);
-        CT const cos_d_theta_m = cos(d_theta_m);
-        CT const sin2_theta_m = math::sqr(sin_theta_m);
-        CT const cos2_theta_m = math::sqr(cos_theta_m);
-        CT const sin2_d_theta_m = math::sqr(sin_d_theta_m);
-        CT const cos2_d_theta_m = math::sqr(cos_d_theta_m);
-        CT const sin_d_lambda_m = sin(d_lambda_m);
-        CT const sin2_d_lambda_m = math::sqr(sin_d_lambda_m);
-
-        CT const H = cos2_theta_m - sin2_d_theta_m;
-        CT const L = sin2_d_theta_m + H * sin2_d_lambda_m;
-        CT const cos_d = CT(1) - CT(2) * L;
-        CT const d = acos(cos_d);
-        CT const sin_d = sin(d);
-
-        CT const one_minus_L = CT(1) - L;
-
-        if ( math::equals(sin_d, CT(0))
-          || math::equals(L, CT(0))
-          || math::equals(one_minus_L, CT(0)) )
-        {
-            result.set(CT(0), CT(0));
-            return result;
-        }
-
-        CT const U = CT(2) * sin2_theta_m * cos2_d_theta_m / one_minus_L;
-        CT const V = CT(2) * sin2_d_theta_m * cos2_theta_m / L;
-        CT const X = U + V;
-        CT const Y = U - V;
-        CT const T = d / sin_d;
-        //CT const D = CT(4) * math::sqr(T);
-        //CT const E = CT(2) * cos_d;
-        //CT const A = D * E;
-        //CT const B = CT(2) * D;
-        //CT const C = T - (A - E) / CT(2);
-    
-        if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
-        {
-            //CT const n1 = X * (A + C*X);
-            //CT const n2 = Y * (B + E*Y);
-            //CT const n3 = D*X*Y;
-
-            //CT const f_sqr = math::sqr(f);
-            //CT const f_sqr_per_64 = f_sqr / CT(64);
-
-            CT const delta1d = f * (T*X-Y) / CT(4);
-            //CT const delta2d = f_sqr_per_64 * (n1 - n2 + n3);
-
-            CT const a = get_radius<0>(spheroid);
-
-            result.distance = a * sin_d * (T - delta1d);
-            //double S2 = a * sin_d * (T - delta1d + delta2d);
-        }
-        else
-        {
-            result.distance = CT(0);
-        }
-    
-        if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) )
-        {
-            // NOTE: if both cos_latX == 0 then below we'd have 0 * INF
-            // it's a situation when the endpoints are on the poles +-90 deg
-            // in this case the azimuth could either be 0 or +-pi
-            // but above always 0 is returned
-
-            // may also be used to calculate distance21
-            //CT const D = CT(4) * math::sqr(T);
-            CT const E = CT(2) * cos_d;
-            //CT const A = D * E;
-            //CT const B = CT(2) * D;
-            // may also be used to calculate distance21
-            CT const f_sqr = math::sqr(f);
-            CT const f_sqr_per_64 = f_sqr / CT(64);
-
-            CT const F = CT(2)*Y-E*(CT(4)-X);
-            //CT const M = CT(32)*T-(CT(20)*T-A)*X-(B+CT(4))*Y;
-            CT const G = f*T/CT(2) + f_sqr_per_64;
-            CT const tan_d_lambda = tan(d_lambda);
-            CT const Q = -(F*G*tan_d_lambda) / CT(4);
-
-            CT const d_lambda_p = (d_lambda + Q) / CT(2);
-            CT const tan_d_lambda_p = tan(d_lambda_p);
-
-            CT const v = atan2(cos_d_theta_m, sin_theta_m * tan_d_lambda_p);
-            CT const u = atan2(-sin_d_theta_m, cos_theta_m * tan_d_lambda_p);
-
-            CT const pi = math::pi<CT>();
-            CT alpha1 = v + u;
-            if ( alpha1 > pi )
-            {
-                alpha1 -= CT(2) * pi;
-            }
-
-            result.azimuth = alpha1;
-        }
-        else
-        {
-            result.azimuth = CT(0);
-        }
-
-        return result;
-    }
-};
-
-}}} // namespace boost::geometry::detail
-
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_THOMAS_INVERSE_HPP
diff --git a/boost/geometry/algorithms/detail/vincenty_direct.hpp b/boost/geometry/algorithms/detail/vincenty_direct.hpp
deleted file mode 100644
index 1c47a0f68d..0000000000
--- a/boost/geometry/algorithms/detail/vincenty_direct.hpp
+++ /dev/null
@@ -1,177 +0,0 @@
-// Boost.Geometry
-
-// 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.
-
-// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
-
-// Use, modification and distribution is subject to the Boost Software License,
-// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
-// http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_DIRECT_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_DIRECT_HPP
-
-
-#include <boost/math/constants/constants.hpp>
-
-#include <boost/geometry/core/radius.hpp>
-#include <boost/geometry/core/srs.hpp>
-
-#include <boost/geometry/util/math.hpp>
-
-#include <boost/geometry/algorithms/detail/flattening.hpp>
-
-
-#ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS
-#define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000
-#endif
-
-
-namespace boost { namespace geometry { namespace detail
-{
-
-template <typename T>
-struct result_direct
-{
-    void set(T const& lo2, T const& la2)
-    {
-        lon2 = lo2;
-        lat2 = la2;
-    }
-    T lon2;
-    T lat2;
-};
-
-/*!
-\brief The solution of the direct problem of geodesics on latlong coordinates, after Vincenty, 1975
-\author See
-    - http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
-    - http://www.icsm.gov.au/gda/gdav2.3.pdf
-\author Adapted from various implementations to get it close to the original document
-    - http://www.movable-type.co.uk/scripts/LatLongVincenty.html
-    - http://exogen.case.edu/projects/geopy/source/geopy.distance.html
-    - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
-
-*/
-template <typename CT>
-struct vincenty_direct
-{
-    typedef result_direct<CT> result_type;
-
-public:
-    template <typename T, typename Dist, typename Azi, typename Spheroid>
-    static inline result_type apply(T const& lo1,
-                                    T const& la1,
-                                    Dist const& distance,
-                                    Azi const& azimuth12,
-                                    Spheroid const& spheroid)
-    {
-        result_type result;
-
-        CT const lon1 = lo1;
-        CT const lat1 = la1;
-
-        if ( math::equals(distance, Dist(0)) || distance < Dist(0) )
-        {
-            result.set(lon1, lat1);
-            return result;
-        }
-
-        CT const radius_a = CT(get_radius<0>(spheroid));
-        CT const radius_b = CT(get_radius<2>(spheroid));
-        CT const flattening = geometry::detail::flattening<CT>(spheroid);
-
-        CT const sin_azimuth12 = sin(azimuth12);
-        CT const cos_azimuth12 = cos(azimuth12);
-
-        // U: reduced latitude, defined by tan U = (1-f) tan phi
-        CT const one_min_f = CT(1) - flattening;
-        CT const tan_U1 = one_min_f * tan(lat1);
-        CT const sigma1 = atan2(tan_U1, cos_azimuth12); // (1)
-
-        // may be calculated from tan using 1 sqrt()
-        CT const U1 = atan(tan_U1);
-        CT const sin_U1 = sin(U1);
-        CT const cos_U1 = cos(U1);
-
-        CT const sin_alpha = cos_U1 * sin_azimuth12; // (2)
-        CT const sin_alpha_sqr = math::sqr(sin_alpha);
-        CT const cos_alpha_sqr = CT(1) - sin_alpha_sqr;
-
-        CT const b_sqr = radius_b * radius_b;
-        CT const u_sqr = cos_alpha_sqr * (radius_a * radius_a - b_sqr) / b_sqr;
-        CT const A = CT(1) + (u_sqr/CT(16384)) * (CT(4096) + u_sqr*(CT(-768) + u_sqr*(CT(320) - u_sqr*CT(175)))); // (3)
-        CT const B = (u_sqr/CT(1024))*(CT(256) + u_sqr*(CT(-128) + u_sqr*(CT(74) - u_sqr*CT(47)))); // (4)
-
-        CT s_div_bA = distance / (radius_b * A);
-        CT sigma = s_div_bA; // (7)
-
-        CT previous_sigma;
-        CT sin_sigma;
-        CT cos_sigma;
-        CT cos_2sigma_m;
-        CT cos_2sigma_m_sqr;
-
-        int counter = 0; // robustness
-
-        do
-        {
-            previous_sigma = sigma;
-
-            CT const two_sigma_m = CT(2) * sigma1 + sigma; // (5)
-
-            sin_sigma = sin(sigma);
-            cos_sigma = cos(sigma);
-            CT const sin_sigma_sqr = math::sqr(sin_sigma);
-            cos_2sigma_m = cos(two_sigma_m);
-            cos_2sigma_m_sqr = math::sqr(cos_2sigma_m);
-
-            CT const delta_sigma = B * sin_sigma * (cos_2sigma_m
-                                        + (B/CT(4)) * ( cos_sigma * (CT(-1) + CT(2)*cos_2sigma_m_sqr)
-                                            - (B/CT(6) * cos_2sigma_m * (CT(-3)+CT(4)*sin_sigma_sqr) * (CT(-3)+CT(4)*cos_2sigma_m_sqr)) )); // (6)
-
-            sigma = s_div_bA + delta_sigma; // (7)
-
-            ++counter; // robustness
-
-        } while ( geometry::math::abs(previous_sigma - sigma) > CT(1e-12)
-               //&& geometry::math::abs(sigma) < pi
-               && counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness
-
-        {
-            result.lat2
-                = atan2( sin_U1 * cos_sigma + cos_U1 * sin_sigma * cos_azimuth12,
-                         one_min_f * math::sqrt(sin_alpha_sqr + math::sqr(sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azimuth12))); // (8)
-        }
-
-        {
-            CT const lambda = atan2( sin_sigma * sin_azimuth12,
-                                     cos_U1 * cos_sigma - sin_U1 * sin_sigma * cos_azimuth12); // (9)
-            CT const C = (flattening/CT(16)) * cos_alpha_sqr * ( CT(4) + flattening * ( CT(4) - CT(3) * cos_alpha_sqr ) ); // (10)
-            CT const L = lambda - (CT(1) - C) * flattening * sin_alpha
-                            * ( sigma + C * sin_sigma * ( cos_2sigma_m + C * cos_sigma * ( CT(-1) + CT(2) * cos_2sigma_m_sqr ) ) ); // (11)
-
-            result.lon2 = lon1 + L;
-        }
-
-        return result;
-    }
-
-    /*
-    inline CT azimuth21() const
-    {
-        // NOTE: signs of X and Y are different than in the original paper
-        return is_distance_zero ?
-               CT(0) :
-               atan2(-sin_alpha, sin_U1 * sin_sigma - cos_U1 * cos_sigma * cos_azimuth12); // (12)
-    }
-    */
-};
-
-}}} // namespace boost::geometry::detail
-
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_DIRECT_HPP
diff --git a/boost/geometry/algorithms/detail/vincenty_inverse.hpp b/boost/geometry/algorithms/detail/vincenty_inverse.hpp
deleted file mode 100644
index fe05e95932..0000000000
--- a/boost/geometry/algorithms/detail/vincenty_inverse.hpp
+++ /dev/null
@@ -1,204 +0,0 @@
-// Boost.Geometry
-
-// 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.
-
-// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
-
-// Use, modification and distribution is subject to the Boost Software License,
-// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
-// http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_INVERSE_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_INVERSE_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/algorithms/detail/flattening.hpp>
-#include <boost/geometry/algorithms/detail/result_inverse.hpp>
-
-
-#ifndef BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS
-#define BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS 1000
-#endif
-
-
-namespace boost { namespace geometry { namespace detail
-{
-
-/*!
-\brief The solution of the inverse problem of geodesics on latlong coordinates, after Vincenty, 1975
-\author See
-    - http://www.ngs.noaa.gov/PUBS_LIB/inverse.pdf
-    - http://www.icsm.gov.au/gda/gdav2.3.pdf
-\author Adapted from various implementations to get it close to the original document
-    - http://www.movable-type.co.uk/scripts/LatLongVincenty.html
-    - http://exogen.case.edu/projects/geopy/source/geopy.distance.html
-    - http://futureboy.homeip.net/fsp/colorize.fsp?fileName=navigation.frink
-
-*/
-template <typename CT, bool EnableDistance, bool EnableAzimuth>
-struct vincenty_inverse
-{
-    typedef result_inverse<CT> result_type;
-
-public:
-    template <typename T1, typename T2, typename Spheroid>
-    static inline result_type apply(T1 const& lon1,
-                                    T1 const& lat1,
-                                    T2 const& lon2,
-                                    T2 const& lat2,
-                                    Spheroid const& spheroid)
-    {
-        result_type result;
-
-        if (math::equals(lat1, lat2) && math::equals(lon1, lon2))
-        {
-            result.set(CT(0), CT(0));
-            return result;
-        }
-
-        CT const c1 = 1;
-        CT const c2 = 2;
-        CT const c3 = 3;
-        CT const c4 = 4;
-        CT const c16 = 16;
-        CT const c_e_12 = CT(1e-12);
-
-        CT const pi = geometry::math::pi<CT>();
-        CT const two_pi = c2 * pi;
-
-        // lambda: difference in longitude on an auxiliary sphere
-        CT L = lon2 - lon1;
-        CT lambda = L;
-
-        if (L < -pi) L += two_pi;
-        if (L > pi) L -= two_pi;
-
-        CT const radius_a = CT(get_radius<0>(spheroid));
-        CT const radius_b = CT(get_radius<2>(spheroid));
-        CT const flattening = geometry::detail::flattening<CT>(spheroid);
-
-        // U: reduced latitude, defined by tan U = (1-f) tan phi
-        CT const one_min_f = c1 - flattening;
-        CT const tan_U1 = one_min_f * tan(lat1); // above (1)
-        CT const tan_U2 = one_min_f * tan(lat2); // above (1)
-
-        // calculate sin U and cos U using trigonometric identities
-        CT const temp_den_U1 = math::sqrt(c1 + math::sqr(tan_U1));
-        CT const temp_den_U2 = math::sqrt(c1 + math::sqr(tan_U2));
-        // cos = 1 / sqrt(1 + tan^2)
-        CT const cos_U1 = c1 / temp_den_U1;
-        CT const cos_U2 = c1 / temp_den_U2;
-        // sin = tan / sqrt(1 + tan^2)
-        CT const sin_U1 = tan_U1 / temp_den_U1;
-        CT const sin_U2 = tan_U2 / temp_den_U2;
-
-        // calculate sin U and cos U directly
-        //CT const U1 = atan(tan_U1);
-        //CT const U2 = atan(tan_U2);
-        //cos_U1 = cos(U1);
-        //cos_U2 = cos(U2);
-        //sin_U1 = tan_U1 * cos_U1; // sin(U1);
-        //sin_U2 = tan_U2 * cos_U2; // sin(U2);
-
-        CT previous_lambda;
-        CT sin_lambda;
-        CT cos_lambda;
-        CT sin_sigma;
-        CT sin_alpha;
-        CT cos2_alpha;
-        CT cos2_sigma_m;
-        CT sigma;
-
-        int counter = 0; // robustness
-
-        do
-        {
-            previous_lambda = lambda; // (13)
-            sin_lambda = sin(lambda);
-            cos_lambda = cos(lambda);
-            sin_sigma = math::sqrt(math::sqr(cos_U2 * sin_lambda) + math::sqr(cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda)); // (14)
-            CT cos_sigma = sin_U1 * sin_U2 + cos_U1 * cos_U2 * cos_lambda; // (15)
-            sin_alpha = cos_U1 * cos_U2 * sin_lambda / sin_sigma; // (17)
-            cos2_alpha = c1 - math::sqr(sin_alpha);
-            cos2_sigma_m = math::equals(cos2_alpha, 0) ? 0 : cos_sigma - c2 * sin_U1 * sin_U2 / cos2_alpha; // (18)
-
-            CT C = flattening/c16 * cos2_alpha * (c4 + flattening * (c4 - c3 * cos2_alpha)); // (10)
-            sigma = atan2(sin_sigma, cos_sigma); // (16)
-            lambda = L + (c1 - C) * flattening * sin_alpha *
-                (sigma + C * sin_sigma * ( cos2_sigma_m + C * cos_sigma * (-c1 + c2 * math::sqr(cos2_sigma_m)))); // (11)
-
-            ++counter; // robustness
-
-        } while ( geometry::math::abs(previous_lambda - lambda) > c_e_12
-               && geometry::math::abs(lambda) < pi
-               && counter < BOOST_GEOMETRY_DETAIL_VINCENTY_MAX_STEPS ); // robustness
-    
-        if ( BOOST_GEOMETRY_CONDITION(EnableDistance) )
-        {
-            // Oops getting hard here
-            // (again, problem is that ttmath cannot divide by doubles, which is OK)
-            CT const c1 = 1;
-            CT const c2 = 2;
-            CT const c3 = 3;
-            CT const c4 = 4;
-            CT const c6 = 6;
-            CT const c47 = 47;
-            CT const c74 = 74;
-            CT const c128 = 128;
-            CT const c256 = 256;
-            CT const c175 = 175;
-            CT const c320 = 320;
-            CT const c768 = 768;
-            CT const c1024 = 1024;
-            CT const c4096 = 4096;
-            CT const c16384 = 16384;
-
-            //CT sqr_u = cos2_alpha * (math::sqr(radius_a) - math::sqr(radius_b)) / math::sqr(radius_b); // above (1)
-            CT sqr_u = cos2_alpha * ( math::sqr(radius_a / radius_b) - c1 ); // above (1)
-
-            CT A = c1 + sqr_u/c16384 * (c4096 + sqr_u * (-c768 + sqr_u * (c320 - c175 * sqr_u))); // (3)
-            CT B = sqr_u/c1024 * (c256 + sqr_u * ( -c128 + sqr_u * (c74 - c47 * sqr_u))); // (4)
-            CT delta_sigma = B * sin_sigma * ( cos2_sigma_m + (B/c4) * (cos(sigma)* (-c1 + c2 * cos2_sigma_m)
-                - (B/c6) * cos2_sigma_m * (-c3 + c4 * math::sqr(sin_sigma)) * (-c3 + c4 * cos2_sigma_m))); // (6)
-
-            result.distance = radius_b * A * (sigma - delta_sigma); // (19)
-        }
-        else
-        {
-            result.distance = CT(0);
-        }
-    
-        if ( BOOST_GEOMETRY_CONDITION(EnableAzimuth) )
-        {
-            result.azimuth = atan2(cos_U2 * sin_lambda, cos_U1 * sin_U2 - sin_U1 * cos_U2 * cos_lambda); // (20)
-        }
-        else
-        {
-            result.azimuth = CT(0);
-        }
-
-        return result;
-    }
-
-//    inline CT azimuth21() const
-//    {
-//        // NOTE: signs of X and Y are different than in the original paper
-//        atan2(-cos_U1 * sin_lambda, sin_U1 * cos_U2 - cos_U1 * sin_U2 * cos_lambda); // (21)
-//    }
-};
-
-}}} // namespace boost::geometry::detail
-
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_VINCENTY_INVERSE_HPP