Imported Upstream version 1.61.0

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

16 files changed, 2459 insertions, 1657 deletions

diff --git a/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp b/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp
index 285edfdd6c..03c06c28d1 100644
--- a/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp
+++ b/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp
@@ -72,7 +72,7 @@ inline void append_no_dups_or_spikes(Range& range, Point const& point,
         << geometry::get<0>(point) << ", " << geometry::get<1>(point) << ")"
         << std::endl;
 #endif
-    // The code below thies condition checks all spikes/dups
+    // The code below this condition checks all spikes/dups
     // for geometries >= 3 points.
     // So we have to check the first potential duplicate differently
     if (boost::size(range) == 1
diff --git a/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp b/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp
index 90901dee70..a8171e1482 100644
--- a/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp
+++ b/boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp
@@ -48,7 +48,6 @@ inline void clear_visit_info(Turns& turns)
         {
             op_it->visited.clear();
         }
-        it->discarded = false;
     }
 }
 
@@ -62,6 +61,34 @@ struct backtrack_state
 };
 
 
+enum traverse_error_type
+{
+    traverse_error_none,
+    traverse_error_no_next_ip_at_start,
+    traverse_error_no_next_ip,
+    traverse_error_dead_end_at_start,
+    traverse_error_dead_end,
+    traverse_error_visit_again,
+    traverse_error_endless_loop
+};
+
+inline std::string traverse_error_string(traverse_error_type error)
+{
+    switch (error)
+    {
+        case traverse_error_none : return "";
+        case traverse_error_no_next_ip_at_start : return "No next IP at start";
+        case traverse_error_no_next_ip : return "No next IP";
+        case traverse_error_dead_end_at_start : return "Dead end at start";
+        case traverse_error_dead_end : return "Dead end";
+        case traverse_error_visit_again : return "Visit again";
+        case traverse_error_endless_loop : return "Endless loop";
+        default : return "";
+    }
+    return "";
+}
+
+
 template
 <
     typename Geometry1,
@@ -73,23 +100,28 @@ class backtrack_check_self_intersections
     {
         bool m_checked;
         inline state()
-            : m_checked()
+            : m_checked(true)
         {}
     };
 public :
     typedef state state_type;
 
-    template <typename Operation, typename Rings, typename Ring, typename Turns, typename RobustPolicy>
+    template <typename Operation, typename Rings, typename Ring, typename Turns, typename RobustPolicy, typename Visitor>
     static inline void apply(std::size_t size_at_start,
                 Rings& rings, Ring& ring,
-                Turns& turns, Operation& operation,
-                std::string const& ,
+                Turns& turns,
+                typename boost::range_value<Turns>::type const& turn,
+                Operation& operation,
+                traverse_error_type traverse_error,
                 Geometry1 const& geometry1,
                 Geometry2 const& geometry2,
                 RobustPolicy const& robust_policy,
-                state_type& state
+                state_type& state,
+                Visitor& visitor
                 )
     {
+        visitor.visit_traverse_reject(turns, turn, operation, traverse_error);
+
         state.m_good = false;
 
         // Check self-intersections and throw exception if appropriate
diff --git a/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp b/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp
index dbf3770357..790819917c 100644
--- a/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp
+++ b/boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp
@@ -22,7 +22,9 @@
 #include <boost/geometry/algorithms/convert.hpp>
 #include <boost/geometry/geometries/concepts/check.hpp>
 #include <boost/geometry/util/range.hpp>
-#include <boost/geometry/views/detail/normalized_view.hpp>
+#include <boost/geometry/iterators/ever_circling_iterator.hpp>
+#include <boost/geometry/views/closeable_view.hpp>
+#include <boost/geometry/views/reversible_view.hpp>
 
 
 namespace boost { namespace geometry
@@ -38,25 +40,33 @@ template <typename Range, bool Reverse, typename SegmentIdentifier, typename Poi
 struct copy_segment_point_range
 {
     static inline bool apply(Range const& range,
-            SegmentIdentifier const& seg_id, bool second,
+            SegmentIdentifier const& seg_id, int offset,
             PointOut& point)
     {
-        detail::normalized_view<Range const> view(range);
+        typedef typename closeable_view
+        <
+            Range const,
+            closure<Range>::value
+        >::type cview_type;
 
-        signed_size_type const n = boost::size(view);
-        signed_size_type index = seg_id.segment_index;
-        if (second)
+        typedef typename reversible_view
+        <
+            cview_type const,
+            Reverse ? iterate_reverse : iterate_forward
+        >::type rview_type;
+
+        cview_type cview(range);
+        rview_type view(cview);
+
+        typedef typename boost::range_iterator<rview_type>::type iterator;
+        geometry::ever_circling_iterator<iterator> it(boost::begin(view), boost::end(view),
+                    boost::begin(view) + seg_id.segment_index, true);
+
+        for (signed_size_type i = 0; i < offset; ++i, ++it)
         {
-            index++;
-            if (index >= n)
-            {
-                index = 0;
-            }
         }
 
-        BOOST_GEOMETRY_ASSERT(index >= 0 && index < n);
-
-        geometry::convert(*(boost::begin(view) + index), point);
+        geometry::convert(*it, point);
         return true;
     }
 };
@@ -66,7 +76,7 @@ template <typename Polygon, bool Reverse, typename SegmentIdentifier, typename P
 struct copy_segment_point_polygon
 {
     static inline bool apply(Polygon const& polygon,
-            SegmentIdentifier const& seg_id, bool second,
+            SegmentIdentifier const& seg_id, int offset,
             PointOut& point)
     {
         // Call ring-version with the right ring
@@ -81,7 +91,7 @@ struct copy_segment_point_polygon
                     seg_id.ring_index < 0
                         ? geometry::exterior_ring(polygon)
                         : range::at(geometry::interior_rings(polygon), seg_id.ring_index),
-                    seg_id, second,
+                    seg_id, offset,
                     point
                 );
     }
@@ -92,11 +102,11 @@ template <typename Box, bool Reverse, typename SegmentIdentifier, typename Point
 struct copy_segment_point_box
 {
     static inline bool apply(Box const& box,
-            SegmentIdentifier const& seg_id, bool second,
+            SegmentIdentifier const& seg_id, int offset,
             PointOut& point)
     {
         signed_size_type index = seg_id.segment_index;
-        if (second)
+        for (int i = 0; i < offset; i++)
         {
             index++;
         }
@@ -119,7 +129,7 @@ template
 struct copy_segment_point_multi
 {
     static inline bool apply(MultiGeometry const& multi,
-                             SegmentIdentifier const& seg_id, bool second,
+                             SegmentIdentifier const& seg_id, int offset,
                              PointOut& point)
     {
 
@@ -130,7 +140,7 @@ struct copy_segment_point_multi
             );
 
         // Call the single-version
-        return Policy::apply(range::at(multi, seg_id.multi_index), seg_id, second, point);
+        return Policy::apply(range::at(multi, seg_id.multi_index), seg_id, offset, point);
     }
 };
 
@@ -271,7 +281,7 @@ struct copy_segment_point
  */
 template<bool Reverse, typename Geometry, typename SegmentIdentifier, typename PointOut>
 inline bool copy_segment_point(Geometry const& geometry,
-            SegmentIdentifier const& seg_id, bool second,
+            SegmentIdentifier const& seg_id, int offset,
             PointOut& point_out)
 {
     concept::check<Geometry const>();
@@ -283,7 +293,7 @@ inline bool copy_segment_point(Geometry const& geometry,
             Reverse,
             SegmentIdentifier,
             PointOut
-        >::apply(geometry, seg_id, second, point_out);
+        >::apply(geometry, seg_id, offset, point_out);
 }
 
 
@@ -300,7 +310,7 @@ template
     typename PointOut
 >
 inline bool copy_segment_point(Geometry1 const& geometry1, Geometry2 const& geometry2,
-            SegmentIdentifier const& seg_id, bool second,
+            SegmentIdentifier const& seg_id, int offset,
             PointOut& point_out)
 {
     concept::check<Geometry1 const>();
@@ -317,7 +327,7 @@ inline bool copy_segment_point(Geometry1 const& geometry1, Geometry2 const& geom
                 Reverse1,
                 SegmentIdentifier,
                 PointOut
-            >::apply(geometry1, seg_id, second, point_out);
+            >::apply(geometry1, seg_id, offset, point_out);
     }
     else if (seg_id.source_index == 1)
     {
@@ -328,7 +338,7 @@ inline bool copy_segment_point(Geometry1 const& geometry1, Geometry2 const& geom
                 Reverse2,
                 SegmentIdentifier,
                 PointOut
-            >::apply(geometry2, seg_id, second, point_out);
+            >::apply(geometry2, seg_id, offset, point_out);
     }
     // Exception?
     return false;
@@ -354,10 +364,33 @@ inline bool copy_segment_points(Geometry1 const& geometry1, Geometry2 const& geo
     concept::check<Geometry1 const>();
     concept::check<Geometry2 const>();
 
-    return copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, false, point1)
-        && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, true,  point2);
+    return copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 0, point1)
+        && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 1, point2);
 }
 
+/*!
+    \brief Helper function, copies three points: two from the specified segment
+    (from, to) and the next one
+    \ingroup overlay
+ */
+template
+<
+    bool Reverse1, bool Reverse2,
+    typename Geometry1, typename Geometry2,
+    typename SegmentIdentifier,
+    typename PointOut
+>
+inline bool copy_segment_points(Geometry1 const& geometry1, Geometry2 const& geometry2,
+            SegmentIdentifier const& seg_id,
+            PointOut& point1, PointOut& point2, PointOut& point3)
+{
+    concept::check<Geometry1 const>();
+    concept::check<Geometry2 const>();
+
+    return copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 0, point1)
+        && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 1, point2)
+        && copy_segment_point<Reverse1, Reverse2>(geometry1, geometry2, seg_id, 2, point3);
+}
 
 
 
diff --git a/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp b/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp
index bc84286241..c79c9cf500 100644
--- a/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp
+++ b/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp
@@ -24,17 +24,20 @@
 
 #include <boost/range.hpp>
 
+#include <boost/geometry/iterators/ever_circling_iterator.hpp>
 #include <boost/geometry/algorithms/detail/ring_identifier.hpp>
 #include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp>
 #include <boost/geometry/algorithms/detail/overlay/handle_colocations.hpp>
-#include <boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp>
+#include <boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp>
 #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp>
+#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
 #include <boost/geometry/policies/robustness/robust_type.hpp>
 #include <boost/geometry/strategies/side.hpp>
 #ifdef BOOST_GEOMETRY_DEBUG_ENRICH
 #  include <boost/geometry/algorithms/detail/overlay/check_enrich.hpp>
 #endif
 
+
 namespace boost { namespace geometry
 {
 
@@ -42,183 +45,6 @@ namespace boost { namespace geometry
 namespace detail { namespace overlay
 {
 
-// Wraps "turn_operation" from turn_info.hpp,
-// giving it extra information
-template <typename TurnOperation>
-struct indexed_turn_operation
-{
-    typedef TurnOperation type;
-
-    std::size_t turn_index;
-    std::size_t operation_index;
-    bool discarded;
-    // use pointers to avoid copies, const& is not possible because of usage in vector
-    segment_identifier const* other_seg_id; // segment id of other segment of intersection of two segments
-    TurnOperation const* subject;
-
-    inline indexed_turn_operation(std::size_t ti, std::size_t oi,
-                TurnOperation const& sub,
-                segment_identifier const& oid)
-        : turn_index(ti)
-        , operation_index(oi)
-        , discarded(false)
-        , other_seg_id(&oid)
-        , subject(boost::addressof(sub))
-    {}
-
-};
-
-template <typename IndexedTurnOperation>
-struct remove_discarded
-{
-    inline bool operator()(IndexedTurnOperation const& operation) const
-    {
-        return operation.discarded;
-    }
-};
-
-
-template
-<
-    typename TurnPoints,
-    typename Indexed,
-    typename Geometry1, typename Geometry2,
-    typename RobustPolicy,
-    bool Reverse1, bool Reverse2,
-    typename Strategy
->
-struct sort_on_segment_and_ratio
-{
-    inline sort_on_segment_and_ratio(TurnPoints const& turn_points
-            , Geometry1 const& geometry1
-            , Geometry2 const& geometry2
-            , RobustPolicy const& robust_policy
-            , Strategy const& strategy
-            , bool* clustered)
-        : m_turn_points(turn_points)
-        , m_geometry1(geometry1)
-        , m_geometry2(geometry2)
-        , m_robust_policy(robust_policy)
-        , m_strategy(strategy)
-        , m_clustered(clustered)
-    {
-    }
-
-private :
-
-    TurnPoints const& m_turn_points;
-    Geometry1 const& m_geometry1;
-    Geometry2 const& m_geometry2;
-    RobustPolicy const& m_robust_policy;
-    Strategy const& m_strategy;
-    mutable bool* m_clustered;
-
-    typedef typename geometry::point_type<Geometry1>::type point_type;
-
-    inline bool default_order(Indexed const& left, Indexed const& right) const
-    {
-        // We've nothing to sort on. Take the indexes
-        return left.turn_index < right.turn_index;
-    }
-
-    inline bool consider_relative_order(Indexed const& left,
-                    Indexed const& right) const
-    {
-        point_type pi, pj, ri, rj, si, sj;
-
-        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
-            left.subject->seg_id,
-            pi, pj);
-        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
-            *left.other_seg_id,
-            ri, rj);
-        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
-            *right.other_seg_id,
-            si, sj);
-
-        typedef typename strategy::side::services::default_strategy
-            <
-                typename cs_tag<point_type>::type
-            >::type strategy;
-
-        int const side_rj_p = strategy::apply(pi, pj, rj);
-        int const side_sj_p = strategy::apply(pi, pj, sj);
-
-        // Put the one turning left (1; right == -1) as last
-        if (side_rj_p != side_sj_p)
-        {
-            return side_rj_p < side_sj_p;
-        }
-
-        int const side_sj_r = strategy::apply(ri, rj, sj);
-        int const side_rj_s = strategy::apply(si, sj, rj);
-
-        // If they both turn left: the most left as last
-        // If they both turn right: this is not relevant, but take also here most left
-        if (side_rj_s != side_sj_r)
-        {
-            return side_rj_s < side_sj_r;
-        }
-
-        return default_order(left, right);
-    }
-
-public :
-
-    // Note that left/right do NOT correspond to m_geometry1/m_geometry2
-    // but to the "indexed_turn_operation"
-    inline bool operator()(Indexed const& left, Indexed const& right) const
-    {
-        if (! (left.subject->seg_id == right.subject->seg_id))
-        {
-            return left.subject->seg_id < right.subject->seg_id;
-        }
-
-        // Both left and right are located on the SAME segment.
-
-        if (! (left.subject->fraction == right.subject->fraction))
-        {
-            return left.subject->fraction < right.subject->fraction;
-        }
-
-
-        // First check "real" intersection (crosses)
-        // -> distance zero due to precision, solve it by sorting
-        if (m_turn_points[left.turn_index].method == method_crosses
-            && m_turn_points[right.turn_index].method == method_crosses)
-        {
-            return consider_relative_order(left, right);
-        }
-
-        // If that is not the case, cluster it later on.
-        // Indicate that this is necessary.
-        *m_clustered = true;
-
-        return default_order(left, right);
-    }
-};
-
-
-template<typename Turns, typename Operations>
-inline void update_discarded(Turns& turn_points, Operations& operations)
-{
-    // Vice-versa, set discarded to true for discarded operations;
-    // AND set discarded points to true
-    for (typename boost::range_iterator<Operations>::type it = boost::begin(operations);
-         it != boost::end(operations);
-         ++it)
-    {
-        if (turn_points[it->turn_index].discarded)
-        {
-            it->discarded = true;
-        }
-        else if (it->discarded)
-        {
-            turn_points[it->turn_index].discarded = true;
-        }
-    }
-}
-
 
 // Sorts IP-s of this ring on segment-identifier, and if on same segment,
 //  on distance.
@@ -227,148 +53,88 @@ inline void update_discarded(Turns& turn_points, Operations& operations)
 // (might be on another segment)
 template
 <
-    typename IndexType,
     bool Reverse1, bool Reverse2,
-    typename Container,
-    typename TurnPoints,
+    typename Operations,
+    typename Turns,
     typename Geometry1, typename Geometry2,
     typename RobustPolicy,
     typename Strategy
 >
-inline void enrich_sort(Container& operations,
-            TurnPoints& turn_points,
+inline void enrich_sort(Operations& operations,
+            Turns const& turns,
             operation_type for_operation,
-            Geometry1 const& geometry1, Geometry2 const& geometry2,
+            Geometry1 const& geometry1,
+            Geometry2 const& geometry2,
             RobustPolicy const& robust_policy,
             Strategy const& strategy)
 {
-    typedef typename IndexType::type operations_type;
-
-    bool clustered = false;
     std::sort(boost::begin(operations),
-                boost::end(operations),
-                sort_on_segment_and_ratio
-                    <
-                        TurnPoints,
-                        IndexType,
-                        Geometry1, Geometry2,
-                        RobustPolicy,
-                        Reverse1, Reverse2,
-                        Strategy
-                    >(turn_points, geometry1, geometry2, robust_policy, strategy, &clustered));
-
-    // DONT'T discard xx / (for union) ix / ii / (for intersection) ux / uu here
-    // It would give way to "lonely" ui turn points, traveling all
-    // the way round. See #105
-
-    if (clustered)
-    {
-        typedef typename boost::range_iterator<Container>::type nc_iterator;
-        nc_iterator it = boost::begin(operations);
-        nc_iterator begin_cluster = boost::end(operations);
-        for (nc_iterator prev = it++;
-            it != boost::end(operations);
-            prev = it++)
-        {
-            operations_type& prev_op = turn_points[prev->turn_index]
-                .operations[prev->operation_index];
-            operations_type& op = turn_points[it->turn_index]
-                .operations[it->operation_index];
-
-            if (prev_op.seg_id == op.seg_id
-                && (turn_points[prev->turn_index].method != method_crosses
-                    || turn_points[it->turn_index].method != method_crosses)
-                && prev_op.fraction == op.fraction)
-            {
-                if (begin_cluster == boost::end(operations))
-                {
-                    begin_cluster = prev;
-                }
-            }
-            else if (begin_cluster != boost::end(operations))
-            {
-                handle_cluster<IndexType, Reverse1, Reverse2>(begin_cluster, it, turn_points,
-                        for_operation, geometry1, geometry2, robust_policy, strategy);
-                begin_cluster = boost::end(operations);
-            }
-        }
-        if (begin_cluster != boost::end(operations))
-        {
-            handle_cluster<IndexType, Reverse1, Reverse2>(begin_cluster, it, turn_points,
-                    for_operation, geometry1, geometry2, robust_policy, strategy);
-        }
-    }
-
-    update_discarded(turn_points, operations);
+            boost::end(operations),
+            less_by_segment_ratio
+                <
+                    Turns,
+                    typename boost::range_value<Operations>::type,
+                    Geometry1, Geometry2,
+                    RobustPolicy,
+                    Reverse1, Reverse2
+                >(turns, for_operation, geometry1, geometry2, robust_policy));
 }
 
 
-template
-<
-    typename IndexType,
-    typename Container,
-    typename TurnPoints
->
-inline void enrich_discard(Container& operations, TurnPoints& turn_points)
-{
-    update_discarded(turn_points, operations);
-
-    // Then delete discarded operations from vector
-    remove_discarded<IndexType> predicate;
-    operations.erase(
-            std::remove_if(boost::begin(operations),
-                    boost::end(operations),
-                    predicate),
-            boost::end(operations));
-}
-
-template
-<
-    typename IndexType,
-    typename Container,
-    typename TurnPoints,
-    typename Geometry1,
-    typename Geometry2,
-    typename Strategy
->
-inline void enrich_assign(Container& operations,
-            TurnPoints& turn_points,
-            operation_type ,
-            Geometry1 const& , Geometry2 const& ,
-            Strategy const& )
+template <typename Operations, typename Turns>
+inline void enrich_assign(Operations& operations, Turns& turns)
 {
-    typedef typename IndexType::type operations_type;
-    typedef typename boost::range_iterator<Container const>::type iterator_type;
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::turn_operation_type op_type;
+    typedef typename boost::range_iterator<Operations>::type iterator_type;
 
 
     if (operations.size() > 0)
     {
         // Assign travel-to-vertex/ip index for each turning point.
-        // Because IP's are circular, PREV starts at the very last one,
-        // being assigned from the first one.
-        // "next ip on same segment" should not be considered circular.
-        bool first = true;
-        iterator_type it = boost::begin(operations);
-        for (iterator_type prev = it + (boost::size(operations) - 1);
-             it != boost::end(operations);
-             prev = it++)
+        // Iterator "next" is circular
+
+        geometry::ever_circling_range_iterator<Operations const> next(operations);
+        ++next;
+
+        for (iterator_type it = boost::begin(operations);
+             it != boost::end(operations); ++it)
         {
-            operations_type& prev_op
-                    = turn_points[prev->turn_index].operations[prev->operation_index];
-            operations_type& op
-                    = turn_points[it->turn_index].operations[it->operation_index];
-
-            prev_op.enriched.travels_to_ip_index
-                    = static_cast<signed_size_type>(it->turn_index);
-            prev_op.enriched.travels_to_vertex_index
-                    = it->subject->seg_id.segment_index;
-
-            if (! first
-                && prev_op.seg_id.segment_index == op.seg_id.segment_index)
+            turn_type& turn = turns[it->turn_index];
+            op_type& op = turn.operations[it->operation_index];
+
+            // Normal behaviour: next should point at next turn:
+            if (it->turn_index == next->turn_index)
+            {
+                ++next;
+            }
+
+            // Cluster behaviour: next should point after cluster, unless
+            // their seg_ids are not the same
+            while (turn.cluster_id != -1
+                   && it->turn_index != next->turn_index
+                   && turn.cluster_id == turns[next->turn_index].cluster_id
+                   && op.seg_id == turns[next->turn_index].operations[next->operation_index].seg_id)
+            {
+                ++next;
+            }
+
+            turn_type const& next_turn = turns[next->turn_index];
+            op_type const& next_op = next_turn.operations[next->operation_index];
+
+            op.enriched.travels_to_ip_index
+                    = static_cast<signed_size_type>(next->turn_index);
+            op.enriched.travels_to_vertex_index
+                    = next->subject->seg_id.segment_index;
+
+            if (op.seg_id.segment_index == next_op.seg_id.segment_index
+                    && op.fraction < next_op.fraction)
             {
-                prev_op.enriched.next_ip_index = static_cast<signed_size_type>(it->turn_index);
+                // Next turn is located further on same segment
+                // assign next_ip_index
+                // (this is one not circular therefore fraction is considered)
+                op.enriched.next_ip_index = static_cast<signed_size_type>(next->turn_index);
             }
-            first = false;
         }
     }
 
@@ -379,20 +145,20 @@ inline void enrich_assign(Container& operations,
              it != boost::end(operations);
              ++it)
         {
-            operations_type& op = turn_points[it->turn_index]
+            op_type& op = turns[it->turn_index]
                 .operations[it->operation_index];
 
             std::cout << it->turn_index
-                << " meth: " << method_char(turn_points[it->turn_index].method)
-                << " seg: " << op.seg_id
-                << " dst: " << op.fraction // needs define
-                << " op: " << operation_char(turn_points[it->turn_index].operations[0].operation)
-                << operation_char(turn_points[it->turn_index].operations[1].operation)
-                << " dsc: " << (turn_points[it->turn_index].discarded ? "T" : "F")
-                << " ->vtx " << op.enriched.travels_to_vertex_index
-                << " ->ip " << op.enriched.travels_to_ip_index
-                << " ->nxt ip " << op.enriched.next_ip_index
-                //<< " vis: " << visited_char(op.visited)
+                << " cl=" << turns[it->turn_index].cluster_id
+                << " meth=" << method_char(turns[it->turn_index].method)
+                << " seg=" << op.seg_id
+                << " dst=" << op.fraction // needs define
+                << " op=" << operation_char(turns[it->turn_index].operations[0].operation)
+                << operation_char(turns[it->turn_index].operations[1].operation)
+                << " (" << operation_char(op.operation) << ")"
+                << " nxt=" << op.enriched.next_ip_index
+                << " / " << op.enriched.travels_to_ip_index
+                << " [vx " << op.enriched.travels_to_vertex_index << "]"
                 << std::endl;
                 ;
         }
@@ -403,43 +169,57 @@ inline void enrich_assign(Container& operations,
 }
 
 
-template <typename IndexedType, typename TurnPoints, typename MappedVector>
-inline void create_map(TurnPoints const& turn_points, MappedVector& mapped_vector)
+template <typename Turns, typename MappedVector>
+inline void create_map(Turns const& turns,
+        detail::overlay::operation_type for_operation,
+        MappedVector& mapped_vector)
 {
-    typedef typename boost::range_value<TurnPoints>::type turn_point_type;
-    typedef typename turn_point_type::container_type container_type;
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::container_type container_type;
+    typedef typename MappedVector::mapped_type mapped_type;
+    typedef typename boost::range_value<mapped_type>::type indexed_type;
 
     std::size_t index = 0;
-    for (typename boost::range_iterator<TurnPoints const>::type
-            it = boost::begin(turn_points);
-         it != boost::end(turn_points);
+    for (typename boost::range_iterator<Turns const>::type
+            it = boost::begin(turns);
+         it != boost::end(turns);
          ++it, ++index)
     {
-        // Add operations on this ring, but skip discarded ones
-        if (! it->discarded)
+        // Add all (non discarded) operations on this ring
+        // Blocked operations or uu on clusters (for intersection)
+        // should be included, to block potential paths in clusters
+        turn_type const& turn = *it;
+        if (turn.discarded)
         {
-            std::size_t op_index = 0;
-            for (typename boost::range_iterator<container_type const>::type
-                    op_it = boost::begin(it->operations);
-                op_it != boost::end(it->operations);
-                ++op_it, ++op_index)
-            {
-                // We should NOT skip blocked operations here
-                // because they can be relevant for "the other side"
-                // NOT if (op_it->operation != operation_blocked)
-
-                ring_identifier ring_id
-                    (
-                        op_it->seg_id.source_index,
-                        op_it->seg_id.multi_index,
-                        op_it->seg_id.ring_index
-                    );
-                mapped_vector[ring_id].push_back
-                    (
-                        IndexedType(index, op_index, *op_it,
-                            it->operations[1 - op_index].seg_id)
-                    );
-            }
+            continue;
+        }
+
+        if (for_operation == operation_intersection
+            && turn.cluster_id == -1
+            && turn.both(operation_union))
+        {
+            // Only include uu turns if part of cluster (to block potential paths),
+            // otherwise they can block possibly viable paths
+            continue;
+        }
+
+        std::size_t op_index = 0;
+        for (typename boost::range_iterator<container_type const>::type
+                op_it = boost::begin(turn.operations);
+            op_it != boost::end(turn.operations);
+            ++op_it, ++op_index)
+        {
+            ring_identifier const ring_id
+                (
+                    op_it->seg_id.source_index,
+                    op_it->seg_id.multi_index,
+                    op_it->seg_id.ring_index
+                );
+            mapped_vector[ring_id].push_back
+                (
+                    indexed_type(index, op_index, *op_it,
+                        it->operations[1 - op_index].seg_id)
+                );
         }
     }
 }
@@ -453,12 +233,12 @@ inline void create_map(TurnPoints const& turn_points, MappedVector& mapped_vecto
 /*!
 \brief All intersection points are enriched with successor information
 \ingroup overlay
-\tparam TurnPoints type of intersection container
+\tparam Turns type of intersection container
             (e.g. vector of "intersection/turn point"'s)
 \tparam Geometry1 \tparam_geometry
 \tparam Geometry2 \tparam_geometry
 \tparam Strategy side strategy type
-\param turn_points container containing intersectionpoints
+\param turns container containing intersectionpoints
 \param for_operation operation_type (union or intersection)
 \param geometry1 \param_geometry
 \param geometry2 \param_geometry
@@ -469,22 +249,24 @@ template
 <
     bool Reverse1, bool Reverse2,
     overlay_type OverlayType,
-    typename TurnPoints,
+    typename Turns,
+    typename Clusters,
     typename Geometry1, typename Geometry2,
     typename RobustPolicy,
     typename Strategy
 >
-inline void enrich_intersection_points(TurnPoints& turn_points,
+inline void enrich_intersection_points(Turns& turns,
+    Clusters& clusters,
     detail::overlay::operation_type for_operation,
     Geometry1 const& geometry1, Geometry2 const& geometry2,
     RobustPolicy const& robust_policy,
     Strategy const& strategy)
 {
-    typedef typename boost::range_value<TurnPoints>::type turn_point_type;
-    typedef typename turn_point_type::turn_operation_type turn_operation_type;
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::turn_operation_type op_type;
     typedef detail::overlay::indexed_turn_operation
         <
-            turn_operation_type
+            op_type
         > indexed_turn_operation;
 
     typedef std::map
@@ -493,51 +275,27 @@ inline void enrich_intersection_points(TurnPoints& turn_points,
             std::vector<indexed_turn_operation>
         > mapped_vector_type;
 
-    // Iterate through turns and discard uu
-    // and check if there are possible colocations
-    bool check_colocations = false;
-    for (typename boost::range_iterator<TurnPoints>::type
-            it = boost::begin(turn_points);
-         it != boost::end(turn_points);
+    bool const has_colocations
+        = detail::overlay::handle_colocations<Reverse1, Reverse2>(turns,
+        clusters, geometry1, geometry2);
+
+    // Discard none turns, if any
+    for (typename boost::range_iterator<Turns>::type
+            it = boost::begin(turns);
+         it != boost::end(turns);
          ++it)
     {
-        if (it->both(detail::overlay::operation_union))
-        {
-            // Discard  (necessary for a.o. #76). With uu, at all points there
-            // is the risk that rings are being traversed twice or more.
-            // Without uu, all rings having only uu will be untouched
-            // and gathered by assemble
-            it->discarded = true;
-            check_colocations = true;
-        }
-        else if (it->combination(detail::overlay::operation_union,
-                                 detail::overlay::operation_blocked))
-        {
-            check_colocations = true;
-        }
-        else if (OverlayType == overlay_difference
-                 && it->both(detail::overlay::operation_intersection))
-        {
-            // For difference operation (u/u -> i/i)
-            check_colocations = true;
-        }
-        else if (it->both(detail::overlay::operation_none))
+        if (it->both(detail::overlay::operation_none))
         {
             it->discarded = true;
         }
     }
 
-    if (check_colocations)
-    {
-        detail::overlay::handle_colocations<OverlayType>(turn_points);
-    }
-
     // Create a map of vectors of indexed operation-types to be able
     // to sort intersection points PER RING
     mapped_vector_type mapped_vector;
 
-    detail::overlay::create_map<indexed_turn_operation>(turn_points, mapped_vector);
-
+    detail::overlay::create_map(turns, for_operation, mapped_vector);
 
     // No const-iterator; contents of mapped copy is temporary,
     // and changed by enrich
@@ -550,8 +308,10 @@ inline void enrich_intersection_points(TurnPoints& turn_points,
     std::cout << "ENRICH-sort Ring "
         << mit->first << std::endl;
 #endif
-        detail::overlay::enrich_sort<indexed_turn_operation, Reverse1, Reverse2>(mit->second, turn_points, for_operation,
-                    geometry1, geometry2, robust_policy, strategy);
+        detail::overlay::enrich_sort<Reverse1, Reverse2>(
+                    mit->second, turns, for_operation,
+                    geometry1, geometry2,
+                    robust_policy, strategy);
     }
 
     for (typename mapped_vector_type::iterator mit
@@ -560,27 +320,20 @@ inline void enrich_intersection_points(TurnPoints& turn_points,
         ++mit)
     {
 #ifdef BOOST_GEOMETRY_DEBUG_ENRICH
-    std::cout << "ENRICH-discard Ring "
+    std::cout << "ENRICH-assign Ring "
         << mit->first << std::endl;
 #endif
-        detail::overlay::enrich_discard<indexed_turn_operation>(mit->second, turn_points);
+        detail::overlay::enrich_assign(mit->second, turns);
     }
 
-    for (typename mapped_vector_type::iterator mit
-        = mapped_vector.begin();
-        mit != mapped_vector.end();
-        ++mit)
+    if (has_colocations)
     {
-#ifdef BOOST_GEOMETRY_DEBUG_ENRICH
-    std::cout << "ENRICH-assign Ring "
-        << mit->first << std::endl;
-#endif
-        detail::overlay::enrich_assign<indexed_turn_operation>(mit->second, turn_points, for_operation,
-                    geometry1, geometry2, strategy);
+        detail::overlay::assign_startable_in_clusters<Reverse1, Reverse2>(
+                clusters, turns, for_operation, geometry1, geometry2);
     }
 
 #ifdef BOOST_GEOMETRY_DEBUG_ENRICH
-    //detail::overlay::check_graph(turn_points, for_operation);
+    //detail::overlay::check_graph(turns, for_operation);
 #endif
 
 }
diff --git a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
index 5964f3ba44..f5c5cc6a2e 100644
--- a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
@@ -32,6 +32,9 @@ struct enrichment_info
         : travels_to_vertex_index(-1)
         , travels_to_ip_index(-1)
         , next_ip_index(-1)
+        , startable(true)
+        , count_left(0)
+        , count_right(0)
     {}
 
     // vertex to which is free travel after this IP,
@@ -44,6 +47,12 @@ struct enrichment_info
 
     // index of next IP on this segment, -1 if there is no one
     signed_size_type next_ip_index;
+
+    bool startable; // Can be used to start in traverse
+
+    // Counts if polygons left/right of this operation
+    std::size_t count_left;
+    std::size_t count_right;
 };
 
 
diff --git a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
index 6f332ddff2..ad75127fab 100644
--- a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
+++ b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
@@ -16,6 +16,7 @@
 
 #include <boost/range.hpp>
 #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp>
+#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
 #include <boost/geometry/algorithms/detail/ring_identifier.hpp>
 #include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp>
@@ -34,6 +35,29 @@ namespace boost { namespace geometry
 namespace detail { namespace overlay
 {
 
+template <typename SegmentRatio>
+struct segment_fraction
+{
+    segment_identifier seg_id;
+    SegmentRatio fraction;
+
+    segment_fraction(segment_identifier const& id, SegmentRatio const& fr)
+        : seg_id(id)
+        , fraction(fr)
+    {}
+
+    segment_fraction()
+    {}
+
+    bool operator<(segment_fraction<SegmentRatio> const& other) const
+    {
+        return seg_id == other.seg_id
+                ? fraction < other.fraction
+                : seg_id < other.seg_id;
+    }
+
+};
+
 struct turn_operation_index
 {
     turn_operation_index(signed_size_type ti = -1,
@@ -43,22 +67,22 @@ struct turn_operation_index
     {}
 
     signed_size_type turn_index;
-    signed_size_type op_index; // basically only 0,1
+    signed_size_type op_index; // only 0,1
 };
 
 
-template <typename TurnPoints>
+template <typename Turns>
 struct less_by_fraction_and_type
 {
-    inline less_by_fraction_and_type(TurnPoints const& turn_points)
-        : m_turns(turn_points)
+    inline less_by_fraction_and_type(Turns const& turns)
+        : m_turns(turns)
     {
     }
 
     inline bool operator()(turn_operation_index const& left,
                            turn_operation_index const& right) const
     {
-        typedef typename boost::range_value<TurnPoints>::type turn_type;
+        typedef typename boost::range_value<Turns>::type turn_type;
         typedef typename turn_type::turn_operation_type turn_operation_type;
 
         turn_type const& left_turn = m_turns[left.turn_index];
@@ -74,6 +98,29 @@ struct less_by_fraction_and_type
             return left_op.fraction < right_op.fraction;
         }
 
+        // Order xx first - used to discard any following colocated turn
+        bool const left_both_xx = left_turn.both(operation_blocked);
+        bool const right_both_xx = right_turn.both(operation_blocked);
+        if (left_both_xx && ! right_both_xx)
+        {
+            return true;
+        }
+        if (! left_both_xx && right_both_xx)
+        {
+            return false;
+        }
+
+        bool const left_both_uu = left_turn.both(operation_union);
+        bool const right_both_uu = right_turn.both(operation_union);
+        if (left_both_uu && ! right_both_uu)
+        {
+            return true;
+        }
+        if (! left_both_uu && right_both_uu)
+        {
+            return false;
+        }
+
         turn_operation_type const& left_other_op
                 = left_turn.operations[1 - left.op_index];
 
@@ -86,88 +133,219 @@ struct less_by_fraction_and_type
     }
 
 private:
-    TurnPoints const& m_turns;
+    Turns const& m_turns;
 };
 
-template <overlay_type OverlayType, typename TurnPoints, typename OperationVector>
-inline void handle_colocation_cluster(TurnPoints& turn_points,
-       segment_identifier const& current_ring_seg_id,
-       OperationVector const& vec)
+template <typename Operation, typename ClusterPerSegment>
+inline signed_size_type get_cluster_id(Operation const& op, ClusterPerSegment const& cluster_per_segment)
+{
+    typedef typename ClusterPerSegment::key_type segment_fraction_type;
+
+    segment_fraction_type seg_frac(op.seg_id, op.fraction);
+    typename ClusterPerSegment::const_iterator it
+            = cluster_per_segment.find(seg_frac);
+
+    if (it == cluster_per_segment.end())
+    {
+        return -1;
+    }
+    return it->second;
+}
+
+template <typename Operation, typename ClusterPerSegment>
+inline void add_cluster_id(Operation const& op,
+    ClusterPerSegment& cluster_per_segment, signed_size_type id)
 {
-    typedef typename boost::range_value<TurnPoints>::type turn_type;
+    typedef typename ClusterPerSegment::key_type segment_fraction_type;
+
+    segment_fraction_type seg_frac(op.seg_id, op.fraction);
+
+    cluster_per_segment[seg_frac] = id;
+}
+
+template <typename Turn, typename ClusterPerSegment>
+inline signed_size_type add_turn_to_cluster(Turn const& turn,
+        ClusterPerSegment& cluster_per_segment, signed_size_type& cluster_id)
+{
+    signed_size_type cid0 = get_cluster_id(turn.operations[0], cluster_per_segment);
+    signed_size_type cid1 = get_cluster_id(turn.operations[1], cluster_per_segment);
+
+    if (cid0 == -1 && cid1 == -1)
+    {
+        ++cluster_id;
+        add_cluster_id(turn.operations[0], cluster_per_segment, cluster_id);
+        add_cluster_id(turn.operations[1], cluster_per_segment, cluster_id);
+        return cluster_id;
+    }
+    else if (cid0 == -1 && cid1 != -1)
+    {
+        add_cluster_id(turn.operations[0], cluster_per_segment, cid1);
+        return cid1;
+    }
+    else if (cid0 != -1 && cid1 == -1)
+    {
+        add_cluster_id(turn.operations[1], cluster_per_segment, cid0);
+        return cid0;
+    }
+    else if (cid0 == cid1)
+    {
+        // Both already added to same cluster, no action
+        return cid0;
+    }
+
+    // Both operations.seg_id/fraction were already part of any cluster, and
+    // these clusters are not the same. Merge of two clusters is necessary
+    std::cout << " TODO: merge " << cid0 << " and " << cid1 << std::endl;
+    return cid0;
+}
+
+template
+<
+    bool Reverse1, bool Reverse2,
+    typename Turns,
+    typename ClusterPerSegment,
+    typename Operations,
+    typename Geometry1,
+    typename Geometry2
+>
+inline void handle_colocation_cluster(Turns& turns,
+        signed_size_type& cluster_id,
+        ClusterPerSegment& cluster_per_segment,
+        Operations const& operations,
+        Geometry1 const& geometry1, Geometry2 const& geometry2)
+{
+    typedef typename boost::range_value<Turns>::type turn_type;
     typedef typename turn_type::turn_operation_type turn_operation_type;
 
-    std::vector<turn_operation_index>::const_iterator vit = vec.begin();
+    std::vector<turn_operation_index>::const_iterator vit = operations.begin();
 
-    turn_type cluster_turn = turn_points[vit->turn_index];
-    turn_operation_type cluster_op
-            = cluster_turn.operations[vit->op_index];
-    segment_identifier cluster_other_id
-            = cluster_turn.operations[1 - vit->op_index].seg_id;
-    bool const discard_colocated
-            = cluster_turn.both(operation_union)
-            || cluster_turn.combination(operation_blocked, operation_union);
+    turn_operation_index ref_toi = *vit;
+    signed_size_type ref_id = -1;
 
-    for (++vit; vit != vec.end(); ++vit)
+    for (++vit; vit != operations.end(); ++vit)
     {
+        turn_type& ref_turn = turns[ref_toi.turn_index];
+        turn_operation_type const& ref_op
+                = ref_turn.operations[ref_toi.op_index];
+
         turn_operation_index const& toi = *vit;
-        turn_type& turn = turn_points[toi.turn_index];
+        turn_type& turn = turns[toi.turn_index];
         turn_operation_type const& op = turn.operations[toi.op_index];
-        segment_identifier const& other_id
-                = turn.operations[1 - toi.op_index].seg_id;
 
-        if (cluster_op.fraction == op.fraction)
+        BOOST_ASSERT(ref_op.seg_id == op.seg_id);
+
+        if (ref_op.fraction == op.fraction)
         {
-            // Two turns of current ring with same source are colocated,
-            // one is from exterior ring, one from interior ring
-            bool const colocated_ext_int
-                = cluster_other_id.multi_index == other_id.multi_index
-                   && cluster_other_id.ring_index == -1
-                   && other_id.ring_index >= 0;
-
-            // Turn of current interior ring with other interior ring
-            bool const touch_int_int
-                = current_ring_seg_id.ring_index >= 0
-                   && other_id.ring_index >= 0;
-
-            if (discard_colocated && colocated_ext_int)
+            turn_operation_type const& other_op = turn.operations[1 - toi.op_index];
+
+            if (ref_id == -1)
             {
-                // If the two turns on this same segment are a
-                // colocation with two different segments on the
-                // other geometry, of the same polygon but with
-                // the outer (u/u or u/x) and the inner ring (non u/u),
-                // that turn with inner ring should be discarded
-                turn.discarded = true;
-                turn.colocated = true;
+                ref_id = add_turn_to_cluster(ref_turn, cluster_per_segment, cluster_id);
             }
-            else if (cluster_turn.colocated
-                     && touch_int_int
-                     && turn.both(operation_intersection))
+            BOOST_ASSERT(ref_id != -1);
+
+            // ref_turn (both operations) are already added to cluster,
+            // so also "op" is already added to cluster,
+            // We only need to add other_op
+            signed_size_type id = get_cluster_id(other_op, cluster_per_segment);
+            if (id != -1 && id != ref_id)
             {
-                // Two holes touch each other at a point where the
-                // exterior ring also touches
-                turn.discarded = true;
-                turn.colocated = true;
             }
-            else if (OverlayType == overlay_difference
-                     && turn.both(operation_intersection)
-                     && colocated_ext_int)
+            else if (id == -1)
+            {
+                // Add to same cluster
+                add_cluster_id(other_op, cluster_per_segment, ref_id);
+                id = ref_id;
+            }
+
+            // In case of colocated xx turns, all other turns may NOT be
+            // followed at all. xx cannot be discarded (otherwise colocated
+            // turns are followed).
+            if (ref_turn.both(operation_blocked))
             {
-                // For difference (polygon inside out) we need to
-                // discard i/i instead, in case of colocations
                 turn.discarded = true;
-                turn.colocated = true;
+                // We can either set or not set colocated because it is not effective on blocked turns
             }
         }
         else
         {
             // Not on same fraction on this segment
-            // assign for next potential cluster
-            cluster_turn = turn;
-            cluster_op = op;
-            cluster_other_id = other_id;
+            // assign for next
+            ref_toi = toi;
+            ref_id = -1;
         }
+    }
+}
+
+template
+<
+    typename Turns,
+    typename Clusters,
+    typename ClusterPerSegment
+>
+inline void assign_cluster_to_turns(Turns& turns,
+        Clusters& clusters,
+        ClusterPerSegment const& cluster_per_segment)
+{
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::turn_operation_type turn_operation_type;
+    typedef typename ClusterPerSegment::key_type segment_fraction_type;
 
+    signed_size_type turn_index = 0;
+    for (typename boost::range_iterator<Turns>::type it = turns.begin();
+         it != turns.end(); ++it, ++turn_index)
+    {
+        turn_type& turn = *it;
+
+        if (turn.discarded)
+        {
+            // They were processed (to create proper map) but will not be added
+            // This might leave a cluster with only 1 turn, which will be fixed
+            // afterwards
+            continue;
+        }
+
+        for (int i = 0; i < 2; i++)
+        {
+            turn_operation_type const& op = turn.operations[i];
+            segment_fraction_type seg_frac(op.seg_id, op.fraction);
+            typename ClusterPerSegment::const_iterator it = cluster_per_segment.find(seg_frac);
+            if (it != cluster_per_segment.end())
+            {
+                if (turn.cluster_id != -1
+                        && turn.cluster_id != it->second)
+                {
+                    std::cout << " CONFLICT " << std::endl;
+                }
+                turn.cluster_id = it->second;
+                clusters[turn.cluster_id].insert(turn_index);
+            }
+        }
+    }
+}
+
+template
+<
+    typename Turns,
+    typename Clusters
+>
+inline void remove_clusters(Turns& turns, Clusters& clusters)
+{
+    typename Clusters::iterator it = clusters.begin();
+    while (it != clusters.end())
+    {
+        // Hold iterator and increase. We can erase cit, this keeps the
+        // iterator valid (cf The standard associative-container erase idiom)
+        typename Clusters::iterator current_it = it;
+        ++it;
+
+        std::set<signed_size_type> const& turn_indices = current_it->second;
+        if (turn_indices.size() == 1)
+        {
+            signed_size_type turn_index = *turn_indices.begin();
+            turns[turn_index].cluster_id = -1;
+            clusters.erase(current_it);
+        }
     }
 }
 
@@ -179,8 +357,16 @@ inline void handle_colocation_cluster(TurnPoints& turn_points,
 // This function can be extended to replace handle_tangencies: at each
 // colocation incoming and outgoing vectors should be inspected
 
-template <overlay_type OverlayType, typename TurnPoints>
-inline void handle_colocations(TurnPoints& turn_points)
+template
+<
+    bool Reverse1, bool Reverse2,
+    typename Turns,
+    typename Clusters,
+    typename Geometry1,
+    typename Geometry2
+>
+inline bool handle_colocations(Turns& turns, Clusters& clusters,
+        Geometry1 const& geometry1, Geometry2 const& geometry2)
 {
     typedef std::map
         <
@@ -195,9 +381,9 @@ inline void handle_colocations(TurnPoints& turn_points)
     map_type map;
 
     int index = 0;
-    for (typename boost::range_iterator<TurnPoints>::type
-            it = boost::begin(turn_points);
-         it != boost::end(turn_points);
+    for (typename boost::range_iterator<Turns>::type
+            it = boost::begin(turns);
+         it != boost::end(turns);
          ++it, ++index)
     {
         map[it->operations[0].seg_id].push_back(turn_operation_index(index, 0));
@@ -220,27 +406,43 @@ inline void handle_colocations(TurnPoints& turn_points)
 
     if (! colocations)
     {
-        return;
+        return false;
     }
 
     // Sort all vectors, per same segment
-    less_by_fraction_and_type<TurnPoints> less(turn_points);
+    less_by_fraction_and_type<Turns> less(turns);
     for (typename map_type::iterator it = map.begin();
          it != map.end(); ++it)
     {
         std::sort(it->second.begin(), it->second.end(), less);
     }
 
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::segment_ratio_type segment_ratio_type;
+
+    typedef std::map
+        <
+            segment_fraction<segment_ratio_type>,
+            signed_size_type
+        > cluster_per_segment_type;
+
+    cluster_per_segment_type cluster_per_segment;
+    signed_size_type cluster_id = 0;
+
     for (typename map_type::const_iterator it = map.begin();
          it != map.end(); ++it)
     {
-        if (it->second.size() > 1)
+        if (it->second.size() > 1u)
         {
-            handle_colocation_cluster<OverlayType>(turn_points,
-                    it->first, it->second);
+            handle_colocation_cluster<Reverse1, Reverse2>(turns, cluster_id,
+                cluster_per_segment, it->second,
+                geometry1, geometry2);
         }
     }
 
+    assign_cluster_to_turns(turns, clusters, cluster_per_segment);
+    remove_clusters(turns, clusters);
+
 #if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS)
     std::cout << "*** Colocations " << map.size() << std::endl;
     for (typename map_type::const_iterator it = map.begin();
@@ -251,21 +453,123 @@ inline void handle_colocations(TurnPoints& turn_points)
              = it->second.begin(); vit != it->second.end(); ++vit)
         {
             turn_operation_index const& toi = *vit;
-            std::cout << geometry::wkt(turn_points[toi.turn_index].point)
+            std::cout << geometry::wkt(turns[toi.turn_index].point)
                 << std::boolalpha
-                << " discarded=" << turn_points[toi.turn_index].discarded
-                << " colocated=" << turn_points[toi.turn_index].colocated
-                << " " << operation_char(turn_points[toi.turn_index].operations[0].operation)
-                << " "  << turn_points[toi.turn_index].operations[0].seg_id
-                << " "  << turn_points[toi.turn_index].operations[0].fraction
-                << " // " << operation_char(turn_points[toi.turn_index].operations[1].operation)
-                << " "  << turn_points[toi.turn_index].operations[1].seg_id
-                << " "  << turn_points[toi.turn_index].operations[1].fraction
+                << " discarded=" << turns[toi.turn_index].discarded
+                << " colocated=" << turns[toi.turn_index].colocated
+                << " " << operation_char(turns[toi.turn_index].operations[0].operation)
+                << " "  << turns[toi.turn_index].operations[0].seg_id
+                << " "  << turns[toi.turn_index].operations[0].fraction
+                << " // " << operation_char(turns[toi.turn_index].operations[1].operation)
+                << " "  << turns[toi.turn_index].operations[1].seg_id
+                << " "  << turns[toi.turn_index].operations[1].fraction
                 << std::endl;
         }
     }
 #endif // DEBUG
 
+    return true;
+}
+
+
+struct is_turn_index
+{
+    is_turn_index(signed_size_type index)
+        : m_index(index)
+    {}
+
+    template <typename Indexed>
+    inline bool operator()(Indexed const& indexed) const
+    {
+        // Indexed is a indexed_turn_operation<Operation>
+        return indexed.turn_index == m_index;
+    }
+
+    std::size_t m_index;
+};
+
+
+template
+<
+    bool Reverse1, bool Reverse2,
+    typename Turns,
+    typename Clusters,
+    typename Geometry1,
+    typename Geometry2
+>
+inline void assign_startable_in_clusters(Clusters& clusters, Turns& turns,
+        operation_type for_operation,
+        Geometry1 const& geometry1, Geometry2 const& geometry2)
+{
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::point_type point_type;
+    typedef typename turn_type::turn_operation_type turn_operation_type;
+
+    // Define sorter, sorting counter-clockwise such that polygons are on the
+    // right side
+    typedef sort_by_side::side_sorter
+        <
+            Reverse1, Reverse2, point_type, std::less<int>
+        > sbs_type;
+
+    for (typename Clusters::iterator mit = clusters.begin();
+         mit != clusters.end(); ++mit)
+    {
+        std::set<signed_size_type> const& ids = mit->second;
+        if (ids.empty())
+        {
+            continue;
+        }
+
+        sbs_type sbs;
+        point_type turn_point; // should be all the same for all turns in cluster
+
+        bool first = true;
+        for (typename std::set<signed_size_type>::const_iterator sit = ids.begin();
+             sit != ids.end(); ++sit)
+        {
+            signed_size_type turn_index = *sit;
+            turn_type const& turn = turns[turn_index];
+            if (first)
+            {
+                turn_point = turn.point;
+            }
+            for (int i = 0; i < 2; i++)
+            {
+                turn_operation_type const& op = turn.operations[i];
+                sbs.add(op, turn_index, i, geometry1, geometry2, first);
+                first = false;
+            }
+        }
+        sbs.apply(turn_point);
+
+        sbs.find_open();
+
+        // Unset the startable flag for all 'closed' spaces
+        for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
+        {
+            const typename sbs_type::rp& ranked = sbs.m_ranked_points[i];
+            turn_type& turn = turns[ranked.turn_index];
+            turn_operation_type& op = turn.operations[ranked.op_index];
+
+            if (ranked.index != sort_by_side::index_to)
+            {
+                continue;
+            }
+
+            op.enriched.count_left = ranked.left_count;
+            op.enriched.count_right = ranked.right_count;
+
+            if ((for_operation == operation_union
+                    && ranked.left_count != 0)
+             || (for_operation == operation_intersection
+                    && ranked.right_count != 2))
+            {
+                op.enriched.startable = false;
+            }
+        }
+
+    }
 }
 
 
diff --git a/boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp b/boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp
deleted file mode 100644
index 277a223d47..0000000000
--- a/boost/geometry/algorithms/detail/overlay/handle_tangencies.hpp
+++ /dev/null
@@ -1,786 +0,0 @@
-// Boost.Geometry (aka GGL, Generic Geometry Library)
-
-// Copyright (c) 2007-2012 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_HANDLE_TANGENCIES_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TANGENCIES_HPP
-
-#include <algorithm>
-
-#include <boost/geometry/algorithms/detail/ring_identifier.hpp>
-#include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp>
-#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
-#include <boost/geometry/algorithms/detail/recalculate.hpp>
-
-#include <boost/geometry/policies/robustness/robust_point_type.hpp>
-#include <boost/geometry/policies/robustness/segment_ratio_type.hpp>
-#include <boost/geometry/policies/robustness/robust_type.hpp>
-
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-#include <boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp>
-#endif
-
-#include <boost/geometry/geometries/point.hpp>
-#include <boost/geometry/geometries/segment.hpp>
-
-
-// TODO: the approach below should be completely replaced by the new
-// get_left_turns, to keep the outgoing vector which has open space one of its
-// sides.
-
-
-namespace boost { namespace geometry
-{
-
-#ifndef DOXYGEN_NO_DETAIL
-namespace detail { namespace overlay
-{
-
-
-template
-<
-    typename TurnPoints,
-    typename Indexed,
-    typename Geometry1, typename Geometry2,
-    typename RobustPolicy,
-    bool Reverse1, bool Reverse2,
-    typename Strategy
->
-struct sort_in_cluster
-{
-    inline sort_in_cluster(TurnPoints const& turn_points
-            , Geometry1 const& geometry1
-            , Geometry2 const& geometry2
-            , RobustPolicy const& robust_policy
-            , Strategy const& strategy)
-        : m_turn_points(turn_points)
-        , m_geometry1(geometry1)
-        , m_geometry2(geometry2)
-        , m_rescale_policy(robust_policy)
-        , m_strategy(strategy)
-    {}
-
-private :
-
-    TurnPoints const& m_turn_points;
-    Geometry1 const& m_geometry1;
-    Geometry2 const& m_geometry2;
-    RobustPolicy const& m_rescale_policy;
-    Strategy const& m_strategy;
-
-    typedef typename Indexed::type turn_operation_type;
-    typedef typename geometry::point_type<Geometry1>::type point_type;
-
-    typedef typename geometry::robust_point_type
-    <
-        point_type,
-        RobustPolicy
-    >::type robust_point_type;
-
-    inline bool default_order(Indexed const& left, Indexed const& right) const
-    {
-        // We've nothing to sort on. Take the indexes
-        return left.turn_index < right.turn_index;
-    }
-
-    // Still necessary in some situations,
-    // for example #case_102_multi, #case_107_multi, #case_recursive_boxes_3
-    inline void get_situation_map(Indexed const& left, Indexed const& right,
-                              robust_point_type& pi_rob, robust_point_type& pj_rob,
-                              robust_point_type& ri_rob, robust_point_type& rj_rob,
-                              robust_point_type& si_rob, robust_point_type& sj_rob) const
-    {
-        point_type pi, pj, ri, rj, si, sj;
-
-        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
-            left.subject->seg_id,
-            pi, pj);
-        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
-            *left.other_seg_id,
-            ri, rj);
-        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
-            *right.other_seg_id,
-            si, sj);
-
-        geometry::recalculate(pi_rob, pi, m_rescale_policy);
-        geometry::recalculate(pj_rob, pj, m_rescale_policy);
-        geometry::recalculate(ri_rob, ri, m_rescale_policy);
-        geometry::recalculate(rj_rob, rj, m_rescale_policy);
-        geometry::recalculate(si_rob, si, m_rescale_policy);
-        geometry::recalculate(sj_rob, sj, m_rescale_policy);
-    }
-
-#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO
-    // This method was still called but did no effect whatsoever on the results,
-    // with or without robustness-rescaling.
-    // Probable cause: we rescale in this file ourselves, ignoring passed policy
-    // TODO: check this more.
-    // Besides this, it currently does not compile for yet unknown reasons
-    // (does not find specialization for segment_ratio_type)
-    // It is currently only called from the Unit Test "multi_intersection.cpp"
-
-    // Determine how p/r and p/s are located.
-    inline void overlap_info(
-        robust_point_type const& pi, robust_point_type const& pj,
-        robust_point_type const& ri, robust_point_type const& rj,
-        robust_point_type const& si, robust_point_type const& sj,
-        bool& pr_overlap, bool& ps_overlap, bool& rs_overlap) const
-    {
-        // Determine how p/r and p/s are located.
-        // One of them is coming from opposite direction.
-
-        typedef segment_intersection_points
-                <
-                    point_type,
-                    typename segment_ratio_type
-                    <
-                        point_type, RobustPolicy
-                    >::type
-                > intersection_return_type;
-
-        typedef strategy::intersection::relate_cartesian_segments
-            <
-                policies::relate::segments_intersection_points
-                    <
-                        intersection_return_type
-                    >
-            > policy;
-
-        typedef model::referring_segment<robust_point_type const> segment_type;
-        segment_type p(pi, pj);
-        segment_type r(ri, rj);
-        segment_type s(si, sj);
-
-        // Get the intersection point (or two points)
-        intersection_return_type pr = policy::apply(p, r, m_rescale_policy, pi, pj, ri, rj);
-        intersection_return_type ps = policy::apply(p, s, m_rescale_policy, pi, pj, si, sj);
-        intersection_return_type rs = policy::apply(r, s, m_rescale_policy, ri, rj, si, sj);
-
-        // Check on overlap
-        pr_overlap = pr.count == 2;
-        ps_overlap = ps.count == 2;
-        rs_overlap = rs.count == 2;
-    }
-#endif
-
-
-#ifdef BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES
-    inline void debug_consider(int order, Indexed const& left,
-            Indexed const& right, std::string const& header,
-            bool skip = true,
-            std::string const& extra = "", bool ret = false
-        ) const
-    {
-        if (skip) return;
-
-        std::cout << "Case: " << header << " for " << left.turn_index << " / " << right.turn_index << std::endl;
-
-        robust_point_type pi, pj, ri, rj, si, sj;
-        get_situation_map(left, right, pi, pj, ri, rj, si, sj);
-
-#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO
-        bool prc = false, psc = false, rsc = false;
-        overlap_info(pi, pj, ri, rj, si, sj, prc, psc, rsc);
-#endif
-
-        int const side_ri_p = m_strategy.apply(pi, pj, ri);
-        int const side_rj_p = m_strategy.apply(pi, pj, rj);
-        int const side_si_p = m_strategy.apply(pi, pj, si);
-        int const side_sj_p = m_strategy.apply(pi, pj, sj);
-        int const side_si_r = m_strategy.apply(ri, rj, si);
-        int const side_sj_r = m_strategy.apply(ri, rj, sj);
-
-#ifdef BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES_MORE
-        std::cout << " Segment p:" << geometry::wkt(pi) << " .. " << geometry::wkt(pj) << std::endl;
-        std::cout << " Segment r:" << geometry::wkt(ri) << " .. " << geometry::wkt(rj) << std::endl;
-        std::cout << " Segment s:" << geometry::wkt(si) << " .. " << geometry::wkt(sj) << std::endl;
-
-        std::cout << " r//p: " << side_ri_p << " / " << side_rj_p << std::endl;
-        std::cout << " s//p: " << side_si_p << " / " << side_sj_p << std::endl;
-        std::cout << " s//r: " << side_si_r << " / " << side_sj_r << std::endl;
-#endif
-
-        std::cout << header
-                //<< " order: " << order
-                << " ops: " << operation_char(left.subject->operation)
-                    << "/" << operation_char(right.subject->operation)
-                << " ri//p: " << side_ri_p
-                << " si//p: " << side_si_p
-                << " si//r: " << side_si_r
-#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO
-                << " cnts: " << int(prc) << ","  << int(psc) << "," << int(rsc)
-#endif
-                //<< " idx: " << left.turn_index << "/" << right.turn_index
-                ;
-
-        if (! extra.empty())
-        {
-            std::cout << " " << extra << " " << (ret ? "true" : "false");
-        }
-        std::cout << std::endl;
-    }
-#else
-    inline void debug_consider(int, Indexed const& ,
-            Indexed const& , std::string const& ,
-            bool = true,
-            std::string const& = "", bool = false
-        ) const
-    {}
-#endif
-
-
-    // ux/ux
-    inline bool consider_ux_ux(Indexed const& left,
-            Indexed const& right
-            , std::string const& // header
-        ) const
-    {
-        bool ret = left.turn_index < right.turn_index;
-
-        // In combination of u/x, x/u: take first union, then blocked.
-        // Solves #88, #61, #56, #80
-        if (left.subject->operation == operation_union
-            && right.subject->operation == operation_blocked)
-        {
-            ret = true;
-        }
-        else if (left.subject->operation == operation_blocked
-            && right.subject->operation == operation_union)
-        {
-            ret = false;
-        }
-        else
-        {
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-            std::cout << "ux/ux unhandled" << std::endl;
-#endif
-        }
-
-        //debug_consider(0, left, right, header, false, "-> return ", ret);
-
-        return ret;
-    }
-
-    inline bool consider_iu_ux(Indexed const& left,
-            Indexed const& right,
-            int order // 1: iu first, -1: ux first
-            , std::string const& // header
-        ) const
-    {
-        bool ret = false;
-
-        if (left.subject->operation == operation_union
-            && right.subject->operation == operation_union)
-        {
-            ret = order == 1;
-        }
-        else if (left.subject->operation == operation_union
-            && right.subject->operation == operation_blocked)
-        {
-            ret = true;
-        }
-        else if (right.subject->operation == operation_union
-            && left.subject->operation == operation_blocked)
-        {
-            ret = false;
-        }
-        else if (left.subject->operation == operation_union)
-        {
-            ret = true;
-        }
-        else if (right.subject->operation == operation_union)
-        {
-            ret = false;
-        }
-        else
-        {
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-            // this still happens in the traverse.cpp test
-            std::cout << " iu/ux unhandled" << std::endl;
-#endif
-            ret = order == 1;
-        }
-
-        //debug_consider(0, left, right, header, false, "-> return", ret);
-        return ret;
-    }
-
-    inline bool consider_iu_ix(Indexed const& left,
-            Indexed const& right,
-            int order // 1: iu first, -1: ix first
-            , std::string const& // header
-        ) const
-    {
-        //debug_consider(order, left, right, header, false, "iu/ix");
-
-        return left.subject->operation == operation_intersection
-                && right.subject->operation == operation_intersection ? order == 1
-            : left.subject->operation == operation_intersection ? false
-            : right.subject->operation == operation_intersection ? true
-            : order == 1;
-    }
-
-    inline bool consider_ix_ix(Indexed const& left, Indexed const& right
-            , std::string const& // header
-            ) const
-    {
-        // Take first intersection, then blocked.
-        if (left.subject->operation == operation_intersection
-            && right.subject->operation == operation_blocked)
-        {
-            return true;
-        }
-        else if (left.subject->operation == operation_blocked
-            && right.subject->operation == operation_intersection)
-        {
-            return false;
-        }
-
-        // Default case, should not occur
-
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-        std::cout << "ix/ix unhandled" << std::endl;
-#endif
-        //debug_consider(0, left, right, header, false, "-> return", ret);
-
-        return default_order(left, right);
-    }
-
-
-    inline bool consider_iu_iu(Indexed const& left, Indexed const& right,
-                    std::string const& header, bool redo = false) const
-    {
-        //debug_consider(0, left, right, header);
-
-        // In general, order it like "union, intersection".
-        if (left.subject->operation == operation_intersection
-            && right.subject->operation == operation_union)
-        {
-            //debug_consider(0, left, right, header, false, "i,u", false);
-            return false;
-        }
-        else if (left.subject->operation == operation_union
-            && right.subject->operation == operation_intersection)
-        {
-            //debug_consider(0, left, right, header, false, "u,i", true);
-            return true;
-        }
-
-        robust_point_type pi, pj, ri, rj, si, sj;
-        get_situation_map(left, right, pi, pj, ri, rj, si, sj);
-
-        int const side_ri_p = m_strategy.apply(pi, pj, ri);
-        int const side_si_p = m_strategy.apply(pi, pj, si);
-        int const side_si_r = m_strategy.apply(ri, rj, si);
-
-        // Both located at same side (#58, pie_21_7_21_0_3)
-        if (side_ri_p * side_si_p == 1 && side_si_r != 0)
-        {
-            if (left.subject->operation == operation_union
-                && right.subject->operation == operation_union)
-            {
-                int const side_ri_s = m_strategy.apply(si, sj, ri);
-                if (side_si_r == side_ri_s)
-                {
-                    return default_order(left, right);
-                }
-
-                // Take the most left one
-                bool const ret = side_si_r == 1;
-                //debug_consider(0, left, right, header, false, "same side", ret);
-                return ret;
-            }
-        }
-
-
-        // Coming from opposite sides (#59, #99)
-        if (side_ri_p * side_si_p == -1)
-        {
-            bool ret = false;
-
-            {
-                ret = side_ri_p == 1; // #100
-                debug_consider(0, left, right, header, false, "opp.", ret);
-                return ret;
-            }
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-            std::cout << " iu/iu coming from opposite unhandled" << std::endl;
-#endif
-        }
-
-#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO
-        // We need EXTRA information here: are p/r/s overlapping?
-        bool pr_ov = false, ps_ov = false, rs_ov = false;
-        overlap_info(pi, pj, ri, rj, si, sj, pr_ov, ps_ov, rs_ov);
-#else
-        // std::cout << "Boost.Geometry: skipping overlap_info" << std::endl;
-#endif
-
-        // One coming from right (#83,#90)
-        // One coming from left (#90, #94, #95)
-        if (side_si_r != 0 && (side_ri_p != 0 || side_si_p != 0))
-        {
-            int const side_ri_s = m_strategy.apply(si, sj, ri);
-            if (side_si_r == side_ri_s)
-            {
-                return default_order(left, right);
-            }
-
-            bool ret = false;
-
-#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO
-            if (pr_ov || ps_ov)
-            {
-                int r = side_ri_p != 0 ? side_ri_p : side_si_p;
-                ret = r * side_si_r == 1;
-            }
-            else
-#endif
-            {
-                ret = side_si_r == 1;
-            }
-
-            debug_consider(0, left, right, header, false, "left or right", ret);
-            return ret;
-        }
-
-        // All aligned (#92, #96)
-        if (side_ri_p == 0 && side_si_p == 0 && side_si_r == 0)
-        {
-            // One of them is coming from opposite direction.
-
-            // Take the one NOT overlapping
-            bool ret = false;
-            bool found = false;
-#if BOOST_GEOMETRY_HANDLE_TANGENCIES_WITH_OVERLAP_INFO
-            if (pr_ov && ! ps_ov)
-            {
-                ret = true;
-                found = true;
-            }
-            else if (!pr_ov && ps_ov)
-            {
-                ret = false;
-                found = true;
-            }
-#endif
-
-            debug_consider(0, left, right, header, false, "aligned", ret);
-            if (found)
-            {
-                return ret;
-            }
-        }
-
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-        std::cout << " iu/iu unhandled" << std::endl;
-        debug_consider(0, left, right, header, false, "unhandled", left.turn_index < right.turn_index);
-#endif
-        if (! redo)
-        {
-            // In some cases behaviour is not symmetrical. TODO: fix this properly
-            //   OR: alternatively we might consider calling all these functions one-way anyway
-            return ! consider_iu_iu(right, left, header, true);
-        }
-
-        return default_order(left, right);
-    }
-
-    inline bool consider_ii(Indexed const& left, Indexed const& right,
-                    std::string const& header) const
-    {
-        debug_consider(0, left, right, header);
-
-        robust_point_type pi, pj, ri, rj, si, sj;
-        get_situation_map(left, right, pi, pj, ri, rj, si, sj);
-
-        int const side_ri_p = m_strategy.apply(pi, pj, ri);
-        int const side_si_p = m_strategy.apply(pi, pj, si);
-
-        // Two other points are (mostly) lying both right of the considered segment
-        // Take the most left one
-        int const side_si_r = m_strategy.apply(ri, rj, si);
-        if (side_ri_p == -1
-            && side_si_p == -1
-            && side_si_r != 0)
-        {
-            bool const ret = side_si_r != 1;
-            return ret;
-        }
-        return default_order(left, right);
-    }
-
-
-public :
-    inline bool operator()(Indexed const& left, Indexed const& right) const
-    {
-        if ((m_turn_points[left.turn_index].discarded || left.discarded)
-            && (m_turn_points[right.turn_index].discarded || right.discarded))
-        {
-            return default_order(left, right);
-        }
-        else if (m_turn_points[left.turn_index].discarded || left.discarded)
-        {
-            // Be careful to sort discarded first, then all others
-            return true;
-        }
-        else if (m_turn_points[right.turn_index].discarded || right.discarded)
-        {
-            // See above so return false here such that right (discarded)
-            // is sorted before left (not discarded)
-            return false;
-        }
-        else if (m_turn_points[left.turn_index].combination(operation_blocked, operation_union)
-                && m_turn_points[right.turn_index].combination(operation_blocked, operation_union))
-        {
-            // ux/ux
-            return consider_ux_ux(left, right, "ux/ux");
-        }
-        else if (m_turn_points[left.turn_index].both(operation_union)
-            && m_turn_points[right.turn_index].both(operation_union))
-        {
-            // uu/uu, Order is arbitrary
-            // Note: uu/uu is discarded now before so this point will
-            //       not be reached.
-            return default_order(left, right);
-        }
-        else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_union)
-                && m_turn_points[right.turn_index].combination(operation_intersection, operation_union))
-        {
-            return consider_iu_iu(left, right, "iu/iu");
-        }
-        else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_blocked)
-                && m_turn_points[right.turn_index].combination(operation_intersection, operation_blocked))
-        {
-            return consider_ix_ix(left, right, "ix/ix");
-        }
-        else if (m_turn_points[left.turn_index].both(operation_intersection)
-                && m_turn_points[right.turn_index].both(operation_intersection))
-        {
-            return consider_ii(left, right, "ii/ii");
-        }
-        else  if (m_turn_points[left.turn_index].combination(operation_union, operation_blocked)
-                && m_turn_points[right.turn_index].combination(operation_intersection, operation_union))
-        {
-            return consider_iu_ux(left, right, -1, "ux/iu");
-        }
-        else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_union)
-                && m_turn_points[right.turn_index].combination(operation_union, operation_blocked))
-        {
-            return consider_iu_ux(left, right, 1, "iu/ux");
-        }
-        else  if (m_turn_points[left.turn_index].combination(operation_intersection, operation_blocked)
-                && m_turn_points[right.turn_index].combination(operation_intersection, operation_union))
-        {
-            return consider_iu_ix(left, right, 1, "ix/iu");
-        }
-        else if (m_turn_points[left.turn_index].combination(operation_intersection, operation_union)
-                && m_turn_points[right.turn_index].combination(operation_intersection, operation_blocked))
-        {
-            return consider_iu_ix(left, right, -1, "iu/ix");
-        }
-        else if (m_turn_points[left.turn_index].method != method_equal
-            && m_turn_points[right.turn_index].method == method_equal
-            )
-        {
-            // If one of them was EQUAL or CONTINUES, it should always come first
-            return false;
-        }
-        else if (m_turn_points[left.turn_index].method == method_equal
-            && m_turn_points[right.turn_index].method != method_equal
-            )
-        {
-            return true;
-        }
-
-        // Now we have no clue how to sort.
-
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-        std::cout << " Consider: " << operation_char(m_turn_points[left.turn_index].operations[0].operation)
-                << operation_char(m_turn_points[left.turn_index].operations[1].operation)
-                << "/" << operation_char(m_turn_points[right.turn_index].operations[0].operation)
-                << operation_char(m_turn_points[right.turn_index].operations[1].operation)
-                << " " << " Take " << left.turn_index << " < " << right.turn_index
-                << std::endl;
-#endif
-
-        return default_order(left, right);
-    }
-};
-
-
-
-template
-<
-    typename IndexType,
-    typename Iterator,
-    typename TurnPoints,
-    typename Geometry1,
-    typename Geometry2,
-    typename Strategy
->
-inline void inspect_cluster(Iterator begin_cluster, Iterator end_cluster,
-            TurnPoints& turn_points,
-            operation_type ,
-            Geometry1 const& , Geometry2 const& ,
-            Strategy const& )
-{
-    int count = 0;
-
-    // Make an analysis about all occuring cases here.
-    std::map<std::pair<operation_type, operation_type>, int> inspection;
-    for (Iterator it = begin_cluster; it != end_cluster; ++it)
-    {
-        operation_type first = turn_points[it->turn_index].operations[0].operation;
-        operation_type second = turn_points[it->turn_index].operations[1].operation;
-        if (first > second)
-        {
-            std::swap(first, second);
-        }
-        inspection[std::make_pair(first, second)]++;
-        count++;
-    }
-
-
-    bool keep_cc = false;
-
-    // Decide about which is going to be discarded here.
-    if (inspection[std::make_pair(operation_union, operation_union)] == 1
-        && inspection[std::make_pair(operation_continue, operation_continue)] == 1)
-    {
-        // In case of uu/cc, discard the uu, that indicates a tangency and
-        // inclusion would disturb the (e.g.) cc-cc-cc ordering
-        // NOTE: uu is now discarded anyhow.
-        keep_cc = true;
-    }
-    else if (count == 2
-        && inspection[std::make_pair(operation_intersection, operation_intersection)] == 1
-        && inspection[std::make_pair(operation_union, operation_intersection)] == 1)
-    {
-        // In case of ii/iu, discard the iu. The ii should always be visited,
-        // Because (in case of not discarding iu) correctly ordering of ii/iu appears impossible
-        for (Iterator it = begin_cluster; it != end_cluster; ++it)
-        {
-            if (turn_points[it->turn_index].combination(operation_intersection, operation_union))
-            {
-                it->discarded = true;
-            }
-        }
-    }
-
-    // Discard any continue turn, unless it is the only thing left
-    //    (necessary to avoid cc-only rings, all being discarded
-    //    e.g. traversal case #75)
-    int nd_count= 0, cc_count = 0;
-    for (Iterator it = begin_cluster; it != end_cluster; ++it)
-    {
-        if (! it->discarded)
-        {
-            nd_count++;
-            if (turn_points[it->turn_index].both(operation_continue))
-            {
-                cc_count++;
-            }
-        }
-    }
-
-    if (nd_count == cc_count)
-    {
-        keep_cc = true;
-    }
-
-    if (! keep_cc)
-    {
-        for (Iterator it = begin_cluster; it != end_cluster; ++it)
-        {
-            if (turn_points[it->turn_index].both(operation_continue))
-            {
-                it->discarded = true;
-            }
-        }
-    }
-}
-
-
-template
-<
-    typename IndexType,
-    bool Reverse1, bool Reverse2,
-    typename Iterator,
-    typename TurnPoints,
-    typename Geometry1,
-    typename Geometry2,
-    typename RobustPolicy,
-    typename Strategy
->
-inline void handle_cluster(Iterator begin_cluster, Iterator end_cluster,
-            TurnPoints& turn_points,
-            operation_type for_operation,
-            Geometry1 const& geometry1, Geometry2 const& geometry2,
-            RobustPolicy& robust_policy,
-            Strategy const& strategy)
-{
-    // First inspect and (possibly) discard rows
-    inspect_cluster<IndexType>(begin_cluster, end_cluster, turn_points,
-            for_operation, geometry1, geometry2, strategy);
-
-
-    // Then sort this range (discarded rows will be ordered first and will be removed in enrich_assign)
-    std::sort(begin_cluster, end_cluster,
-                sort_in_cluster
-                    <
-                        TurnPoints,
-                        IndexType,
-                        Geometry1, Geometry2,
-                        RobustPolicy,
-                        Reverse1, Reverse2,
-                        Strategy
-                    >(turn_points, geometry1, geometry2, robust_policy, strategy));
-
-#if defined(BOOST_GEOMETRY_DEBUG_HANDLE_TANGENCIES)
-    typedef typename IndexType::type operations_type;
-    operations_type const& op = turn_points[begin_cluster->turn_index].operations[begin_cluster->operation_index];
-    std::cout << std::endl << "Clustered points on equal distance " << op.fraction << std::endl;
-
-    std::cout << "->Indexes ";
-    for (Iterator it = begin_cluster; it != end_cluster; ++it)
-    {
-        std::cout << " " << it->turn_index;
-    }
-    std::cout << std::endl << "->Methods: ";
-    for (Iterator it = begin_cluster; it != end_cluster; ++it)
-    {
-        std::cout << " " << method_char(turn_points[it->turn_index].method);
-    }
-    std::cout << std::endl << "->Operations: ";
-    for (Iterator it = begin_cluster; it != end_cluster; ++it)
-    {
-        std::cout << " " << operation_char(turn_points[it->turn_index].operations[0].operation)
-            << operation_char(turn_points[it->turn_index].operations[1].operation);
-    }
-    std::cout << std::endl << "->Discarded: ";
-    for (Iterator it = begin_cluster; it != end_cluster; ++it)
-    {
-        std::cout << " " << (it->discarded ? "true" : "false");
-    }
-    std::cout << std::endl;
-        //<< "\tOn segments: "    << prev_op.seg_id  << " / "  << prev_op.other_id
-        //<< " and "  << op.seg_id << " / " << op.other_id
-        //<< geometry::distance(turn_points[prev->turn_index].point, turn_points[it->turn_index].point)
-#endif
-
-}
-
-
-}} // namespace detail::overlay
-#endif //DOXYGEN_NO_DETAIL
-
-
-}} // namespace boost::geometry
-
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TANGENCIES_HPP
diff --git a/boost/geometry/algorithms/detail/overlay/handle_touch.hpp b/boost/geometry/algorithms/detail/overlay/handle_touch.hpp
new file mode 100644
index 0000000000..1febefc83a
--- /dev/null
+++ b/boost/geometry/algorithms/detail/overlay/handle_touch.hpp
@@ -0,0 +1,336 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2015 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_HANDLE_TOUCH_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TOUCH_HPP
+
+#include <cstddef>
+
+#include <boost/range.hpp>
+
+#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
+#include <boost/geometry/geometries/concepts/check.hpp>
+#include <boost/geometry/algorithms/detail/ring_identifier.hpp>
+#include <boost/geometry/algorithms/detail/overlay/segment_identifier.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail { namespace overlay
+{
+
+
+template <typename Turns, typename Visitor>
+class handle_touch_uu
+{
+private :
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename boost::range_iterator<Turns>::type turn_iterator;
+    typedef typename boost::range_iterator<Turns const>::type turn_const_iterator;
+
+    typedef typename boost::range_iterator
+        <
+            typename turn_type::container_type const
+        >::type operation_const_iterator;
+
+public :
+
+    handle_touch_uu(Visitor& visitor)
+        : m_visitor(visitor)
+    {}
+
+    inline void apply(detail::overlay::operation_type operation, Turns& turns)
+    {
+        if (! has_uu(turns))
+        {
+            // Performance - if there is no u/u at all, nothing to be done
+            return;
+        }
+
+        // Iterate through all u/u points
+        int turn_index = 0;
+        for (turn_iterator it = boost::begin(turns);
+             it != boost::end(turns);
+             ++it, ++turn_index)
+        {
+            turn_type& turn = *it;
+            if (! turn.both(operation_union))
+            {
+                continue;
+            }
+
+            m_visitor.print("handle_touch uu:", turns, turn_index);
+
+            bool const traverse = turn_should_be_traversed(turns, turn, turn_index);
+            bool const start = traverse
+                               && turn_should_be_startable(turns, turn, turn_index);
+            m_visitor.print("handle_touch, ready ", turns, turn_index);
+//                          << std::boolalpha
+//                          << traverse << " " << start
+
+            if (traverse)
+            {
+                // Indicate the sources should switch here to create
+                // separate rings (outer ring / inner ring)
+                turn.switch_source = true;
+            }
+            // TODO: this is often not correct, fix this
+            turn.operations[0].enriched.startable = start;
+            turn.operations[1].enriched.startable = start;
+        }
+    }
+
+private :
+
+    // Generic utility to be moved somewhere else
+    static inline
+    ring_identifier ring_id_from_seg_id(const segment_identifier& seg_id)
+    {
+        return ring_identifier(seg_id.source_index,
+                               seg_id.multi_index,
+                               seg_id.ring_index);
+    }
+
+    static inline
+    ring_identifier ring_id_from_op(const turn_type& turn,
+                                    int operation_index)
+    {
+        return ring_id_from_seg_id(turn.operations[operation_index].seg_id);
+    }
+
+    static inline bool in_range(const Turns& turns, signed_size_type index)
+    {
+        signed_size_type const turns_size =
+                static_cast<signed_size_type>(boost::size(turns));
+        return index >= 0 && index < turns_size;
+    }
+
+    static inline bool has_uu(const Turns& turns)
+    {
+        for (turn_const_iterator it = boost::begin(turns);
+            it != boost::end(turns);
+            ++it)
+        {
+            const turn_type& turn = *it;
+            if (turn.both(operation_union))
+            {
+                return true;
+            }
+        }
+        return false;
+    }
+
+    static inline
+    bool turn_should_be_startable(const Turns& turns,
+                                  const turn_type& uu_turn,
+                                  signed_size_type uu_turn_index)
+    {
+        return turn_startable(turns, uu_turn, 0, uu_turn_index)
+            || turn_startable(turns, uu_turn, 1, uu_turn_index);
+
+    }
+
+    static inline
+    bool turn_startable(const Turns& turns,
+                        const turn_type& uu_turn,
+                        std::size_t op_index,
+                        signed_size_type original_turn_index,
+                        std::size_t iteration = 0)
+    {
+        if (iteration >= boost::size(turns))
+        {
+            // Defensive check to avoid infinite recursion
+            return false;
+        }
+
+        signed_size_type const index
+                = uu_turn.operations[op_index].enriched.travels_to_ip_index;
+        if (index == original_turn_index)
+        {
+            // Completely traveled, having u/u only, via this op_index
+            return true;
+        }
+
+        if (! in_range(turns, index))
+        {
+            return false;
+        }
+
+        const turn_type& new_turn = turns[index];
+
+        if (new_turn.operations[0].enriched.startable)
+        {
+            // Already selectable - no need to select u/u turn too
+            return false;
+        }
+
+        // If this u/u turn is traversed normally (without skipping), sources are switched
+        return turn_startable(turns, new_turn, 1 - op_index,
+                              original_turn_index, iteration + 1);
+    }
+
+    inline bool turn_should_be_traversed(const Turns& turns,
+                                  const turn_type& uu_turn,
+                                  signed_size_type uu_turn_index)
+    {
+        return turn_should_be_traversed(turns, uu_turn, uu_turn_index, 0)
+            || turn_should_be_traversed(turns, uu_turn, uu_turn_index, 1);
+    }
+
+    inline bool turn_should_be_traversed(const Turns& turns,
+                                  const turn_type& uu_turn,
+                                  signed_size_type uu_turn_index,
+                                  int uu_operation_index)
+    {
+        // Suppose this is a u/u turn between P and Q
+        // Examine all other turns on P and check if Q can be reached
+        // Use one of the operations and check if you can reach the other
+        signed_size_type const to_turn_index
+            = uu_turn.operations[uu_operation_index].enriched.travels_to_ip_index;
+        if (! in_range(turns, to_turn_index))
+        {
+            return false;
+        }
+
+        m_visitor.print("Examine:", turns, to_turn_index);
+        ring_identifier const other_ring_id
+            = ring_id_from_op(uu_turn, 1 - uu_operation_index);
+
+        bool complete = false;
+        return can_reach(complete, turns, turns[to_turn_index], uu_operation_index,
+                         other_ring_id, uu_turn_index, to_turn_index);
+    }
+
+    inline bool can_reach(bool& complete, const Turns& turns,
+                                 const turn_type& turn,
+                                 signed_size_type uu_operation_index,
+                                 const ring_identifier& target_ring_id,
+                                 signed_size_type uu_turn_index,
+                                 signed_size_type to_turn_index,
+                                 std::size_t iteration = 0)
+    {
+        if (complete)
+        {
+            return false;
+        }
+
+        if (turn.cluster_id >= 0)
+        {
+            // Clustered turns are yet not supported
+            return false;
+        }
+
+        if (iteration >= boost::size(turns))
+        {
+            m_visitor.print("Too much iterations");
+            // Defensive check to avoid infinite recursion
+            return false;
+        }
+
+        if (uu_operation_index != -1 && turn.both(operation_union))
+        {
+            // If we end up in a u/u turn, check the way how, for this operation
+            m_visitor.print("Via u/u");
+            return can_reach_via(complete, turns, uu_operation_index,
+                                 turn.operations[uu_operation_index],
+                                 target_ring_id,
+                                 uu_turn_index, to_turn_index, iteration);
+        }
+        else
+        {
+            // Check if specified ring can be reached via one of both operations
+            return can_reach_via(complete, turns, 0, turn.operations[0], target_ring_id,
+                                 uu_turn_index, to_turn_index, iteration)
+                || can_reach_via(complete, turns, 1, turn.operations[1], target_ring_id,
+                                 uu_turn_index, to_turn_index, iteration);
+        }
+    }
+
+    template <typename Operation>
+    inline bool can_reach_via(bool& complete, const Turns& turns,
+            signed_size_type operation_index,
+            const Operation& operation,
+            const ring_identifier& target_ring_id,
+            signed_size_type uu_turn_index,
+            signed_size_type to_turn_index,
+            std::size_t iteration = 0)
+    {
+        if (operation.operation != operation_union
+            && operation.operation != operation_continue)
+        {
+            return false;
+        }
+
+        signed_size_type const index = operation.enriched.travels_to_ip_index;
+        if (index == to_turn_index)
+        {
+            m_visitor.print("Dead end at", turns, index);
+            // Completely traveled, the target is not found
+            return false;
+        }
+        if (index == uu_turn_index)
+        {
+            // End up where trial was started
+            m_visitor.print("Travel complete at", turns, index);
+            complete = true;
+            return false;
+        }
+
+        if (! in_range(turns, index))
+        {
+            return false;
+        }
+
+        m_visitor.print("Now to", turns, index, operation_index);
+        const turn_type& new_turn = turns[index];
+
+        if (new_turn.both(operation_union))
+        {
+            ring_identifier const ring_id = ring_id_from_op(new_turn, operation_index);
+            if (ring_id == target_ring_id)
+            {
+                m_visitor.print("Found (at u/u)!");
+                return true;
+            }
+        }
+        else
+        {
+            ring_identifier const ring_id1 = ring_id_from_op(new_turn, 0);
+            ring_identifier const ring_id2 = ring_id_from_op(new_turn, 1);
+            if (ring_id1 == target_ring_id || ring_id2 == target_ring_id)
+            {
+                m_visitor.print("Found!");
+                return true;
+            }
+        }
+
+        // Recursively check this turn
+        return can_reach(complete, turns, new_turn, operation_index, target_ring_id,
+                         uu_turn_index, to_turn_index, iteration + 1);
+    }
+
+private :
+    Visitor m_visitor;
+};
+
+template <typename Turns, typename Visitor>
+inline void handle_touch(detail::overlay::operation_type operation,
+                         Turns& turns, Visitor& visitor)
+{
+    handle_touch_uu<Turns, Visitor> handler(visitor);
+    handler.apply(operation, turns);
+}
+
+}} // namespace detail::overlay
+#endif // DOXYGEN_NO_DETAIL
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_HANDLE_TOUCH_HPP
diff --git a/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp b/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp
new file mode 100644
index 0000000000..21868a2939
--- /dev/null
+++ b/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp
@@ -0,0 +1,203 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2007-2015 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_SORT_ON_SEGMENT_RATIO_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_ON_SEGMENT_RATIO_HPP
+
+#include <cstddef>
+#include <algorithm>
+#include <map>
+#include <set>
+#include <vector>
+
+#include <boost/range.hpp>
+
+#include <boost/geometry/algorithms/detail/overlay/copy_segment_point.hpp>
+#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
+#include <boost/geometry/strategies/side.hpp>
+
+namespace boost { namespace geometry
+{
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail { namespace overlay
+{
+
+// Wraps "turn_operation" from turn_info.hpp,
+// giving it extra information, necessary for sorting
+template <typename TurnOperation>
+struct indexed_turn_operation
+{
+    typedef TurnOperation type;
+
+    std::size_t turn_index;
+    std::size_t operation_index;
+    bool skip;
+    // use pointers to avoid copies, const& is not possible because of usage in vector
+    segment_identifier const* other_seg_id; // segment id of other segment of intersection of two segments
+    TurnOperation const* subject;
+
+    inline indexed_turn_operation(std::size_t ti, std::size_t oi,
+                TurnOperation const& sub,
+                segment_identifier const& oid)
+        : turn_index(ti)
+        , operation_index(oi)
+        , skip(false)
+        , other_seg_id(&oid)
+        , subject(boost::addressof(sub))
+    {}
+
+};
+
+template
+<
+    typename Turns,
+    typename Indexed,
+    typename Geometry1, typename Geometry2,
+    typename RobustPolicy,
+    bool Reverse1, bool Reverse2
+>
+struct less_by_segment_ratio
+{
+    inline less_by_segment_ratio(Turns const& turns
+            , operation_type for_operation
+            , Geometry1 const& geometry1
+            , Geometry2 const& geometry2
+            , RobustPolicy const& robust_policy)
+        : m_turns(turns)
+        , m_for_operation(for_operation)
+        , m_geometry1(geometry1)
+        , m_geometry2(geometry2)
+        , m_robust_policy(robust_policy)
+    {
+    }
+
+private :
+
+    Turns const& m_turns;
+    operation_type m_for_operation;
+    Geometry1 const& m_geometry1;
+    Geometry2 const& m_geometry2;
+    RobustPolicy const& m_robust_policy;
+
+    typedef typename geometry::point_type<Geometry1>::type point_type;
+
+    inline bool default_order(Indexed const& left, Indexed const& right) const
+    {
+        // We've nothing to sort on. Take the indexes
+        return left.turn_index < right.turn_index;
+    }
+
+    inline bool consider_relative_order(Indexed const& left,
+                    Indexed const& right) const
+    {
+        point_type pi, pj, ri, rj, si, sj;
+
+        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
+            left.subject->seg_id,
+            pi, pj);
+        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
+            *left.other_seg_id,
+            ri, rj);
+        geometry::copy_segment_points<Reverse1, Reverse2>(m_geometry1, m_geometry2,
+            *right.other_seg_id,
+            si, sj);
+
+        typedef typename strategy::side::services::default_strategy
+            <
+                typename cs_tag<point_type>::type
+            >::type strategy;
+
+        int const side_rj_p = strategy::apply(pi, pj, rj);
+        int const side_sj_p = strategy::apply(pi, pj, sj);
+
+        // Put the one turning left (1; right == -1) as last
+        if (side_rj_p != side_sj_p)
+        {
+            return side_rj_p < side_sj_p;
+        }
+
+        int const side_sj_r = strategy::apply(ri, rj, sj);
+        int const side_rj_s = strategy::apply(si, sj, rj);
+
+        // If they both turn left: the most left as last
+        // If they both turn right: this is not relevant, but take also here most left
+        if (side_rj_s != side_sj_r)
+        {
+            return side_rj_s < side_sj_r;
+        }
+
+        return default_order(left, right);
+    }
+
+
+public :
+
+    // Note that left/right do NOT correspond to m_geometry1/m_geometry2
+    // but to the "indexed_turn_operation"
+    inline bool operator()(Indexed const& left, Indexed const& right) const
+    {
+        if (! (left.subject->seg_id == right.subject->seg_id))
+        {
+            return left.subject->seg_id < right.subject->seg_id;
+        }
+
+        // Both left and right are located on the SAME segment.
+
+        if (! (left.subject->fraction == right.subject->fraction))
+        {
+            return left.subject->fraction < right.subject->fraction;
+        }
+
+
+        typedef typename boost::range_value<Turns>::type turn_type;
+        turn_type const& left_turn = m_turns[left.turn_index];
+        turn_type const& right_turn = m_turns[right.turn_index];
+
+        // First check "real" intersection (crosses)
+        // -> distance zero due to precision, solve it by sorting
+        if (left_turn.method == method_crosses
+            && right_turn.method == method_crosses)
+        {
+            return consider_relative_order(left, right);
+        }
+
+        bool const left_both_xx = left_turn.both(operation_blocked);
+        bool const right_both_xx = right_turn.both(operation_blocked);
+        if (left_both_xx && ! right_both_xx)
+        {
+            return true;
+        }
+        if (! left_both_xx && right_both_xx)
+        {
+            return false;
+        }
+
+        bool const left_both_uu = left_turn.both(operation_union);
+        bool const right_both_uu = right_turn.both(operation_union);
+        if (left_both_uu && ! right_both_uu)
+        {
+            return true;
+        }
+        if (! left_both_uu && right_both_uu)
+        {
+            return false;
+        }
+
+        return default_order(left, right);
+    }
+};
+
+
+}} // namespace detail::overlay
+#endif //DOXYGEN_NO_DETAIL
+
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_ON_SEGMENT_RATIO_HPP
diff --git a/boost/geometry/algorithms/detail/overlay/overlay.hpp b/boost/geometry/algorithms/detail/overlay/overlay.hpp
index 6eb0b8864c..c3ecaa0b01 100644
--- a/boost/geometry/algorithms/detail/overlay/overlay.hpp
+++ b/boost/geometry/algorithms/detail/overlay/overlay.hpp
@@ -26,6 +26,7 @@
 #include <boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp>
 #include <boost/geometry/algorithms/detail/overlay/enrichment_info.hpp>
 #include <boost/geometry/algorithms/detail/overlay/get_turns.hpp>
+#include <boost/geometry/algorithms/detail/overlay/handle_touch.hpp>
 #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp>
 #include <boost/geometry/algorithms/detail/overlay/traverse.hpp>
 #include <boost/geometry/algorithms/detail/overlay/traversal_info.hpp>
@@ -59,28 +60,50 @@ namespace detail { namespace overlay
 {
 
 
-template <typename TurnPoints, typename TurnInfoMap>
-inline void get_ring_turn_info(TurnInfoMap& turn_info_map,
-        TurnPoints const& turn_points)
+//! Default visitor for overlay, doing nothing
+struct overlay_null_visitor
 {
-    typedef typename boost::range_value<TurnPoints>::type turn_point_type;
-    typedef typename turn_point_type::container_type container_type;
+    void print(char const* ) {}
 
-    for (typename boost::range_iterator<TurnPoints const>::type
-            it = boost::begin(turn_points);
-         it != boost::end(turn_points);
+    template <typename Turns>
+    void print(char const* , Turns const& , int) {}
+
+    template <typename Turns>
+    void print(char const* , Turns const& , int , int ) {}
+
+    template <typename Turns>
+    void visit_turns(int , Turns const& ) {}
+
+    template <typename Clusters, typename Turns>
+    void visit_clusters(Clusters const& , Turns const& ) {}
+
+    template <typename Turns, typename Turn, typename Operation>
+    void visit_traverse(Turns const& , Turn const& , Operation const& , char const*)
+    {}
+
+    template <typename Turns, typename Turn, typename Operation>
+    void visit_traverse_reject(Turns const& , Turn const& , Operation const& , traverse_error_type )
+    {}
+};
+
+template <typename Turns, typename TurnInfoMap>
+inline void get_ring_turn_info(TurnInfoMap& turn_info_map, Turns const& turns)
+{
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::container_type container_type;
+
+    for (typename boost::range_iterator<Turns const>::type
+            it = boost::begin(turns);
+         it != boost::end(turns);
          ++it)
     {
-        typename boost::range_value<TurnPoints>::type const& turn_info = *it;
-        bool both_uu = turn_info.both(operation_union);
-        bool skip = (turn_info.discarded || both_uu)
-            && ! turn_info.any_blocked()
-            && ! turn_info.both(operation_intersection)
-            ;
+        typename boost::range_value<Turns>::type const& turn_info = *it;
 
-        if (! both_uu && turn_info.colocated)
+        if (turn_info.discarded
+            && ! turn_info.any_blocked()
+            && ! turn_info.colocated)
         {
-            skip = true;
+            continue;
         }
 
         for (typename boost::range_iterator<container_type const>::type
@@ -88,21 +111,13 @@ inline void get_ring_turn_info(TurnInfoMap& turn_info_map,
             op_it != boost::end(turn_info.operations);
             ++op_it)
         {
-            ring_identifier ring_id
+            ring_identifier const ring_id
                 (
                     op_it->seg_id.source_index,
                     op_it->seg_id.multi_index,
                     op_it->seg_id.ring_index
                 );
-
-            if (! skip)
-            {
-                turn_info_map[ring_id].has_normal_turn = true;
-            }
-            else if (both_uu)
-            {
-                turn_info_map[ring_id].has_uu_turn = true;
-            }
+            turn_info_map[ring_id].has_normal_turn = true;
         }
     }
 }
@@ -164,12 +179,13 @@ template
 >
 struct overlay
 {
-    template <typename RobustPolicy, typename OutputIterator, typename Strategy>
+    template <typename RobustPolicy, typename OutputIterator, typename Strategy, typename Visitor>
     static inline OutputIterator apply(
                 Geometry1 const& geometry1, Geometry2 const& geometry2,
                 RobustPolicy const& robust_policy,
                 OutputIterator out,
-                Strategy const& )
+                Strategy const& ,
+                Visitor& visitor)
     {
         bool const is_empty1 = geometry::is_empty(geometry1);
         bool const is_empty2 = geometry::is_empty(geometry2);
@@ -193,14 +209,23 @@ struct overlay
             point_type,
             typename geometry::segment_ratio_type<point_type, RobustPolicy>::type
         > turn_info;
-        typedef std::deque<turn_info> container_type;
+        typedef std::deque<turn_info> turn_container_type;
 
         typedef std::deque
             <
                 typename geometry::ring_type<GeometryOut>::type
             > ring_container_type;
 
-        container_type turn_points;
+        // Define the clusters, mapping cluster_id -> turns
+        typedef std::map
+            <
+                signed_size_type,
+                std::set<signed_size_type>
+            > cluster_type;
+
+        cluster_type clusters;
+
+        turn_container_type turns;
 
 #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
 std::cout << "get turns" << std::endl;
@@ -210,20 +235,42 @@ std::cout << "get turns" << std::endl;
             <
                 Reverse1, Reverse2,
                 detail::overlay::assign_null_policy
-            >(geometry1, geometry2, robust_policy, turn_points, policy);
+            >(geometry1, geometry2, robust_policy, turns, policy);
+
+        visitor.visit_turns(1, turns);
+
+        static const operation_type op_type
+                = OverlayType == overlay_union
+                  ? geometry::detail::overlay::operation_union
+                  : geometry::detail::overlay::operation_intersection;
 
 #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
 std::cout << "enrich" << std::endl;
 #endif
         typename Strategy::side_strategy_type side_strategy;
-        geometry::enrich_intersection_points<Reverse1, Reverse2, OverlayType>(turn_points,
-                OverlayType == overlay_union
-                    ? geometry::detail::overlay::operation_union
-                    : geometry::detail::overlay::operation_intersection,
+        geometry::enrich_intersection_points<Reverse1, Reverse2, OverlayType>(turns,
+                clusters, op_type,
                     geometry1, geometry2,
                     robust_policy,
                     side_strategy);
 
+        visitor.visit_turns(2, turns);
+
+        visitor.visit_clusters(clusters, turns);
+
+
+#if 0
+        // TODO: does not work always correctly, move to traverse and fix
+        if (op_type == geometry::detail::overlay::operation_union)
+        {
+            #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+            std::cout << "handle_touch" << std::endl;
+            #endif
+
+            handle_touch(op_type, turns, visitor);
+        }
+#endif
+
 #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
 std::cout << "traverse" << std::endl;
 #endif
@@ -231,18 +278,17 @@ std::cout << "traverse" << std::endl;
         // Note that these rings are always in clockwise order, even in CCW polygons,
         // and are marked as "to be reversed" below
         ring_container_type rings;
-        traverse<Reverse1, Reverse2, Geometry1, Geometry2>::apply
+        traverse<Reverse1, Reverse2, Geometry1, Geometry2, op_type>::apply
                 (
                     geometry1, geometry2,
-                    OverlayType == overlay_union
-                        ? geometry::detail::overlay::operation_union
-                        : geometry::detail::overlay::operation_intersection,
                     robust_policy,
-                    turn_points, rings
+                    turns, rings,
+                    clusters,
+                    visitor
                 );
 
         std::map<ring_identifier, ring_turn_info> turn_info_per_ring;
-        get_ring_turn_info(turn_info_per_ring, turn_points);
+        get_ring_turn_info(turn_info_per_ring, turns);
 
         typedef ring_properties
         <
@@ -272,6 +318,17 @@ std::cout << "traverse" << std::endl;
 
         return add_rings<GeometryOut>(selected_ring_properties, geometry1, geometry2, rings, out);
     }
+
+    template <typename RobustPolicy, typename OutputIterator, typename Strategy>
+    static inline OutputIterator apply(
+                Geometry1 const& geometry1, Geometry2 const& geometry2,
+                RobustPolicy const& robust_policy,
+                OutputIterator out,
+                Strategy const& strategy)
+    {
+        overlay_null_visitor visitor;
+        return apply(geometry1, geometry2, robust_policy, out, strategy, visitor);
+    }
 };
 
 
diff --git a/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp b/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp
index e77a163dd5..5447c8813e 100644
--- a/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp
+++ b/boost/geometry/algorithms/detail/overlay/segment_identifier.hpp
@@ -38,6 +38,7 @@ struct segment_identifier
         , multi_index(-1)
         , ring_index(-1)
         , segment_index(-1)
+        , piece_index(-1)
     {}
 
     inline segment_identifier(signed_size_type src,
@@ -48,6 +49,7 @@ struct segment_identifier
         , multi_index(mul)
         , ring_index(rin)
         , segment_index(seg)
+        , piece_index(-1)
     {}
 
     inline bool operator<(segment_identifier const& other) const
@@ -85,6 +87,9 @@ struct segment_identifier
     signed_size_type multi_index;
     signed_size_type ring_index;
     signed_size_type segment_index;
+
+    // For buffer - todo: move this to buffer-only
+    signed_size_type piece_index;
 };
 
 
diff --git a/boost/geometry/algorithms/detail/overlay/select_rings.hpp b/boost/geometry/algorithms/detail/overlay/select_rings.hpp
index d18e012b2d..1b3cd866d7 100644
--- a/boost/geometry/algorithms/detail/overlay/select_rings.hpp
+++ b/boost/geometry/algorithms/detail/overlay/select_rings.hpp
@@ -35,13 +35,11 @@ namespace detail { namespace overlay
 
 struct ring_turn_info
 {
-    bool has_uu_turn;
     bool has_normal_turn;
     bool within_other;
 
     ring_turn_info()
-        : has_uu_turn(false)
-        , has_normal_turn(false)
+        : has_normal_turn(false)
         , within_other(false)
     {}
 };
@@ -174,7 +172,7 @@ struct decide<overlay_union>
 {
     static bool include(ring_identifier const& , ring_turn_info const& info)
     {
-        return info.has_uu_turn || ! info.within_other;
+        return ! info.within_other;
     }
 
     static bool reversed(ring_identifier const& , ring_turn_info const& )
@@ -190,19 +188,16 @@ struct decide<overlay_difference>
     {
         // Difference: A - B
 
-        // If this is A (source_index=0) and there is only a u/u turn,
-        // then the ring is inside B
-        // If this is B (source_index=1) and there is only a u/u turn,
-        // then the ring is NOT inside A
+        // If this is A (source_index=0), then the ring is inside B
+        // If this is B (source_index=1), then the ring is NOT inside A
 
         // If this is A and the ring is within the other geometry,
         // then we should NOT include it.
         // If this is B then we SHOULD include it.
 
-        bool const is_first = id.source_index == 0;
-        bool const within_other = info.within_other
-            || (is_first && info.has_uu_turn);
-        return is_first ? ! within_other : within_other;
+        return id.source_index == 0
+            ? ! info.within_other
+            : info.within_other;
     }
 
     static bool reversed(ring_identifier const& id, ring_turn_info const& info)
@@ -219,7 +214,7 @@ struct decide<overlay_intersection>
 {
     static bool include(ring_identifier const& , ring_turn_info const& info)
     {
-        return ! info.has_uu_turn && info.within_other;
+        return info.within_other;
     }
 
     static bool reversed(ring_identifier const& , ring_turn_info const& )
@@ -266,19 +261,16 @@ inline void update_ring_selection(Geometry1 const& geometry1,
             continue;
         }
 
-        if (! info.has_uu_turn)
+        // Check if the ring is within the other geometry, by taking
+        // a point lying on the ring
+        switch(id.source_index)
         {
-            // Check if the ring is within the other geometry, by taking
-            // a point lying on the ring
-            switch(id.source_index)
-            {
-                case 0 :
-                    info.within_other = geometry::within(it->second.point, geometry2);
-                    break;
-                case 1 :
-                    info.within_other = geometry::within(it->second.point, geometry1);
-                    break;
-            }
+            case 0 :
+                info.within_other = geometry::within(it->second.point, geometry2);
+                break;
+            case 1 :
+                info.within_other = geometry::within(it->second.point, geometry1);
+                break;
         }
 
         if (decide<OverlayType>::include(id, info))
diff --git a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
new file mode 100644
index 0000000000..894cddab87
--- /dev/null
+++ b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
@@ -0,0 +1,434 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2015 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_SORT_BY_SIDE_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP
+
+#include <algorithm>
+#include <vector>
+
+#include <boost/geometry/algorithms/detail/direction_code.hpp>
+#include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
+#include <boost/geometry/strategies/side.hpp>
+
+namespace boost { namespace geometry
+{
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail { namespace overlay { namespace sort_by_side
+{
+
+enum index_type { index_unknown = -1, index_from = 0, index_to = 1 };
+
+// Point-wrapper, adding some properties
+template <typename Point>
+struct ranked_point
+{
+    ranked_point()
+        : main_rank(0)
+        , turn_index(-1)
+        , op_index(-1)
+        , index(index_unknown)
+        , left_count(0)
+        , right_count(0)
+        , operation(operation_none)
+    {}
+
+    ranked_point(const Point& p, signed_size_type ti, signed_size_type oi,
+                 index_type i, operation_type op, segment_identifier sid)
+        : point(p)
+        , main_rank(0)
+        , turn_index(ti)
+        , op_index(oi)
+        , index(i)
+        , left_count(0)
+        , right_count(0)
+        , operation(op)
+        , seg_id(sid)
+    {}
+
+    Point point;
+    std::size_t main_rank;
+    signed_size_type turn_index;
+    signed_size_type op_index;
+    index_type index;
+    std::size_t left_count;
+    std::size_t right_count;
+    operation_type operation;
+    segment_identifier seg_id;
+};
+
+struct less_by_turn_index
+{
+    template <typename T>
+    inline bool operator()(const T& first, const T& second) const
+    {
+        return first.turn_index == second.turn_index
+            ? first.index < second.index
+            : first.turn_index < second.turn_index
+            ;
+    }
+};
+
+struct less_by_index
+{
+    template <typename T>
+    inline bool operator()(const T& first, const T& second) const
+    {
+        // First order by from/to
+        if (first.index != second.index)
+        {
+            return first.index < second.index;
+        }
+        // All the same, order by turn index (we might consider length too)
+        return first.turn_index < second.turn_index;
+    }
+};
+
+struct less_false
+{
+    template <typename T>
+    inline bool operator()(const T&, const T& ) const
+    {
+        return false;
+    }
+};
+
+template <typename Point, typename LessOnSame, typename Compare>
+struct less_by_side
+{
+    typedef typename strategy::side::services::default_strategy
+        <
+            typename cs_tag<Point>::type
+        >::type side;
+
+    less_by_side(const Point& p1, const Point& p2)
+        : m_p1(p1)
+        , m_p2(p2)
+    {}
+
+    template <typename T>
+    inline bool operator()(const T& first, const T& second) const
+    {
+        LessOnSame on_same;
+        Compare compare;
+
+        int const side_first = side::apply(m_p1, m_p2, first.point);
+        int const side_second = side::apply(m_p1, m_p2, second.point);
+
+        if (side_first == 0 && side_second == 0)
+        {
+            // Both collinear. They might point into different directions: <------*------>
+            // If so, order the one going backwards as the very first.
+
+            int const first_code = direction_code(m_p1, m_p2, first.point);
+            int const second_code = direction_code(m_p1, m_p2, second.point);
+
+            // Order by code, backwards first, then forward.
+            return first_code != second_code
+                ? first_code < second_code
+                : on_same(first, second)
+                ;
+        }
+        else if (side_first == 0
+                && direction_code(m_p1, m_p2, first.point) == -1)
+        {
+            // First collinear and going backwards.
+            // Order as the very first, so return always true
+            return true;
+        }
+        else if (side_second == 0
+            && direction_code(m_p1, m_p2, second.point) == -1)
+        {
+            // Second is collinear and going backwards
+            // Order as very last, so return always false
+            return false;
+        }
+
+        // They are not both collinear
+
+        if (side_first != side_second)
+        {
+            return compare(side_first, side_second);
+        }
+
+        // They are both left, both right, and/or both collinear (with each other and/or with p1,p2)
+        // Check mutual side
+        int const side_second_wrt_first = side::apply(m_p2, first.point, second.point);
+
+        if (side_second_wrt_first == 0)
+        {
+            return on_same(first, second);
+        }
+
+        int const side_first_wrt_second = -side_second_wrt_first;
+
+        // Both are on same side, and not collinear
+        // Union: return true if second is right w.r.t. first, so -1,
+        // so other is 1. union has greater as compare functor
+        // Intersection: v.v.
+        return compare(side_first_wrt_second, side_second_wrt_first);
+    }
+
+private :
+    Point m_p1, m_p2;
+};
+
+template <bool Reverse1, bool Reverse2, typename Point, typename Compare>
+struct side_sorter
+{
+    typedef ranked_point<Point> rp;
+
+    inline void set_origin(Point const& origin)
+    {
+        m_origin = origin;
+    }
+
+    template <typename Operation, typename Geometry1, typename Geometry2>
+    void add(Operation const& op, signed_size_type turn_index, signed_size_type op_index,
+            Geometry1 const& geometry1,
+            Geometry2 const& geometry2,
+            bool is_origin)
+    {
+        Point point1, point2, point3;
+        geometry::copy_segment_points<Reverse1, Reverse2>(geometry1, geometry2,
+                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, index_from, op.operation, op.seg_id));
+        m_ranked_points.push_back(rp(point_to, turn_index, op_index, index_to, op.operation, op.seg_id));
+
+        if (is_origin)
+        {
+            m_origin = point1;
+        }
+    }
+
+    void apply(Point const& turn_point)
+    {
+        // We need three compare functors:
+        // 1) to order clockwise (union) or counter clockwise (intersection)
+        // 2) to order by side, resulting in unique ranks for all points
+        // 3) to order by side, resulting in non-unique ranks
+        //    to give colinear points
+
+        // Sort by side and assign rank
+        less_by_side<Point, less_by_index, Compare> less_unique(m_origin, turn_point);
+        less_by_side<Point, less_false, Compare> less_non_unique(m_origin, turn_point);
+
+        std::sort(m_ranked_points.begin(), m_ranked_points.end(), less_unique);
+
+        std::size_t colinear_rank = 0;
+        for (std::size_t i = 0; i < m_ranked_points.size(); i++)
+        {
+            if (i > 0
+                && less_non_unique(m_ranked_points[i - 1], m_ranked_points[i]))
+            {
+                // It is not collinear
+                colinear_rank++;
+            }
+
+            m_ranked_points[i].main_rank = colinear_rank;
+        }
+    }
+
+    template <signed_size_type segment_identifier::*Member, typename Map>
+    void find_open_generic(Map& handled)
+    {
+        for (std::size_t i = 0; i < m_ranked_points.size(); i++)
+        {
+            const rp& ranked = m_ranked_points[i];
+            if (ranked.index != index_from)
+            {
+                continue;
+            }
+
+            signed_size_type const& index = ranked.seg_id.*Member;
+            if (! handled[index])
+            {
+                find_polygons_for_source<Member>(index, i);
+                handled[index] = true;
+            }
+        }
+    }
+
+    void find_open()
+    {
+        // TODO: we might pass Buffer as overlay_type, instead on the fly below
+        bool one_source = true;
+        for (std::size_t i = 0; i < m_ranked_points.size(); i++)
+        {
+            const rp& ranked = m_ranked_points[i];
+            signed_size_type const& src = ranked.seg_id.source_index;
+            if (src != 0)
+            {
+                one_source = false;
+                break;
+            }
+        }
+
+        if (one_source)
+        {
+            // by multi index
+            std::map<signed_size_type, bool> handled;
+            find_open_generic
+                <
+                    &segment_identifier::piece_index
+                >(handled);
+        }
+        else
+        {
+            // by source (there should only source 0,1) TODO assert this
+            bool handled[2] = {false, false};
+            find_open_generic
+                <
+                    &segment_identifier::source_index
+                >(handled);
+        }
+    }
+
+    void reverse()
+    {
+        if (m_ranked_points.empty())
+        {
+            return;
+        }
+
+        int const last = 1 + m_ranked_points.back().main_rank;
+
+        // Move iterator after main_rank==0
+        bool has_first = false;
+        typename container_type::iterator it = m_ranked_points.begin() + 1;
+        for (; it != m_ranked_points.end() && it->main_rank == 0; ++it)
+        {
+            has_first = true;
+        }
+
+        if (has_first)
+        {
+            // Reverse first part (having main_rank == 0), if any,
+            // but skip the very first row
+            std::reverse(m_ranked_points.begin() + 1, it);
+            for (typename container_type::iterator fit = m_ranked_points.begin();
+                 fit != it; ++fit)
+            {
+                BOOST_ASSERT(fit->main_rank == 0);
+            }
+        }
+
+        // Reverse the rest (main rank > 0)
+        std::reverse(it, m_ranked_points.end());
+        for (; it != m_ranked_points.end(); ++it)
+        {
+            BOOST_ASSERT(it->main_rank > 0);
+            it->main_rank = last - it->main_rank;
+        }
+    }
+
+//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;
+        return result >= m_ranked_points.size() ? 0 : result;
+    }
+
+    //! member is pointer to member (source_index or multi_index)
+    template <signed_size_type segment_identifier::*Member>
+    std::size_t move(signed_size_type member_index, std::size_t index) const
+    {
+        std::size_t result = move(index);
+        while (m_ranked_points[result].seg_id.*Member != member_index)
+        {
+            result = move(result);
+        }
+        return result;
+    }
+
+    void assign_ranks(std::size_t min_rank, std::size_t max_rank, int side_index)
+    {
+        for (std::size_t i = 0; i < m_ranked_points.size(); i++)
+        {
+            rp& ranked = m_ranked_points[i];
+            // Suppose there are 8 ranks, if min=4,max=6: assign 4,5,6
+            // if min=5,max=2: assign from 5,6,7,1,2
+            bool const in_range
+                = max_rank >= min_rank
+                ? ranked.main_rank >= min_rank && ranked.main_rank <= max_rank
+                : ranked.main_rank >= min_rank || ranked.main_rank <= max_rank
+                ;
+
+            if (in_range)
+            {
+                if (side_index == 1)
+                {
+                    ranked.left_count++;
+                }
+                else if (side_index == 2)
+                {
+                    ranked.right_count++;
+                }
+            }
+        }
+    }
+
+    template <signed_size_type segment_identifier::*Member>
+    void find_polygons_for_source(signed_size_type the_index,
+                std::size_t start_index)
+    {
+        int state = 1; // 'closed', because start_index is "from", arrives at the turn
+        std::size_t last_from_rank = m_ranked_points[start_index].main_rank;
+        std::size_t previous_rank = m_ranked_points[start_index].main_rank;
+
+        for (std::size_t index = move<Member>(the_index, start_index);
+             ;
+             index = move<Member>(the_index, index))
+        {
+            rp& ranked = m_ranked_points[index];
+
+            if (ranked.main_rank != previous_rank && state == 0)
+            {
+                assign_ranks(last_from_rank, previous_rank - 1, 1);
+                assign_ranks(last_from_rank + 1, previous_rank, 2);
+            }
+
+            if (index == start_index)
+            {
+                return;
+            }
+
+            if (ranked.index == index_from)
+            {
+                last_from_rank = ranked.main_rank;
+                state++;
+            }
+            else if (ranked.index == index_to)
+            {
+                state--;
+            }
+
+            previous_rank = ranked.main_rank;
+        }
+    }
+};
+
+
+}}} // namespace detail::overlay::sort_by_side
+#endif //DOXYGEN_NO_DETAIL
+
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_SORT_BY_SIDE_HPP
diff --git a/boost/geometry/algorithms/detail/overlay/traverse.hpp b/boost/geometry/algorithms/detail/overlay/traverse.hpp
index 45e15d13d0..a8f4232124 100644
--- a/boost/geometry/algorithms/detail/overlay/traverse.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traverse.hpp
@@ -15,6 +15,7 @@
 
 #include <boost/geometry/algorithms/detail/overlay/backtrack_check_si.hpp>
 #include <boost/geometry/algorithms/detail/overlay/copy_segments.hpp>
+#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
 #include <boost/geometry/algorithms/num_points.hpp>
 #include <boost/geometry/core/access.hpp>
@@ -65,186 +66,740 @@ inline void debug_traverse(Turn const& , Operation, const char*)
 #endif
 
 
-template <typename Info, typename Turn>
-inline void set_visited_for_continue(Info& info, Turn const& turn)
+//! Metafunction to define side_order (clockwise, ccw) by operation_type
+template <operation_type OpType>
+struct side_compare {};
+
+template <>
+struct side_compare<operation_union>
 {
-    // On "continue", set "visited" for ALL directions
-    if (turn.operation == detail::overlay::operation_continue)
+    typedef std::greater<int> type;
+};
+
+template <>
+struct side_compare<operation_intersection>
+{
+    typedef std::less<int> type;
+};
+
+
+template
+<
+    bool Reverse1,
+    bool Reverse2,
+    operation_type OperationType,
+    typename Geometry1,
+    typename Geometry2,
+    typename Turns,
+    typename Clusters,
+    typename RobustPolicy,
+    typename Visitor,
+    typename Backtrack
+>
+struct traversal
+{
+    typedef typename side_compare<OperationType>::type side_compare_type;
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::turn_operation_type turn_operation_type;
+
+    typedef typename geometry::point_type<Geometry1>::type point_type;
+    typedef sort_by_side::side_sorter
+        <
+            Reverse1, Reverse2,
+            point_type, side_compare_type
+        > sbs_type;
+
+    inline traversal(Geometry1 const& geometry1, Geometry2 const& geometry2,
+            Turns& turns, Clusters const& clusters,
+            RobustPolicy const& robust_policy, Visitor& visitor)
+        : m_geometry1(geometry1)
+        , m_geometry2(geometry2)
+        , m_turns(turns)
+        , m_clusters(clusters)
+        , m_robust_policy(robust_policy)
+        , m_visitor(visitor)
+        , m_has_uu(false)
+        , m_has_only_uu(true)
+        , m_switch_at_uu(true)
+    {}
+
+
+    inline bool select_source(signed_size_type turn_index,
+                              segment_identifier const& seg_id1,
+                              segment_identifier const& seg_id2)
     {
-        for (typename boost::range_iterator
-            <
-                typename Info::container_type
-            >::type it = boost::begin(info.operations);
-            it != boost::end(info.operations);
-            ++it)
+        if (OperationType == operation_intersection)
         {
-            if (it->visited.none())
+            // For intersections always switch sources
+            return seg_id1.source_index != seg_id2.source_index;
+        }
+        else if (OperationType == operation_union)
+        {
+            // For uu, only switch sources if indicated
+            turn_type const& turn = m_turns[turn_index];
+
+            // TODO: pass this information
+            bool const is_buffer
+                    = turn.operations[0].seg_id.source_index
+                      == turn.operations[1].seg_id.source_index;
+
+            if (is_buffer)
             {
-                it->visited.set_visited();
+                // 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;
             }
+
+            // Temporarily use m_switch_at_uu, which does not solve all cases,
+            // but the majority of the more simple cases, making the interior
+            // rings valid
+            return m_switch_at_uu // turn.switch_source
+                    ? seg_id1.source_index != seg_id2.source_index
+                    : seg_id1.source_index == seg_id2.source_index;
         }
+        return false;
     }
-}
 
+    inline
+    signed_size_type get_next_turn_index(turn_operation_type const& op) const
+    {
+        return op.enriched.next_ip_index == -1
+                ? op.enriched.travels_to_ip_index
+                : op.enriched.next_ip_index;
+    }
 
-template
-<
-    bool Reverse1, bool Reverse2,
-    typename GeometryOut,
-    typename G1,
-    typename G2,
-    typename Turns,
-    typename IntersectionInfo,
-    typename RobustPolicy
->
-inline bool assign_next_ip(G1 const& g1, G2 const& g2,
-            Turns& turns,
-            typename boost::range_iterator<Turns>::type& ip,
-            GeometryOut& current_output,
-            IntersectionInfo& info,
-            segment_identifier& seg_id,
-            RobustPolicy const& robust_policy)
-{
-    info.visited.set_visited();
-    set_visited_for_continue(*ip, info);
+    inline bool traverse_possible(signed_size_type turn_index) const
+    {
+        if (turn_index == -1)
+        {
+            return false;
+        }
+
+        turn_type const& turn = m_turns[turn_index];
+
+        // It is not a dead end if there is an operation to continue, or of
+        // there is a cluster (assuming for now we can get out of the cluster)
+        return turn.cluster_id >= 0
+            || turn.has(OperationType)
+            || turn.has(operation_continue);
+    }
 
-    // If there is no next IP on this segment
-    if (info.enriched.next_ip_index < 0)
+    inline bool select_operation(turn_type& turn,
+                signed_size_type start_turn_index,
+                segment_identifier const& seg_id,
+                int& selected_op_index)
     {
-        if (info.enriched.travels_to_vertex_index < 0
-            || info.enriched.travels_to_ip_index < 0)
+        if (turn.discarded)
         {
             return false;
         }
 
-        BOOST_GEOMETRY_ASSERT(info.enriched.travels_to_vertex_index >= 0);
-        BOOST_GEOMETRY_ASSERT(info.enriched.travels_to_ip_index >= 0);
+        bool result = false;
+
+        typename turn_operation_type::comparable_distance_type
+                max_remaining_distance = 0;
 
-        if (info.seg_id.source_index == 0)
+        selected_op_index = -1;
+        for (int i = 0; i < 2; i++)
         {
-            geometry::copy_segments<Reverse1>(g1, info.seg_id,
-                    info.enriched.travels_to_vertex_index,
-                    robust_policy,
-                    current_output);
+            turn_operation_type const& op = turn.operations[i];
+            if (op.visited.started())
+            {
+                selected_op_index = i;
+                return true;
+            }
+
+            signed_size_type const next_turn_index = get_next_turn_index(op);
+
+            // In some cases there are two alternatives.
+            // For "ii", take the other one (alternate)
+            //           UNLESS the other one is already visited
+            // For "uu", take the same one (see above);
+            // For "cc", take either one, but if there is a starting one,
+            //           take that one. If next is dead end, skip that one.
+            if (   (op.operation == operation_continue
+                    && traverse_possible(next_turn_index)
+                    && ! result)
+                || (op.operation == OperationType
+                    && ! op.visited.finished()
+                    && (! result
+                        || select_source(next_turn_index, op.seg_id, seg_id)
+                        )
+                    )
+                )
+            {
+                if (op.operation == operation_continue)
+                {
+                    max_remaining_distance = op.remaining_distance;
+                }
+                selected_op_index = i;
+                debug_traverse(turn, op, " Candidate");
+                result = true;
+            }
+
+            if (op.operation == operation_continue && result)
+            {
+                if (next_turn_index == start_turn_index)
+                {
+                    selected_op_index = i;
+                    debug_traverse(turn, op, " Candidate override (start)");
+                }
+                else if (op.remaining_distance > max_remaining_distance)
+                {
+                    max_remaining_distance = op.remaining_distance;
+                    selected_op_index = i;
+                    debug_traverse(turn, op, " Candidate override (remaining)");
+                }
+            }
         }
-        else
+
+        if (result)
         {
-            geometry::copy_segments<Reverse2>(g2, info.seg_id,
-                    info.enriched.travels_to_vertex_index,
-                    robust_policy,
-                    current_output);
+           debug_traverse(turn, turn.operations[selected_op_index], "  Accepted");
         }
-        seg_id = info.seg_id;
-        ip = boost::begin(turns) + info.enriched.travels_to_ip_index;
+
+        return result;
     }
-    else
+
+    inline bool select_from_cluster(signed_size_type& turn_index,
+        int& op_index, signed_size_type start_turn_index,
+        sbs_type const& sbs, bool allow_pass_rank)
     {
-        ip = boost::begin(turns) + info.enriched.next_ip_index;
-        seg_id = info.seg_id;
+        bool const is_union = OperationType == operation_union;
+        bool const is_intersection = OperationType == operation_intersection;
+
+        std::size_t selected_rank = 0;
+        std::size_t min_rank = 0;
+        bool result = false;
+        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 (result && ranked_point.main_rank > selected_rank)
+            {
+                return result;
+            }
+
+            turn_type const& ranked_turn = m_turns[ranked_point.turn_index];
+            turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.op_index];
+
+            if (result && ranked_op.visited.finalized())
+            {
+                // One of the arcs in the same direction as the selected result
+                // is already traversed.
+                return false;
+            }
+
+            if (! allow_pass_rank && ranked_op.visited.finalized())
+            {
+                // Skip this one, go to next
+                min_rank = ranked_point.main_rank;
+                continue;
+            }
+
+            if (ranked_point.index == sort_by_side::index_to
+                && (ranked_point.main_rank > min_rank
+                    || ranked_turn.both(operation_continue)))
+            {
+                if ((is_union
+                     && 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)
+                    {
+                        // Don't override - only override if arrive at start
+                        continue;
+                    }
+
+                    turn_index = ranked_point.turn_index;
+                    op_index = ranked_point.op_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.op_index;
+                    }
+
+                    result = true;
+                    selected_rank = ranked_point.main_rank;
+                }
+                else if (! allow_pass_rank)
+                {
+                    return result;
+                }
+            }
+        }
+        return result;
     }
 
-    detail::overlay::append_no_dups_or_spikes(current_output, ip->point,
-        robust_policy);
+    inline bool select_turn_from_cluster(signed_size_type& turn_index,
+            int& op_index, signed_size_type start_turn_index,
+            point_type const& point)
+    {
+        bool const is_union = OperationType == operation_union;
 
-    return true;
-}
+        turn_type const& turn = m_turns[turn_index];
+        BOOST_ASSERT(turn.cluster_id >= 0);
 
+        typename Clusters::const_iterator mit = m_clusters.find(turn.cluster_id);
+        BOOST_ASSERT(mit != m_clusters.end());
 
-inline bool select_source(operation_type operation,
-                          signed_size_type source1,
-                          signed_size_type source2)
-{
-    return (operation == operation_intersection && source1 != source2)
-        || (operation == operation_union && source1 == source2)
-        ;
-}
+        std::set<signed_size_type> const& ids = mit->second;
 
+        sbs_type sbs;
+        sbs.set_origin(point);
 
-template
-<
-    typename Turn,
-    typename Iterator
->
-inline bool select_next_ip(operation_type operation,
-            Turn& turn,
-            segment_identifier const& seg_id,
-            Iterator& selected)
-{
-    if (turn.discarded)
+        for (typename std::set<signed_size_type>::const_iterator sit = ids.begin();
+             sit != ids.end(); ++sit)
+        {
+            signed_size_type cluster_turn_index = *sit;
+            turn_type const& cluster_turn = m_turns[cluster_turn_index];
+            if (cluster_turn.discarded)
+            {
+                // Defensive check, discarded turns should not be in cluster
+                continue;
+            }
+
+            for (int i = 0; i < 2; i++)
+            {
+                sbs.add(cluster_turn.operations[i], cluster_turn_index, i,
+                        m_geometry1, m_geometry2, false);
+            }
+        }
+
+        sbs.apply(turn.point);
+
+        int open_count = 0;
+        if (is_union)
+        {
+            // Check how many open spaces there are.
+            // TODO: might be moved to sbs itself, though it also uses turns
+
+            std::size_t last_rank = 0;
+            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.main_rank > last_rank
+                    && ranked_point.index == sort_by_side::index_to)
+                {
+                    turn_type const& ranked_turn = m_turns[ranked_point.turn_index];
+                    turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.op_index];
+                    if (ranked_op.enriched.count_left == 0
+                         && ranked_op.enriched.count_right > 0)
+                    {
+                        open_count++;
+                        last_rank = ranked_point.main_rank;
+                    }
+                }
+            }
+        }
+
+        bool allow = false;
+        if (open_count > 1)
+        {
+            sbs.reverse();
+            allow = true;
+        }
+
+        return select_from_cluster(turn_index, op_index, start_turn_index, sbs, allow);
+    }
+
+    inline void change_index_for_self_turn(signed_size_type& to_vertex_index,
+                turn_type const& start_turn,
+                turn_operation_type const& start_op,
+                int start_op_index)
     {
-        return false;
+        turn_operation_type const& other_op
+                = start_turn.operations[1 - start_op_index];
+        if (start_op.seg_id.source_index != other_op.seg_id.source_index)
+        {
+            // Not a buffer/self-turn
+            return;
+        }
+
+        // It travels to itself, can happen. If this is a buffer, it can
+        // sometimes travel to itself in the following configuration:
+        //
+        // +---->--+
+        // |       |
+        // |   +---*----+ *: one turn, with segment index 2/7
+        // |   |   |    |
+        // |   +---C    | C: closing point (start/end)
+        // |            |
+        // +------------+
+        //
+        // If it starts on segment 2 and travels to itself on segment 2, that
+        // should be corrected to 7 because that is the shortest path
+        //
+        // Also a uu turn (touching with another buffered ring) might have this
+        // apparent configuration, but there it should
+        // always travel the whole ring
+
+        bool const correct
+                = ! start_turn.both(operation_union)
+                  && start_op.seg_id.segment_index == to_vertex_index;
+
+#if defined(BOOST_GEOMETRY_DEBUG_TRAVERSE)
+        std::cout << " WARNING: self-buffer "
+                  << " correct=" << correct
+                  << " turn=" << operation_char(start_turn.operations[0].operation)
+                  << operation_char(start_turn.operations[1].operation)
+                  << " start=" << start_op.seg_id.segment_index
+                  << " from=" << to_vertex_index
+                  << " to=" << other_op.enriched.travels_to_vertex_index
+                  << std::endl;
+#endif
+
+        if (correct)
+        {
+            to_vertex_index = other_op.enriched.travels_to_vertex_index;
+        }
     }
 
-    bool has_tp = false;
+    template <typename Ring>
+    inline traverse_error_type travel_to_next_turn(signed_size_type start_turn_index,
+                int start_op_index,
+                signed_size_type& turn_index,
+                int& op_index,
+                segment_identifier& seg_id,
+                Ring& current_ring,
+                bool is_start)
+    {
+        int const previous_op_index = op_index;
+        signed_size_type const previous_turn_index = turn_index;
+        turn_type& previous_turn = m_turns[turn_index];
+        turn_operation_type& previous_op = previous_turn.operations[op_index];
 
-    typedef typename std::iterator_traits
-    <
-        Iterator
-    >::value_type operation_type;
+        // If there is no next IP on this segment
+        if (previous_op.enriched.next_ip_index < 0)
+        {
+            if (previous_op.enriched.travels_to_vertex_index < 0
+                || previous_op.enriched.travels_to_ip_index < 0)
+            {
+                return is_start
+                        ? traverse_error_no_next_ip_at_start
+                        : traverse_error_no_next_ip;
+            }
+
+            signed_size_type to_vertex_index = previous_op.enriched.travels_to_vertex_index;
+
+            if (is_start &&
+                    previous_op.enriched.travels_to_ip_index == start_turn_index)
+            {
+                change_index_for_self_turn(to_vertex_index, previous_turn,
+                    previous_op, start_op_index);
+            }
+
+            if (previous_op.seg_id.source_index == 0)
+            {
+                geometry::copy_segments<Reverse1>(m_geometry1,
+                        previous_op.seg_id, to_vertex_index,
+                        m_robust_policy, current_ring);
+            }
+            else
+            {
+                geometry::copy_segments<Reverse2>(m_geometry2,
+                        previous_op.seg_id, to_vertex_index,
+                        m_robust_policy, current_ring);
+            }
+            seg_id = previous_op.seg_id;
+            turn_index = previous_op.enriched.travels_to_ip_index;
+        }
+        else
+        {
+            turn_index = previous_op.enriched.next_ip_index;
+            seg_id = previous_op.seg_id;
+        }
+
+        // turn_index is not yet finally selected, can change for clusters
+        bool const has_cluster = m_turns[turn_index].cluster_id >= 0;
+        if (has_cluster)
+        {
+
+            if (! select_turn_from_cluster(turn_index, op_index,
+                    start_turn_index, current_ring.back()))
+            {
+                return is_start
+                    ? traverse_error_no_next_ip_at_start
+                    : traverse_error_no_next_ip;
+            }
+
+            if (is_start && turn_index == previous_turn_index)
+            {
+                op_index = previous_op_index;
+            }
+        }
+
+        turn_type& current_turn = m_turns[turn_index];
+        detail::overlay::append_no_dups_or_spikes(current_ring, current_turn.point,
+            m_robust_policy);
+
+        if (is_start)
+        {
+            // Register the start
+            previous_op.visited.set_started();
+            m_visitor.visit_traverse(m_turns, previous_turn, previous_op, "Start");
+        }
+
+        if (! has_cluster)
+        {
+            if (! select_operation(current_turn,
+                            start_turn_index,
+                            seg_id,
+                            op_index))
+            {
+                return is_start
+                    ? traverse_error_dead_end_at_start
+                    : traverse_error_dead_end;
+            }
+        }
+
+        turn_operation_type& op = current_turn.operations[op_index];
+        if (op.visited.finalized() || op.visited.visited())
+        {
+            return traverse_error_visit_again;
+        }
 
-    typename operation_type::comparable_distance_type
-            max_remaining_distance = 0;
+        // Register the visit
+        set_visited(current_turn, op);
+        m_visitor.visit_traverse(m_turns, current_turn, op, "Visit");
 
-    selected = boost::end(turn.operations);
-    for (Iterator it = boost::begin(turn.operations);
-        it != boost::end(turn.operations);
-        ++it)
+        return traverse_error_none;
+    }
+
+    inline void finalize_visit_info()
     {
-        if (it->visited.started())
+        for (typename boost::range_iterator<Turns>::type
+            it = boost::begin(m_turns);
+            it != boost::end(m_turns);
+            ++it)
         {
-            selected = it;
-            //std::cout << " RETURN";
-            return true;
+            turn_type& turn = *it;
+            for (int i = 0; i < 2; i++)
+            {
+                turn_operation_type& op = turn.operations[i];
+                op.visited.finalize();
+            }
         }
+    }
 
-        // In some cases there are two alternatives.
-        // For "ii", take the other one (alternate)
-        //           UNLESS the other one is already visited
-        // For "uu", take the same one (see above);
-        // For "cc", take either one, but if there is a starting one,
-        //           take that one.
-        if (   (it->operation == operation_continue
-                && (! has_tp || it->visited.started()
-                    )
-                )
-            || (it->operation == operation
-                && ! it->visited.finished()
-                && (! has_tp
-                    || select_source(operation,
-                            it->seg_id.source_index, seg_id.source_index)
-                    )
-                )
-            )
+    inline void set_visited(turn_type& turn, turn_operation_type& op)
+    {
+        // On "continue", set "visited" for ALL directions in this turn
+        if (op.operation == detail::overlay::operation_continue)
         {
-            if (it->operation == operation_continue)
+            for (int i = 0; i < 2; i++)
             {
-                max_remaining_distance = it->remaining_distance;
+                turn_operation_type& op = turn.operations[i];
+                if (op.visited.none())
+                {
+                    op.visited.set_visited();
+                }
             }
-            selected = it;
-            debug_traverse(turn, *it, " Candidate");
-            has_tp = true;
+        }
+        else
+        {
+            op.visited.set_visited();
+        }
+    }
+
+
+    template <typename Ring>
+    inline traverse_error_type traverse(Ring& ring,
+            signed_size_type start_turn_index, int start_op_index)
+    {
+        turn_type const& start_turn = m_turns[start_turn_index];
+        turn_operation_type& start_op = m_turns[start_turn_index].operations[start_op_index];
+
+        detail::overlay::append_no_dups_or_spikes(ring, start_turn.point,
+            m_robust_policy);
+
+        signed_size_type current_turn_index = start_turn_index;
+        int current_op_index = start_op_index;
+        segment_identifier current_seg_id;
+
+        traverse_error_type error = travel_to_next_turn(start_turn_index,
+                    start_op_index,
+                    current_turn_index, current_op_index, current_seg_id,
+                    ring, true);
+
+        if (error != traverse_error_none)
+        {
+            // This is not necessarily a problem, it happens for clustered turns
+            // which are "build in" or otherwise point inwards
+            return error;
+        }
+
+        if (current_turn_index == start_turn_index)
+        {
+            start_op.visited.set_finished();
+            m_visitor.visit_traverse(m_turns, m_turns[current_turn_index], start_op, "Early finish");
+            return traverse_error_none;
         }
 
-        if (it->operation == operation_continue && has_tp)
+        std::size_t const max_iterations = 2 + 2 * m_turns.size();
+        for (std::size_t i = 0; i <= max_iterations; i++)
         {
-            if (it->remaining_distance > max_remaining_distance)
+            // We assume clockwise polygons only, non self-intersecting, closed.
+            // However, the input might be different, and checking validity
+            // is up to the library user.
+
+            // Therefore we make here some sanity checks. If the input
+            // violates the assumptions, the output polygon will not be correct
+            // but the routine will stop and output the current polygon, and
+            // will continue with the next one.
+
+            // Below three reasons to stop.
+            error = travel_to_next_turn(start_turn_index, start_op_index,
+                    current_turn_index, current_op_index, current_seg_id,
+                    ring, false);
+
+            if (error != traverse_error_none)
             {
-                max_remaining_distance = it->remaining_distance;
-                selected = it;
-                debug_traverse(turn, *it, " Candidate override");
+                return error;
+            }
+
+            if (current_turn_index == start_turn_index
+                    && current_op_index == start_op_index)
+            {
+                start_op.visited.set_finished();
+                m_visitor.visit_traverse(m_turns, start_turn, start_op, "Finish");
+                return traverse_error_none;
             }
         }
+
+        return traverse_error_endless_loop;
     }
 
-    if (has_tp)
+    template <typename Rings>
+    void traverse_with_operation(turn_type const& start_turn,
+            std::size_t turn_index, int op_index,
+            Rings& rings, std::size_t& finalized_ring_size,
+            typename Backtrack::state_type& state)
     {
-       debug_traverse(turn, *selected, "  Accepted");
+        typedef typename boost::range_value<Rings>::type ring_type;
+
+        turn_operation_type const& start_op = start_turn.operations[op_index];
+
+        if (! start_op.visited.none()
+            || ! start_op.enriched.startable
+            || start_op.visited.rejected()
+            || ! (start_op.operation == OperationType
+                || start_op.operation == detail::overlay::operation_continue))
+        {
+            return;
+        }
+
+        ring_type ring;
+        traverse_error_type traverse_error = traverse(ring, turn_index, op_index);
+
+        if (traverse_error == traverse_error_none)
+        {
+            std::size_t const min_num_points
+                    = core_detail::closure::minimum_ring_size
+                            <
+                                geometry::closure<ring_type>::value
+                            >::value;
+
+            if (geometry::num_points(ring) >= min_num_points)
+            {
+                clean_closing_dups_and_spikes(ring, m_robust_policy);
+                rings.push_back(ring);
+
+                finalize_visit_info();
+                finalized_ring_size++;
+            }
+        }
+        else
+        {
+            Backtrack::apply(
+                finalized_ring_size,
+                rings, ring, m_turns, start_turn,
+                m_turns[turn_index].operations[op_index],
+                traverse_error,
+                m_geometry1, m_geometry2, m_robust_policy,
+                state, m_visitor);
+        }
     }
 
+    template <typename Rings>
+    void iterate(Rings& rings, std::size_t& finalized_ring_size,
+                 typename Backtrack::state_type& state,
+                 int pass)
+    {
+        if (pass == 1)
+        {
+            if (OperationType == operation_intersection)
+            {
+                // Second pass currently only used for uu
+                return;
+            }
+            if (! m_has_uu)
+            {
+                // There is no uu found in first pass
+                return;
+            }
+            if (m_has_only_uu)
+            {
+                m_switch_at_uu = false;
+            }
+        }
+
+        // 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];
 
-    return has_tp;
-}
+            if (start_turn.discarded || start_turn.blocked())
+            {
+                // Skip discarded and blocked turns
+                continue;
+            }
+            if (OperationType == 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;
+                    }
+                }
+                else
+                {
+                    m_has_only_uu = false;
+                }
+            }
+
+            for (int op_index = 0; op_index < 2; op_index++)
+            {
+                traverse_with_operation(start_turn, turn_index, op_index,
+                        rings, finalized_ring_size, state);
+            }
+        }
+    }
 
+private :
+    Geometry1 const& m_geometry1;
+    Geometry2 const& m_geometry2;
+    Turns& m_turns;
+    Clusters const& m_clusters;
+    RobustPolicy const& m_robust_policy;
+    Visitor& m_visitor;
+
+    // Next members are only used for operation union
+    bool m_has_uu;
+    bool m_has_only_uu;
+    bool m_switch_at_uu;
+};
 
 
 /*!
@@ -256,185 +811,45 @@ template
     bool Reverse1, bool Reverse2,
     typename Geometry1,
     typename Geometry2,
+    operation_type OperationType,
     typename Backtrack = backtrack_check_self_intersections<Geometry1, Geometry2>
 >
 class traverse
 {
 public :
-    template <typename RobustPolicy, typename Turns, typename Rings>
+    template
+    <
+        typename RobustPolicy,
+        typename Turns,
+        typename Rings,
+        typename Visitor,
+        typename Clusters
+    >
     static inline void apply(Geometry1 const& geometry1,
                 Geometry2 const& geometry2,
-                detail::overlay::operation_type operation,
                 RobustPolicy const& robust_policy,
-                Turns& turns, Rings& rings)
+                Turns& turns, Rings& rings,
+                Clusters const& clusters,
+                Visitor& visitor)
     {
-        typedef typename boost::range_value<Rings>::type ring_type;
-        typedef typename boost::range_iterator<Turns>::type turn_iterator;
-        typedef typename boost::range_value<Turns>::type turn_type;
-        typedef typename boost::range_iterator
+        traversal
             <
-                typename turn_type::container_type
-            >::type turn_operation_iterator_type;
+                Reverse1, Reverse2, OperationType,
+                Geometry1, Geometry2,
+                Turns, Clusters,
+                RobustPolicy, Visitor,
+                Backtrack
+            > trav(geometry1, geometry2, turns, clusters,
+                   robust_policy, visitor);
 
-        std::size_t const min_num_points
-                = core_detail::closure::minimum_ring_size
-                        <
-                            geometry::closure<ring_type>::value
-                        >::value;
-
-        std::size_t size_at_start = boost::size(rings);
+        std::size_t finalized_ring_size = boost::size(rings);
 
         typename Backtrack::state_type state;
-        do
-        {
-            state.reset();
 
-            // Iterate through all unvisited points
-            for (turn_iterator it = boost::begin(turns);
-                state.good() && it != boost::end(turns);
-                ++it)
-            {
-                // Skip discarded ones
-                if (! (it->discarded || ! it->selectable_start || it->blocked()))
-                {
-                    for (turn_operation_iterator_type iit = boost::begin(it->operations);
-                        state.good() && iit != boost::end(it->operations);
-                        ++iit)
-                    {
-                        if (iit->visited.none()
-                            && ! iit->visited.rejected()
-                            && (iit->operation == operation
-                                || iit->operation == detail::overlay::operation_continue)
-                            )
-                        {
-                            set_visited_for_continue(*it, *iit);
-
-                            ring_type current_output;
-                            detail::overlay::append_no_dups_or_spikes(current_output,
-                                it->point, robust_policy);
-
-                            turn_iterator current = it;
-                            turn_operation_iterator_type current_iit = iit;
-                            segment_identifier current_seg_id;
-
-                            if (! detail::overlay::assign_next_ip<Reverse1, Reverse2>(
-                                        geometry1, geometry2,
-                                        turns,
-                                        current, current_output,
-                                        *iit, current_seg_id,
-                                        robust_policy))
-                            {
-                                Backtrack::apply(
-                                    size_at_start,
-                                    rings, current_output, turns, *current_iit,
-                                    "No next IP",
-                                    geometry1, geometry2, robust_policy, state);
-                            }
-
-                            if (! detail::overlay::select_next_ip(
-                                            operation,
-                                            *current,
-                                            current_seg_id,
-                                            current_iit))
-                            {
-                                Backtrack::apply(
-                                    size_at_start,
-                                    rings, current_output, turns, *iit,
-                                    "Dead end at start",
-                                    geometry1, geometry2, robust_policy, state);
-                            }
-                            else
-                            {
-
-                                iit->visited.set_started();
-                                detail::overlay::debug_traverse(*it, *iit, "-> Started");
-                                detail::overlay::debug_traverse(*current, *current_iit, "Selected  ");
-
-
-                                typename boost::range_size<Turns>::type i = 0;
-
-                                while (current_iit != iit && state.good())
-                                {
-                                    if (current_iit->visited.visited())
-                                    {
-                                        // It visits a visited node again, without passing the start node.
-                                        // This makes it suspicious for endless loops
-                                        Backtrack::apply(
-                                            size_at_start,
-                                            rings,  current_output, turns, *iit,
-                                            "Visit again",
-                                            geometry1, geometry2, robust_policy, state);
-                                    }
-                                    else
-                                    {
-
-
-                                        // We assume clockwise polygons only, non self-intersecting, closed.
-                                        // However, the input might be different, and checking validity
-                                        // is up to the library user.
-
-                                        // Therefore we make here some sanity checks. If the input
-                                        // violates the assumptions, the output polygon will not be correct
-                                        // but the routine will stop and output the current polygon, and
-                                        // will continue with the next one.
-
-                                        // Below three reasons to stop.
-                                        detail::overlay::assign_next_ip<Reverse1, Reverse2>(
-                                            geometry1, geometry2,
-                                            turns, current, current_output,
-                                            *current_iit, current_seg_id,
-                                            robust_policy);
-
-                                        if (! detail::overlay::select_next_ip(
-                                                    operation,
-                                                    *current,
-                                                    current_seg_id,
-                                                    current_iit))
-                                        {
-                                            // Should not occur in valid (non-self-intersecting) polygons
-                                            // Should not occur in self-intersecting polygons without spikes
-                                            // Might occur in polygons with spikes
-                                            Backtrack::apply(
-                                                size_at_start,
-                                                rings,  current_output, turns, *iit,
-                                                "Dead end",
-                                                geometry1, geometry2, robust_policy, state);
-                                        }
-                                        else
-                                        {
-                                            detail::overlay::debug_traverse(*current, *current_iit, "Selected  ");
-                                        }
-
-                                        if (i++ > 2 + 2 * turns.size())
-                                        {
-                                            // Sanity check: there may be never more loops
-                                            // than turn points.
-                                            // Turn points marked as "ii" can be visited twice.
-                                            Backtrack::apply(
-                                                size_at_start,
-                                                rings,  current_output, turns, *iit,
-                                                "Endless loop",
-                                                geometry1, geometry2, robust_policy, state);
-                                        }
-                                    }
-                                }
-
-                                if (state.good())
-                                {
-                                    iit->visited.set_finished();
-                                    detail::overlay::debug_traverse(*current, *iit, "->Finished");
-                                    if (geometry::num_points(current_output) >= min_num_points)
-                                    {
-                                        clean_closing_dups_and_spikes(current_output, robust_policy);
-                                        rings.push_back(current_output);
-                                    }
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        } while (! state.good());
+        for (int pass = 0; pass < 2; pass++)
+        {
+            trav.iterate(rings, finalized_ring_size, state, pass);
+        }
     }
 };
 
diff --git a/boost/geometry/algorithms/detail/overlay/turn_info.hpp b/boost/geometry/algorithms/detail/overlay/turn_info.hpp
index 1ac77d5796..699997fc38 100644
--- a/boost/geometry/algorithms/detail/overlay/turn_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/turn_info.hpp
@@ -58,6 +58,8 @@ enum method_type
 template <typename Point, typename SegmentRatio>
 struct turn_operation
 {
+    typedef SegmentRatio segment_ratio_type;
+
     operation_type operation;
     segment_identifier seg_id;
     SegmentRatio fraction;
@@ -91,23 +93,25 @@ template
 struct turn_info
 {
     typedef Point point_type;
+    typedef SegmentRatio segment_ratio_type;
     typedef Operation turn_operation_type;
     typedef Container container_type;
 
     Point point;
     method_type method;
+    int cluster_id;
     bool discarded;
-    bool selectable_start; // Can be used as starting-turn in traverse
     bool colocated;
-
+    bool switch_source; // For u/u turns which can either switch or not
 
     Container operations;
 
     inline turn_info()
         : method(method_none)
+        , cluster_id(-1)
         , discarded(false)
-        , selectable_start(true)
         , colocated(false)
+        , switch_source(false)
     {}
 
     inline bool both(operation_type type) const
diff --git a/boost/geometry/algorithms/detail/overlay/visit_info.hpp b/boost/geometry/algorithms/detail/overlay/visit_info.hpp
index 4284a801a1..9e1e6b9056 100644
--- a/boost/geometry/algorithms/detail/overlay/visit_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/visit_info.hpp
@@ -28,11 +28,13 @@ private :
 
     int m_visit_code;
     bool m_rejected;
+    bool m_final;
 
 public:
     inline visit_info()
         : m_visit_code(0)
         , m_rejected(false)
+        , m_final(false)
     {}
 
     inline void set_visited() { m_visit_code = VISITED; }
@@ -49,15 +51,24 @@ public:
     inline bool started() const { return m_visit_code == STARTED; }
     inline bool finished() const { return m_visit_code == FINISHED; }
     inline bool rejected() const { return m_rejected; }
+    inline bool finalized() const { return m_final; }
 
     inline void clear()
     {
-        if (! rejected())
+        if (! m_rejected && ! m_final)
         {
             m_visit_code = NONE;
         }
     }
 
+    inline void finalize()
+    {
+        if (visited() || started() || finished() )
+        {
+            m_final = true;
+        }
+    }
+
 #ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
     friend std::ostream& operator<<(std::ostream &os, visit_info const& v)
     {