Imported Upstream version 1.67.0upstream/1.67.0

40 files changed, 1575 insertions, 1425 deletions

diff --git a/boost/geometry/algorithms/area.hpp b/boost/geometry/algorithms/area.hpp
index 18aea24036..c6e237e7cc 100644
--- a/boost/geometry/algorithms/area.hpp
+++ b/boost/geometry/algorithms/area.hpp
@@ -3,9 +3,10 @@
 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
-// This file was modified by Oracle on 2017.
-// Modifications copyright (c) 2017 Oracle and/or its affiliates.
+// This file was modified by Oracle on 2017, 2018.
+// Modifications copyright (c) 2017-2018 Oracle and/or its affiliates.
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
@@ -19,7 +20,6 @@
 #define BOOST_GEOMETRY_ALGORITHMS_AREA_HPP
 
 #include <boost/concept_check.hpp>
-#include <boost/mpl/if.hpp>
 #include <boost/range/functions.hpp>
 #include <boost/range/metafunctions.hpp>
 
@@ -43,7 +43,9 @@
 #include <boost/geometry/algorithms/detail/multi_sum.hpp>
 
 #include <boost/geometry/strategies/area.hpp>
+#include <boost/geometry/strategies/area_result.hpp>
 #include <boost/geometry/strategies/default_area_result.hpp>
+#include <boost/geometry/strategies/default_strategy.hpp>
 
 #include <boost/geometry/strategies/concepts/area_concept.hpp>
 
@@ -56,6 +58,7 @@
 namespace boost { namespace geometry
 {
 
+
 #ifndef DOXYGEN_NO_DETAIL
 namespace detail { namespace area
 {
@@ -83,10 +86,10 @@ template
 struct ring_area
 {
     template <typename Ring, typename Strategy>
-    static inline typename Strategy::return_type
+    static inline typename area_result<Ring, Strategy>::type
     apply(Ring const& ring, Strategy const& strategy)
     {
-        BOOST_CONCEPT_ASSERT( (geometry::concepts::AreaStrategy<Strategy>) );
+        BOOST_CONCEPT_ASSERT( (geometry::concepts::AreaStrategy<Ring, Strategy>) );
         assert_dimension<Ring, 2>();
 
         // Ignore warning (because using static method sometimes) on strategy
@@ -98,7 +101,7 @@ struct ring_area
         if (boost::size(ring)
                 < core_detail::closure::minimum_ring_size<Closure>::value)
         {
-            return typename Strategy::return_type();
+            return typename area_result<Ring, Strategy>::type();
         }
 
         typedef typename reversible_view<Ring const, Direction>::type rview_type;
@@ -110,7 +113,7 @@ struct ring_area
 
         rview_type rview(ring);
         view_type view(rview);
-        typename Strategy::state_type state;
+        typename Strategy::template state<Ring> state;
         iterator_type it = boost::begin(view);
         iterator_type end = boost::end(view);
 
@@ -144,9 +147,13 @@ template
 struct area : detail::calculate_null
 {
     template <typename Strategy>
-    static inline typename Strategy::return_type apply(Geometry const& geometry, Strategy const& strategy)
+    static inline typename area_result<Geometry, Strategy>::type
+        apply(Geometry const& geometry, Strategy const& strategy)
     {
-        return calculate_null::apply<typename Strategy::return_type>(geometry, strategy);
+        return calculate_null::apply
+            <
+                typename area_result<Geometry, Strategy>::type
+            >(geometry, strategy);
     }
 };
 
@@ -170,10 +177,11 @@ template <typename Polygon>
 struct area<Polygon, polygon_tag> : detail::calculate_polygon_sum
 {
     template <typename Strategy>
-    static inline typename Strategy::return_type apply(Polygon const& polygon, Strategy const& strategy)
+    static inline typename area_result<Polygon, Strategy>::type
+        apply(Polygon const& polygon, Strategy const& strategy)
     {
         return calculate_polygon_sum::apply<
-            typename Strategy::return_type,
+            typename area_result<Polygon, Strategy>::type,
             detail::area::ring_area
                 <
                     order_as_direction<geometry::point_order<Polygon>::value>::value,
@@ -188,12 +196,12 @@ template <typename MultiGeometry>
 struct area<MultiGeometry, multi_polygon_tag> : detail::multi_sum
 {
     template <typename Strategy>
-    static inline typename Strategy::return_type
+    static inline typename area_result<MultiGeometry, Strategy>::type
     apply(MultiGeometry const& multi, Strategy const& strategy)
     {
         return multi_sum::apply
                <
-                   typename Strategy::return_type,
+                   typename area_result<MultiGeometry, Strategy>::type,
                    area<typename boost::range_value<MultiGeometry>::type>
                >(multi, strategy);
     }
@@ -204,39 +212,73 @@ struct area<MultiGeometry, multi_polygon_tag> : detail::multi_sum
 #endif // DOXYGEN_NO_DISPATCH
 
 
-namespace resolve_variant {
+namespace resolve_strategy
+{
+
+struct area
+{
+    template <typename Geometry, typename Strategy>
+    static inline typename area_result<Geometry, Strategy>::type
+    apply(Geometry const& geometry, Strategy const& strategy)
+    {
+        return dispatch::area<Geometry>::apply(geometry, strategy);
+    }
+
+    template <typename Geometry>
+    static inline typename area_result<Geometry>::type
+    apply(Geometry const& geometry, default_strategy)
+    {
+        typedef typename strategy::area::services::default_strategy
+            <
+                typename cs_tag<Geometry>::type
+            >::type strategy_type;
+
+        return dispatch::area<Geometry>::apply(geometry, strategy_type());
+    }
+};
+
+
+} // namespace resolve_strategy
+
+
+namespace resolve_variant
+{
 
 template <typename Geometry>
 struct area
 {
     template <typename Strategy>
-    static inline typename Strategy::return_type apply(Geometry const& geometry,
-                                                       Strategy const& strategy)
+    static inline typename area_result<Geometry, Strategy>::type
+        apply(Geometry const& geometry, Strategy const& strategy)
     {
-        return dispatch::area<Geometry>::apply(geometry, strategy);
+        return resolve_strategy::area::apply(geometry, strategy);
     }
 };
 
 template <BOOST_VARIANT_ENUM_PARAMS(typename T)>
 struct area<boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)> >
 {
+    typedef boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)> variant_type;
+
     template <typename Strategy>
-    struct visitor: boost::static_visitor<typename Strategy::return_type>
+    struct visitor
+        : boost::static_visitor<typename area_result<variant_type, Strategy>::type>
     {
         Strategy const& m_strategy;
 
         visitor(Strategy const& strategy): m_strategy(strategy) {}
 
         template <typename Geometry>
-        typename Strategy::return_type operator()(Geometry const& geometry) const
+        typename area_result<variant_type, Strategy>::type
+        operator()(Geometry const& geometry) const
         {
             return area<Geometry>::apply(geometry, m_strategy);
         }
     };
 
     template <typename Strategy>
-    static inline typename Strategy::return_type
-    apply(boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)> const& geometry,
+    static inline typename area_result<variant_type, Strategy>::type
+    apply(variant_type const& geometry,
           Strategy const& strategy)
     {
         return boost::apply_visitor(visitor<Strategy>(strategy), geometry);
@@ -268,22 +310,14 @@ and Geographic as well.
 \qbk{[area] [area_output]}
 */
 template <typename Geometry>
-inline typename default_area_result<Geometry>::type area(Geometry const& geometry)
+inline typename area_result<Geometry>::type
+area(Geometry const& geometry)
 {
     concepts::check<Geometry const>();
 
-    // TODO put this into a resolve_strategy stage
-    //      (and take the return type from resolve_variant)
-    typedef typename point_type<Geometry>::type point_type;
-    typedef typename strategy::area::services::default_strategy
-        <
-            typename cs_tag<point_type>::type,
-            point_type
-        >::type strategy_type;
-
     // detail::throw_on_empty_input(geometry);
 
-    return resolve_variant::area<Geometry>::apply(geometry, strategy_type());
+    return resolve_variant::area<Geometry>::apply(geometry, default_strategy());
 }
 
 /*!
@@ -301,19 +335,19 @@ inline typename default_area_result<Geometry>::type area(Geometry const& geometr
 \qbk{
 [include reference/algorithms/area.qbk]
 
+[heading Available Strategies]
+\* [link geometry.reference.strategies.strategy_area_cartesian Cartesian]
+\* [link geometry.reference.strategies.strategy_area_spherical Spherical]
+\* [link geometry.reference.strategies.strategy_area_geographic Geographic]
+
 [heading Example]
 [area_with_strategy]
 [area_with_strategy_output]
-
-[heading Available Strategies]
-\* [link geometry.reference.strategies.strategy_area_surveyor Surveyor (cartesian)]
-\* [link geometry.reference.strategies.strategy_area_spherical Spherical]
-[/link geometry.reference.strategies.strategy_area_geographic Geographic]
 }
  */
 template <typename Geometry, typename Strategy>
-inline typename Strategy::return_type area(
-        Geometry const& geometry, Strategy const& strategy)
+inline typename area_result<Geometry, Strategy>::type
+area(Geometry const& geometry, Strategy const& strategy)
 {
     concepts::check<Geometry const>();
 
diff --git a/boost/geometry/algorithms/convert.hpp b/boost/geometry/algorithms/convert.hpp
index 6a8ba1acb6..6ccaa0dd4b 100644
--- a/boost/geometry/algorithms/convert.hpp
+++ b/boost/geometry/algorithms/convert.hpp
@@ -5,6 +5,10 @@
 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
 // Copyright (c) 2014 Adam Wulkiewicz, Lodz, Poland.
 
+// This file was modified by Oracle on 2017.
+// Modifications copyright (c) 2017, Oracle and/or its affiliates.
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
 // Parts of Boost.Geometry are redesigned from Geodan's Geographic Library
 // (geolib/GGL), copyright (c) 1995-2010 Geodan, Amsterdam, the Netherlands.
 
@@ -29,10 +33,8 @@
 
 #include <boost/geometry/arithmetic/arithmetic.hpp>
 #include <boost/geometry/algorithms/not_implemented.hpp>
-#include <boost/geometry/algorithms/append.hpp>
 #include <boost/geometry/algorithms/clear.hpp>
 #include <boost/geometry/algorithms/for_each.hpp>
-#include <boost/geometry/algorithms/detail/assign_values.hpp>
 #include <boost/geometry/algorithms/detail/assign_box_corners.hpp>
 #include <boost/geometry/algorithms/detail/assign_indexed_point.hpp>
 #include <boost/geometry/algorithms/detail/convert_point_to_point.hpp>
@@ -153,8 +155,24 @@ struct range_to_range
             geometry::closure<Range1>::value
         >::type view_type;
 
+    struct default_policy
+    {
+        template <typename Point1, typename Point2>
+        static inline void apply(Point1 const& point1, Point2 & point2)
+        {
+            geometry::detail::conversion::convert_point_to_point(point1, point2);
+        }
+    };
+    
     static inline void apply(Range1 const& source, Range2& destination)
     {
+        apply(source, destination, default_policy());
+    }
+
+    template <typename ConvertPointPolicy>
+    static inline ConvertPointPolicy apply(Range1 const& source, Range2& destination,
+                                           ConvertPointPolicy convert_point)
+    {
         geometry::clear(destination);
 
         rview_type rview(source);
@@ -179,8 +197,12 @@ struct range_to_range
             it != boost::end(view) && i < n;
             ++it, ++i)
         {
-            geometry::append(destination, *it);
+            typename boost::range_value<Range2>::type point;
+            convert_point.apply(*it, point);
+            range::push_back(destination, point);
         }
+
+        return convert_point;
     }
 };
 
diff --git a/boost/geometry/algorithms/correct.hpp b/boost/geometry/algorithms/correct.hpp
index 07e012fafa..21b9977ce5 100644
--- a/boost/geometry/algorithms/correct.hpp
+++ b/boost/geometry/algorithms/correct.hpp
@@ -3,7 +3,7 @@
 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
 // Copyright (c) 2008-2012 Bruno Lalande, Paris, France.
 // Copyright (c) 2009-2012 Mateusz Loskot, London, UK.
-// Copyright (c) 2014 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2014-2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -144,7 +144,7 @@ struct correct_ring
         detail::correct_closure::close_or_open_ring<Ring>::apply(r);
 
         // Check area
-        typedef typename Strategy::return_type area_result_type;
+        typedef typename area_result<Ring, Strategy>::type area_result_type;
         Predicate<area_result_type> predicate;
         area_result_type const zero = 0;
         if (predicate(ring_area_type::apply(r, strategy), zero))
@@ -322,8 +322,7 @@ inline void correct(Geometry& geometry)
 
     typedef typename strategy::area::services::default_strategy
         <
-            typename cs_tag<point_type>::type,
-            point_type
+            typename cs_tag<point_type>::type
         >::type strategy_type;
 
     resolve_variant::correct<Geometry>::apply(geometry, strategy_type());
diff --git a/boost/geometry/algorithms/densify.hpp b/boost/geometry/algorithms/densify.hpp
new file mode 100644
index 0000000000..ed948b3464
--- /dev/null
+++ b/boost/geometry/algorithms/densify.hpp
@@ -0,0 +1,425 @@
+// Boost.Geometry
+
+// Copyright (c) 2017-2018, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#ifndef BOOST_GEOMETRY_ALGORITHMS_DENSIFY_HPP
+#define BOOST_GEOMETRY_ALGORITHMS_DENSIFY_HPP
+
+
+#include <boost/geometry/algorithms/clear.hpp>
+#include <boost/geometry/algorithms/detail/convert_point_to_point.hpp>
+#include <boost/geometry/algorithms/not_implemented.hpp>
+#include <boost/geometry/core/closure.hpp>
+#include <boost/geometry/core/exception.hpp>
+#include <boost/geometry/core/point_type.hpp>
+#include <boost/geometry/core/tag.hpp>
+#include <boost/geometry/core/tags.hpp>
+#include <boost/geometry/strategies/default_strategy.hpp>
+#include <boost/geometry/strategies/densify.hpp>
+#include <boost/geometry/util/condition.hpp>
+#include <boost/geometry/util/range.hpp>
+
+#include <boost/range/size.hpp>
+#include <boost/range/value_type.hpp>
+
+#include <boost/throw_exception.hpp>
+
+
+namespace boost { namespace geometry
+{
+
+
+#ifndef DOXYGEN_NO_DETAIL
+namespace detail { namespace densify
+{
+
+template <typename Range>
+struct push_back_policy
+{
+    typedef typename boost::range_value<Range>::type point_type;
+
+    inline explicit push_back_policy(Range & rng)
+        : m_rng(rng)
+    {}
+
+    inline void apply(point_type const& p)
+    {
+        range::push_back(m_rng, p);
+    }
+
+private:
+    Range & m_rng;
+};
+
+template <typename Range, typename Point>
+inline void convert_and_push_back(Range & range, Point const& p)
+{
+    typename boost::range_value<Range>::type p2;
+    geometry::detail::conversion::convert_point_to_point(p, p2);
+    range::push_back(range, p2);
+}
+
+template <bool AppendLastPoint = true>
+struct densify_range
+{
+    template <typename FwdRng, typename MutRng, typename T, typename Strategy>
+    static inline void apply(FwdRng const& rng, MutRng & rng_out,
+                             T const& len, Strategy const& strategy)
+    {
+        typedef typename boost::range_iterator<FwdRng const>::type iterator_t;
+        typedef typename boost::range_value<FwdRng>::type point_t;
+
+        iterator_t it = boost::begin(rng);
+        iterator_t end = boost::end(rng);
+
+        if (it == end) // empty(rng)
+        {
+            return;
+        }
+            
+        push_back_policy<MutRng> policy(rng_out);
+
+        iterator_t prev = it;
+        for ( ++it ; it != end ; prev = it++)
+        {
+            point_t const& p0 = *prev;
+            point_t const& p1 = *it;
+
+            convert_and_push_back(rng_out, p0);
+
+            strategy.apply(p0, p1, policy, len);
+        }
+
+        if (BOOST_GEOMETRY_CONDITION(AppendLastPoint))
+        {
+            convert_and_push_back(rng_out, *prev); // back(rng)
+        }
+    }
+};
+
+template <bool IsClosed1, bool IsClosed2> // false, X
+struct densify_ring
+{
+    template <typename Geometry, typename GeometryOut, typename T, typename Strategy>
+    static inline void apply(Geometry const& ring, GeometryOut & ring_out,
+                             T const& len, Strategy const& strategy)
+    {
+        geometry::detail::densify::densify_range<true>
+            ::apply(ring, ring_out, len, strategy);
+
+        if (boost::size(ring) <= 1)
+            return;
+
+        typedef typename point_type<Geometry>::type point_t;
+        point_t const& p0 = range::back(ring);
+        point_t const& p1 = range::front(ring);
+
+        push_back_policy<GeometryOut> policy(ring_out);
+
+        strategy.apply(p0, p1, policy, len);
+
+        if (BOOST_GEOMETRY_CONDITION(IsClosed2))
+        {
+            convert_and_push_back(ring_out, p1);
+        }
+    }
+};
+
+template <>
+struct densify_ring<true, true>
+    : densify_range<true>
+{};
+
+template <>
+struct densify_ring<true, false>
+    : densify_range<false>
+{};
+
+
+}} // namespace detail::densify
+#endif // DOXYGEN_NO_DETAIL
+
+
+#ifndef DOXYGEN_NO_DISPATCH
+namespace dispatch
+{
+
+
+template
+<
+    typename Geometry,
+    typename GeometryOut,
+    typename Tag1 = typename tag<Geometry>::type,
+    typename Tag2 = typename tag<GeometryOut>::type
+>
+struct densify
+    : not_implemented<Tag1, Tag2>
+{};
+
+template <typename Geometry, typename GeometryOut>
+struct densify<Geometry, GeometryOut, linestring_tag, linestring_tag>
+    : geometry::detail::densify::densify_range<>
+{};
+
+template <typename Geometry, typename GeometryOut>
+struct densify<Geometry, GeometryOut, multi_linestring_tag, multi_linestring_tag>
+{
+    template <typename T, typename Strategy>
+    static void apply(Geometry const& mls, GeometryOut & mls_out,
+                      T const& len, Strategy const& strategy)
+    {
+        std::size_t count = boost::size(mls);
+        range::resize(mls_out, count);
+
+        for (std::size_t i = 0 ; i < count ; ++i)
+        {
+            geometry::detail::densify::densify_range<>
+                ::apply(range::at(mls, i), range::at(mls_out, i),
+                        len, strategy);
+        }
+    }
+};
+
+template <typename Geometry, typename GeometryOut>
+struct densify<Geometry, GeometryOut, ring_tag, ring_tag>
+    : geometry::detail::densify::densify_ring
+        <
+            geometry::closure<Geometry>::value != geometry::open,
+            geometry::closure<GeometryOut>::value != geometry::open
+        >
+{};
+
+template <typename Geometry, typename GeometryOut>
+struct densify<Geometry, GeometryOut, polygon_tag, polygon_tag>
+{
+    template <typename T, typename Strategy>
+    static void apply(Geometry const& poly, GeometryOut & poly_out,
+                      T const& len, Strategy const& strategy)
+    {
+        apply_ring(exterior_ring(poly), exterior_ring(poly_out),
+                   len, strategy);
+
+        std::size_t count = boost::size(interior_rings(poly));
+        range::resize(interior_rings(poly_out), count);
+
+        for (std::size_t i = 0 ; i < count ; ++i)
+        {
+            apply_ring(range::at(interior_rings(poly), i),
+                       range::at(interior_rings(poly_out), i),
+                       len, strategy);
+        }
+    }
+
+    template <typename Ring, typename RingOut, typename T, typename Strategy>
+    static void apply_ring(Ring const& ring, RingOut & ring_out,
+                           T const& len, Strategy const& strategy)
+    {
+        densify<Ring, RingOut, ring_tag, ring_tag>
+            ::apply(ring, ring_out, len, strategy);
+    }
+};
+
+template <typename Geometry, typename GeometryOut>
+struct densify<Geometry, GeometryOut, multi_polygon_tag, multi_polygon_tag>
+{
+    template <typename T, typename Strategy>
+    static void apply(Geometry const& mpoly, GeometryOut & mpoly_out,
+                      T const& len, Strategy const& strategy)
+    {
+        std::size_t count = boost::size(mpoly);
+        range::resize(mpoly_out, count);
+
+        for (std::size_t i = 0 ; i < count ; ++i)
+        {
+            apply_poly(range::at(mpoly, i),
+                       range::at(mpoly_out, i),
+                       len, strategy);
+        }
+    }
+
+    template <typename Poly, typename PolyOut, typename T, typename Strategy>
+    static void apply_poly(Poly const& poly, PolyOut & poly_out,
+                           T const& len, Strategy const& strategy)
+    {
+        densify<Poly, PolyOut, polygon_tag, polygon_tag>::
+            apply(poly, poly_out, len, strategy);
+    }
+};
+
+
+} // namespace dispatch
+#endif // DOXYGEN_NO_DISPATCH
+
+
+namespace resolve_strategy
+{
+
+struct densify
+{
+    template <typename Geometry, typename Distance, typename Strategy>
+    static inline void apply(Geometry const& geometry,
+                             Geometry& out,
+                             Distance const& max_distance,
+                             Strategy const& strategy)
+    {
+        dispatch::densify<Geometry, Geometry>
+            ::apply(geometry, out, max_distance, strategy);
+    }
+
+    template <typename Geometry, typename Distance>
+    static inline void apply(Geometry const& geometry,
+                             Geometry& out,
+                             Distance const& max_distance,
+                             default_strategy)
+    {
+        typedef typename strategy::densify::services::default_strategy
+            <
+                typename cs_tag<Geometry>::type
+            >::type strategy_type;
+        
+        /*BOOST_CONCEPT_ASSERT(
+            (concepts::DensifyStrategy<strategy_type>)
+        );*/
+
+        apply(geometry, out, max_distance, strategy_type());
+    }
+};
+
+} // namespace resolve_strategy
+
+
+namespace resolve_variant {
+
+template <typename Geometry>
+struct densify
+{
+    template <typename Distance, typename Strategy>
+    static inline void apply(Geometry const& geometry,
+                             Geometry& out,
+                             Distance const& max_distance,
+                             Strategy const& strategy)
+    {
+        resolve_strategy::densify::apply(geometry, out, max_distance, strategy);
+    }
+};
+
+template <BOOST_VARIANT_ENUM_PARAMS(typename T)>
+struct densify<boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)> >
+{
+    template <typename Distance, typename Strategy>
+    struct visitor: boost::static_visitor<void>
+    {
+        Distance const& m_max_distance;
+        Strategy const& m_strategy;
+
+        visitor(Distance const& max_distance, Strategy const& strategy)
+            : m_max_distance(max_distance)
+            , m_strategy(strategy)
+        {}
+
+        template <typename Geometry>
+        void operator()(Geometry const& geometry, Geometry& out) const
+        {
+            densify<Geometry>::apply(geometry, out, m_max_distance, m_strategy);
+        }
+    };
+
+    template <typename Distance, typename Strategy>
+    static inline void
+    apply(boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)> const& geometry,
+          boost::variant<BOOST_VARIANT_ENUM_PARAMS(T)>& out,
+          Distance const& max_distance,
+          Strategy const& strategy)
+    {
+        boost::apply_visitor(
+            visitor<Distance, Strategy>(max_distance, strategy),
+            geometry,
+            out
+        );
+    }
+};
+
+} // namespace resolve_variant
+
+
+/*!
+\brief Densify a geometry using a specified strategy
+\ingroup densify
+\tparam Geometry \tparam_geometry
+\tparam Distance A numerical distance measure
+\tparam Strategy A type fulfilling a DensifyStrategy concept
+\param geometry Input geometry, to be densified
+\param out Output geometry, densified version of the input geometry
+\param max_distance Distance threshold (in units depending on strategy)
+\param strategy Densify strategy to be used for densification
+
+\qbk{distinguish,with strategy}
+\qbk{[include reference/algorithms/densify.qbk]}
+
+\qbk{
+[heading Available Strategies]
+\* [link geometry.reference.strategies.strategy_densify_cartesian Cartesian]
+\* [link geometry.reference.strategies.strategy_densify_spherical Spherical]
+\* [link geometry.reference.strategies.strategy_densify_geographic Geographic]
+
+[heading Example]
+[densify_strategy]
+[densify_strategy_output]
+}
+*/
+template <typename Geometry, typename Distance, typename Strategy>
+inline void densify(Geometry const& geometry,
+                    Geometry& out,
+                    Distance const& max_distance,
+                    Strategy const& strategy)
+{
+    concepts::check<Geometry>();
+
+    if (max_distance <= Distance(0))
+    {
+        BOOST_THROW_EXCEPTION(geometry::invalid_input_exception());
+    }
+
+    geometry::clear(out);
+
+    resolve_variant::densify
+        <
+            Geometry
+        >::apply(geometry, out, max_distance, strategy);
+}
+
+
+/*!
+\brief Densify a geometry
+\ingroup densify
+\tparam Geometry \tparam_geometry
+\tparam Distance A numerical distance measure
+\param geometry Input geometry, to be densified
+\param out Output geometry, densified version of the input geometry
+\param max_distance Distance threshold (in units depending on coordinate system)
+
+\qbk{[include reference/algorithms/densify.qbk]}
+
+\qbk{
+[heading Example]
+[densify]
+[densify_output]
+}
+*/
+template <typename Geometry, typename Distance>
+inline void densify(Geometry const& geometry,
+                    Geometry& out,
+                    Distance const& max_distance)
+{
+    densify(geometry, out, max_distance, default_strategy());
+}
+
+
+}} // namespace boost::geometry
+
+#endif // BOOST_GEOMETRY_ALGORITHMS_DENSIFY_HPP
diff --git a/boost/geometry/algorithms/detail/azimuth.hpp b/boost/geometry/algorithms/detail/azimuth.hpp
index a5863d7d24..21007778bb 100644
--- a/boost/geometry/algorithms/detail/azimuth.hpp
+++ b/boost/geometry/algorithms/detail/azimuth.hpp
@@ -15,18 +15,22 @@
 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_AZIMUTH_HPP
 
+
+#include <boost/geometry/algorithms/not_implemented.hpp>
+
 #include <boost/geometry/core/cs.hpp>
 #include <boost/geometry/core/access.hpp>
 #include <boost/geometry/core/radian_access.hpp>
 #include <boost/geometry/core/tags.hpp>
 
-#include <boost/geometry/util/math.hpp>
-
-#include <boost/geometry/algorithms/not_implemented.hpp>
-
 #include <boost/geometry/formulas/spherical.hpp>
 #include <boost/geometry/formulas/vincenty_inverse.hpp>
 
+#include <boost/geometry/srs/spheroid.hpp>
+
+#include <boost/geometry/util/math.hpp>
+
+
 namespace boost { namespace geometry
 {
 
diff --git a/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp b/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp
index 990081a86c..a969b21cfb 100644
--- a/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp
+++ b/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp
@@ -41,10 +41,6 @@
 
 #include <boost/geometry/views/detail/normalized_view.hpp>
 
-#if defined(BOOST_GEOMETRY_BUFFER_SIMPLIFY_WITH_AX)
-#include <boost/geometry/strategies/cartesian/distance_projected_point_ax.hpp>
-#endif
-
 
 namespace boost { namespace geometry
 {
@@ -67,45 +63,21 @@ inline void simplify_input(Range const& range,
     // sensitive to small scale input features, however the result will
     // look better.
     // It also gets rid of duplicate points
-#if ! defined(BOOST_GEOMETRY_BUFFER_SIMPLIFY_WITH_AX)
-    geometry::simplify(range, simplified, distance.simplify_distance());
-#else
 
-    typedef typename boost::range_value<Range>::type point_type;
-    typedef strategy::distance::detail::projected_point_ax<> ax_type;
-    typedef typename strategy::distance::services::return_type
+    typedef typename geometry::point_type<Range>::type point_type;
+    typedef typename strategy::distance::services::default_strategy
     <
-        strategy::distance::detail::projected_point_ax<>,
-        point_type,
-        point_type
-    >::type return_type;
-
-    typedef strategy::distance::detail::projected_point_ax_less
+        point_tag, segment_tag, point_type
+    >::type ds_strategy_type;
+    typedef strategy::simplify::douglas_peucker
     <
-        return_type
-    > comparator_type;
+        point_type, ds_strategy_type
+    > strategy_type;
+
+    geometry::detail::simplify::simplify_range<2>::apply(range,
+        simplified, distance.simplify_distance(),
+        strategy_type());
 
-    typedef strategy::simplify::detail::douglas_peucker
-    <
-        point_type,
-        strategy::distance::detail::projected_point_ax<>,
-        comparator_type
-    > dp_ax;
-
-    return_type max_distance(distance.simplify_distance() * 2.0,
-                             distance.simplify_distance());
-    comparator_type comparator(max_distance);
-    dp_ax strategy(comparator);
-
-    geometry::simplify(range, simplified, max_distance, strategy);
-#endif
-
-    if (boost::size(simplified) == 2
-        && geometry::equals(geometry::range::front(simplified),
-                geometry::range::back(simplified)))
-    {
-        traits::resize<Range>::apply(simplified, 1);
-    }
 }
 
 
diff --git a/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp b/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp
index 92dcdcc7b0..0374b53a99 100644
--- a/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp
+++ b/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp
@@ -127,6 +127,10 @@ public :
     void visit_traverse_reject(Turns const& , Turn const& , Operation const& ,
             detail::overlay::traverse_error_type )
     {}
+
+    template <typename Rings>
+    void visit_generated_rings(Rings const& )
+    {}
 };
 
 
diff --git a/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp b/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp
index bc9c1ca7fb..ba824243cc 100644
--- a/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp
+++ b/boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2016-2017.
 // Modifications copyright (c) 2016-2017 Oracle and/or its affiliates.
@@ -155,8 +156,14 @@ struct buffered_piece_collection
             robust_point_type
         >::type robust_area_strategy_type;
 
-    typedef typename area_strategy_type::return_type area_result_type;
-    typedef typename robust_area_strategy_type::return_type robust_area_result_type;
+    typedef typename area_strategy_type::template result_type
+        <
+            point_type
+        >::type area_result_type;
+    typedef typename robust_area_strategy_type::template result_type
+        <
+            robust_point_type
+        >::type robust_area_result_type;
 
     typedef typename geometry::rescale_policy_type
         <
@@ -603,7 +610,11 @@ struct buffered_piece_collection
                 if (multi0 == multi1)
                 {
                     const deflate_properties& prop =  properties[multi0];
-                    if (! prop.has_inflated && prop.count < 3)
+
+                    // NOTE: Keep brackets around prop.count
+                    // avoid gcc-bug "parse error in template argument list"
+                    // GCC versions 5.4 and 5.5 (and probably more)
+                    if (! prop.has_inflated && (prop.count) < 3)
                     {
                         // Property is not inflated
                         // Not enough points, this might be caused by <float> where
@@ -1570,10 +1581,8 @@ struct buffered_piece_collection
             }
         }
 
-        // Assign parents, checking orientation but NOT discarding double
-        // negative rings (negative child with negative parent)
-        detail::overlay::assign_parents(offsetted_rings, traversed_rings,
-                selected, m_intersection_strategy, true, false);
+        detail::overlay::assign_parents<overlay_buffer>(offsetted_rings, traversed_rings,
+                selected, m_intersection_strategy);
         return detail::overlay::add_rings<GeometryOutput>(selected, offsetted_rings, traversed_rings, out,
                                                           m_area_strategy);
     }
diff --git a/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp b/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp
index 29e49f9dae..eb6fc02c8c 100644
--- a/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp
+++ b/boost/geometry/algorithms/detail/buffer/turn_in_piece_visitor.hpp
@@ -36,6 +36,8 @@
 
 #if defined(BOOST_GEOMETRY_BUFFER_USE_SIDE_OF_INTERSECTION)
 #include <boost/geometry/strategies/cartesian/side_of_intersection.hpp>
+#else
+#include <boost/geometry/strategies/agnostic/point_in_poly_winding.hpp>
 #endif
 
 
diff --git a/boost/geometry/algorithms/detail/distance/multipoint_to_geometry.hpp b/boost/geometry/algorithms/detail/distance/multipoint_to_geometry.hpp
index ccd3aa462e..8f1b33cfb2 100644
--- a/boost/geometry/algorithms/detail/distance/multipoint_to_geometry.hpp
+++ b/boost/geometry/algorithms/detail/distance/multipoint_to_geometry.hpp
@@ -193,6 +193,7 @@ struct distance
         >
 {};
 
+
 template <typename MultiPoint, typename Linear, typename Strategy>
 struct distance
     <
diff --git a/boost/geometry/algorithms/detail/envelope/segment.hpp b/boost/geometry/algorithms/detail/envelope/segment.hpp
index 97f2fc84a8..06c64bafc8 100644
--- a/boost/geometry/algorithms/detail/envelope/segment.hpp
+++ b/boost/geometry/algorithms/detail/envelope/segment.hpp
@@ -27,7 +27,6 @@
 #include <boost/geometry/core/coordinate_system.hpp>
 #include <boost/geometry/core/coordinate_type.hpp>
 #include <boost/geometry/core/cs.hpp>
-#include <boost/geometry/core/srs.hpp>
 #include <boost/geometry/core/point_type.hpp>
 #include <boost/geometry/core/radian_access.hpp>
 #include <boost/geometry/core/tags.hpp>
diff --git a/boost/geometry/algorithms/detail/extreme_points.hpp b/boost/geometry/algorithms/detail/extreme_points.hpp
index 61e984ee3c..607997813c 100644
--- a/boost/geometry/algorithms/detail/extreme_points.hpp
+++ b/boost/geometry/algorithms/detail/extreme_points.hpp
@@ -3,7 +3,7 @@
 // Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands.
 // Copyright (c) 2008-2013 Bruno Lalande, Paris, France.
 // Copyright (c) 2009-2013 Mateusz Loskot, London, UK.
-// Copyright (c) 2013-2014 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -536,6 +536,24 @@ inline bool extreme_points(Geometry const& geometry,
 }
 
 
+template
+<
+    std::size_t Edge,
+    typename Geometry,
+    typename Extremes,
+    typename Intruders
+>
+inline bool extreme_points(Geometry const& geometry,
+                           Extremes& extremes,
+                           Intruders& intruders)
+{
+    typedef typename strategy::side::services::default_strategy
+            <
+                typename cs_tag<Geometry>::type
+            >::type strategy_type;
+
+    return geometry::extreme_points<Edge>(geometry,extremes, intruders, strategy_type());
+}
 
 }} // namespace boost::geometry
 
diff --git a/boost/geometry/algorithms/detail/is_valid/interface.hpp b/boost/geometry/algorithms/detail/is_valid/interface.hpp
index ee013377c4..e7f5c5783e 100644
--- a/boost/geometry/algorithms/detail/is_valid/interface.hpp
+++ b/boost/geometry/algorithms/detail/is_valid/interface.hpp
@@ -1,6 +1,6 @@
-// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Boost.Geometry
 
-// Copyright (c) 2014-2017, Oracle and/or its affiliates.
+// Copyright (c) 2014-2018, Oracle and/or its affiliates.
 
 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
@@ -18,9 +18,9 @@
 #include <boost/variant/static_visitor.hpp>
 #include <boost/variant/variant_fwd.hpp>
 
-#include <boost/geometry/geometries/concepts/check.hpp>
-
 #include <boost/geometry/algorithms/dispatch/is_valid.hpp>
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/geometries/concepts/check.hpp>
 #include <boost/geometry/policies/is_valid/default_policy.hpp>
 #include <boost/geometry/policies/is_valid/failing_reason_policy.hpp>
 #include <boost/geometry/policies/is_valid/failure_type_policy.hpp>
diff --git a/boost/geometry/algorithms/detail/is_valid/ring.hpp b/boost/geometry/algorithms/detail/is_valid/ring.hpp
index 0b95950430..40698155b5 100644
--- a/boost/geometry/algorithms/detail/is_valid/ring.hpp
+++ b/boost/geometry/algorithms/detail/is_valid/ring.hpp
@@ -1,5 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
+
 // Copyright (c) 2014-2017, Oracle and/or its affiliates.
 
 // Contributed and/or modified by Menelaos Karavelas, on behalf of Oracle
@@ -107,8 +109,6 @@ struct is_properly_oriented
     {
         boost::ignore_unused(visitor);
 
-        typedef typename point_type<Ring>::type point_type;
-
         typedef detail::area::ring_area
             <
                 order_as_direction<geometry::point_order<Ring>::value>::value,
@@ -117,8 +117,8 @@ struct is_properly_oriented
 
         typedef typename Strategy::template area_strategy
             <
-                point_type
-            >::type::return_type area_result_type;
+                Ring
+            >::type::template result_type<Ring>::type area_result_type;
 
         typename ring_area_predicate
             <
@@ -129,7 +129,7 @@ struct is_properly_oriented
         area_result_type const zero = 0;
         area_result_type const area
             = ring_area_type::apply(ring,
-                                    strategy.template get_area_strategy<point_type>());
+                                    strategy.template get_area_strategy<Ring>());
         if (predicate(area, zero))
         {
             return visitor.template apply<no_failure>();
diff --git a/boost/geometry/algorithms/detail/overlay/add_rings.hpp b/boost/geometry/algorithms/detail/overlay/add_rings.hpp
index f64eb0b069..242e30cbcb 100644
--- a/boost/geometry/algorithms/detail/overlay/add_rings.hpp
+++ b/boost/geometry/algorithms/detail/overlay/add_rings.hpp
@@ -1,6 +1,12 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
+
+// This file was modified by Oracle on 2017.
+// Modifications copyright (c) 2017, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
 
 // Use, modification and distribution is subject to the Boost Software License,
 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
@@ -10,8 +16,10 @@
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_ADD_RINGS_HPP
 
 #include <boost/range.hpp>
+#include <boost/throw_exception.hpp>
 
 #include <boost/geometry/core/closure.hpp>
+#include <boost/geometry/core/exception.hpp>
 #include <boost/geometry/algorithms/area.hpp>
 #include <boost/geometry/algorithms/detail/overlay/convert_ring.hpp>
 #include <boost/geometry/algorithms/detail/overlay/get_ring.hpp>
@@ -87,7 +95,10 @@ inline OutputIterator add_rings(SelectionMap const& map,
             add_rings_error_handling error_handling = add_rings_ignore_unordered)
 {
     typedef typename SelectionMap::const_iterator iterator;
-    typedef typename AreaStrategy::return_type area_type;
+    typedef typename AreaStrategy::template result_type
+        <
+            GeometryOut
+        >::type area_type;
 
     area_type const zero = 0;
     std::size_t const min_num_points = core_detail::closure::minimum_ring_size
diff --git a/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp b/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp
deleted file mode 100644
index 3f2aea1b1d..0000000000
--- a/boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp
+++ /dev/null
@@ -1,256 +0,0 @@
-// Boost.Geometry (aka GGL, Generic Geometry Library)
-
-// Copyright (c) 2016 Barend Gehrels, Amsterdam, the Netherlands.
-
-// Use, modification and distribution is subject to the Boost Software License,
-// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
-// http://www.boost.org/LICENSE_1_0.txt)
-
-#ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP
-
-#include <set>
-
-#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
-
-namespace boost { namespace geometry
-{
-
-#ifndef DOXYGEN_NO_DETAIL
-namespace detail { namespace overlay { namespace sort_by_side
-{
-
-struct ring_with_direction
-{
-    ring_identifier ring_id;
-    direction_type direction;
-
-    signed_size_type turn_index;
-    int operation_index;
-    operation_type operation;
-    signed_size_type region_id;
-    bool isolated;
-
-    inline bool operator<(ring_with_direction const& other) const
-    {
-        return this->ring_id != other.ring_id
-                ? this->ring_id < other.ring_id
-                : this->direction < other.direction;
-    }
-
-    ring_with_direction()
-        : direction(dir_unknown)
-        , turn_index(-1)
-        , operation_index(0)
-        , operation(operation_none)
-        , region_id(-1)
-        , isolated(false)
-    {}
-};
-
-struct rank_with_rings
-{
-    // Define a set having a ring, with its direction (from/to). Each ring
-    // arrive at / leaves a cluster only once. TODO: this is not true for
-    // invalid ring. The rank needs to be considered too.
-    typedef std::set<ring_with_direction> container_type;
-    std::size_t rank;
-    container_type rings;
-
-    rank_with_rings()
-        : rank(0)
-    {
-    }
-
-    inline bool all_equal(direction_type dir_type) const
-    {
-        for (container_type::const_iterator it = rings.begin();
-             it != rings.end(); ++it)
-        {
-            if (it->direction != dir_type)
-            {
-                return false;
-            }
-        }
-        return true;
-    }
-
-    inline bool all_to() const
-    {
-        return all_equal(sort_by_side::dir_to);
-    }
-
-    inline bool all_from() const
-    {
-        return all_equal(sort_by_side::dir_from);
-    }
-
-    inline bool has_only(operation_type op) const
-    {
-        for (container_type::const_iterator it = rings.begin();
-             it != rings.end(); ++it)
-        {
-            const ring_with_direction& rwd = *it;
-            if (rwd.operation != op)
-            {
-                return false;
-            }
-        }
-        return true;
-    }
-
-    //! Check if set has both op1 and op2, but no others
-    inline bool has_only_both(operation_type op1, operation_type op2) const
-    {
-        bool has1 = false;
-        bool has2 = false;
-        for (container_type::const_iterator it = rings.begin();
-             it != rings.end(); ++it)
-        {
-            const ring_with_direction& rwd = *it;
-
-            if (rwd.operation == op1) { has1 = true; }
-            else if (rwd.operation == op2) { has2 = true; }
-            else { return false; }
-        }
-        return has1 && has2;
-    }
-
-    inline bool is_isolated() const
-    {
-        for (container_type::const_iterator it = rings.begin();
-             it != rings.end(); ++it)
-        {
-            const ring_with_direction& rwd = *it;
-            if (! rwd.isolated)
-            {
-                return false;
-            }
-        }
-        return true;
-    }
-
-    inline bool has_unique_region_id() const
-    {
-        signed_size_type region_id = -1;
-        for (container_type::const_iterator it = rings.begin();
-             it != rings.end(); ++it)
-        {
-            const ring_with_direction& rwd = *it;
-            if (region_id == -1)
-            {
-                region_id = rwd.region_id;
-            }
-            else if (rwd.region_id != region_id)
-            {
-                return false;
-            }
-        }
-        return true;
-    }
-
-    inline signed_size_type region_id() const
-    {
-        signed_size_type region_id = -1;
-        for (container_type::const_iterator it = rings.begin();
-             it != rings.end(); ++it)
-        {
-            const ring_with_direction& rwd = *it;
-            if (region_id == -1)
-            {
-                region_id = rwd.region_id;
-            }
-            else if (rwd.region_id != region_id)
-            {
-                return -1;
-            }
-        }
-        return region_id;
-    }
-
-    template <typename Turns>
-    inline bool traversable(Turns const& turns) const
-    {
-        typedef typename boost::range_value<Turns>::type turn_type;
-        typedef typename turn_type::turn_operation_type turn_operation_type;
-
-        for (container_type::const_iterator it = rings.begin();
-             it != rings.end(); ++it)
-        {
-            const ring_with_direction& rwd = *it;
-            turn_type const& turn = turns[rwd.turn_index];
-            turn_operation_type const& op = turn.operations[rwd.operation_index];
-
-            // TODO: this is still necessary, but makes it order-dependent
-            // which should not be done.
-
-            // This would obsolete the whole function and should be solved
-            // in a different way
-            if (op.visited.finalized() || op.visited.visited())
-            {
-                return false;
-            }
-        }
-        return true;
-    }
-
-};
-
-template <typename Sbs, typename Turns>
-inline void aggregate_operations(Sbs const& sbs, std::vector<rank_with_rings>& aggregation,
-                                 Turns const& turns,
-                                 operation_type target_operation)
-{
-    typedef typename boost::range_value<Turns>::type turn_type;
-    typedef typename turn_type::turn_operation_type turn_operation_type;
-
-    aggregation.clear();
-    for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
-    {
-        typename Sbs::rp const& ranked_point = sbs.m_ranked_points[i];
-
-        turn_type const& turn = turns[ranked_point.turn_index];
-
-        turn_operation_type const& op = turn.operations[ranked_point.operation_index];
-
-        if (! ((target_operation == operation_union && ranked_point.rank == 0)
-               || op.operation == target_operation
-               || op.operation == operation_continue
-               || (op.operation == operation_blocked && ranked_point.direction == dir_from)))
-        {
-            // Always take rank 0 (because self-turns are blocked)
-            // Don't consider union/blocked (aggregate is only used for intersections)
-            // Blocked is allowed for from
-            continue;
-        }
-
-        if (aggregation.empty() || aggregation.back().rank != ranked_point.rank)
-        {
-            rank_with_rings current;
-            current.rank = ranked_point.rank;
-            aggregation.push_back(current);
-        }
-
-        ring_with_direction rwd;
-        segment_identifier const& sid = ranked_point.seg_id;
-
-        rwd.ring_id = ring_identifier(sid.source_index, sid.multi_index, sid.ring_index);
-        rwd.direction = ranked_point.direction;
-        rwd.turn_index = ranked_point.turn_index;
-        rwd.operation_index = ranked_point.operation_index;
-        rwd.operation = op.operation;
-        rwd.region_id = op.enriched.region_id;
-        rwd.isolated = op.enriched.isolated;
-
-        aggregation.back().rings.insert(rwd);
-    }
-}
-
-
-}}} // namespace detail::overlay::sort_by_side
-#endif //DOXYGEN_NO_DETAIL
-
-
-}} // namespace boost::geometry
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_AGGREGATE_OPERATIONS_HPP
diff --git a/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp b/boost/geometry/algorithms/detail/overlay/append_no_dups_or_spikes.hpp
index fb73840798..724996ae33 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
@@ -102,6 +102,43 @@ inline void append_no_dups_or_spikes(Range& range, Point const& point,
     }
 }
 
+template <typename Range, typename Point, typename SideStrategy, typename RobustPolicy>
+inline void append_no_collinear(Range& range, Point const& point,
+        SideStrategy const& strategy,
+        RobustPolicy const& robust_policy)
+{
+    // Stricter version, not allowing any point in a linear row
+    // (spike, continuation or same point)
+
+    // 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
+        && points_equal_or_close(*(boost::begin(range)), point, robust_policy))
+    {
+        return;
+    }
+
+    traits::push_back<Range>::apply(range, point);
+
+    // If a point is equal, or forming a spike, remove the pen-ultimate point
+    // because this one caused the spike.
+    // If so, the now-new-pen-ultimate point can again cause a spike
+    // (possibly at a corner). So keep doing this.
+    // Besides spikes it will also avoid adding duplicates.
+    while(boost::size(range) >= 3
+            && point_is_collinear(point,
+                *(boost::end(range) - 3),
+                *(boost::end(range) - 2),
+                strategy,
+                robust_policy))
+    {
+        // Use the Concept/traits, so resize and append again
+        traits::resize<Range>::apply(range, boost::size(range) - 2);
+        traits::push_back<Range>::apply(range, point);
+    }
+}
+
 template <typename Range, typename SideStrategy, typename RobustPolicy>
 inline void clean_closing_dups_and_spikes(Range& range,
                 SideStrategy const& strategy,
@@ -137,8 +174,8 @@ inline void clean_closing_dups_and_spikes(Range& range,
         }
 
         // Check if closing point is a spike (this is so if the second point is
-        // considered as a spike w.r.t. the last segment)
-        if (point_is_spike_or_equal(*second, *ultimate, *first, strategy, robust_policy))
+        // considered as collinear w.r.t. the last segment)
+        if (point_is_collinear(*second, *ultimate, *first, strategy, robust_policy))
         {
             range::erase(range, first);
             if (BOOST_GEOMETRY_CONDITION(closed))
diff --git a/boost/geometry/algorithms/detail/overlay/assign_parents.hpp b/boost/geometry/algorithms/detail/overlay/assign_parents.hpp
index c8ce651007..3be5393486 100644
--- a/boost/geometry/algorithms/detail/overlay/assign_parents.hpp
+++ b/boost/geometry/algorithms/detail/overlay/assign_parents.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -210,10 +211,9 @@ struct assign_visitor
 };
 
 
-
-
 template
 <
+    overlay_type OverlayType,
     typename Geometry1, typename Geometry2,
     typename RingCollection,
     typename RingMap,
@@ -223,10 +223,13 @@ inline void assign_parents(Geometry1 const& geometry1,
             Geometry2 const& geometry2,
             RingCollection const& collection,
             RingMap& ring_map,
-            Strategy const& strategy,
-            bool check_for_orientation = false,
-            bool discard_double_negative = false)
+            Strategy const& strategy)
 {
+    static bool const is_difference = OverlayType == overlay_difference;
+    static bool const is_buffer = OverlayType == overlay_buffer;
+    static bool const is_dissolve = OverlayType == overlay_dissolve;
+    static bool const check_for_orientation = is_buffer || is_dissolve;
+
     typedef typename geometry::tag<Geometry1>::type tag1;
     typedef typename geometry::tag<Geometry2>::type tag2;
 
@@ -236,7 +239,7 @@ inline void assign_parents(Geometry1 const& geometry1,
     typedef typename Strategy::template area_strategy
         <
             point_type
-        >::type::return_type area_result_type;
+        >::type::template result_type<point_type>::type area_result_type;
 
     typedef typename RingMap::iterator map_iterator_type;
 
@@ -293,12 +296,14 @@ inline void assign_parents(Geometry1 const& geometry1,
                 return;
             }
 
-            if (count_positive == 1)
+            if (count_positive == 1 && ! is_difference && ! is_dissolve)
             {
                 // Optimization for one outer ring
                 // -> assign this as parent to all others (all interior rings)
                 // In unions, this is probably the most occuring case and gives
                 //    a dramatic improvement (factor 5 for star_comb testcase)
+                // In difference or other cases where interior rings might be
+                // located outside the outer ring, this cannot be done
                 ring_identifier id_of_positive = vector[index_positive].id;
                 ring_info_type& outer = ring_map[id_of_positive];
                 index = 0;
@@ -346,13 +351,13 @@ inline void assign_parents(Geometry1 const& geometry1,
                 bool const pos = math::larger(info.get_area(), 0);
                 bool const parent_pos = math::larger(parent.area, 0);
 
-                bool const double_neg = discard_double_negative && ! pos && ! parent_pos;
+                bool const double_neg = is_dissolve && ! pos && ! parent_pos;
 
                 if ((pos && parent_pos) || double_neg)
                 {
                     // Discard positive inner ring with positive parent
                     // Also, for some cases (dissolve), negative inner ring
-                    // with negative parent shouild be discarded
+                    // with negative parent should be discarded
                     info.discarded = true;
                 }
 
@@ -386,6 +391,7 @@ inline void assign_parents(Geometry1 const& geometry1,
 // Version for one geometry (called by buffer/dissolve)
 template
 <
+    overlay_type OverlayType,
     typename Geometry,
     typename RingCollection,
     typename RingMap,
@@ -394,16 +400,13 @@ template
 inline void assign_parents(Geometry const& geometry,
             RingCollection const& collection,
             RingMap& ring_map,
-            Strategy const& strategy,
-            bool check_for_orientation,
-            bool discard_double_negative)
+            Strategy const& strategy)
 {
     // Call it with an empty geometry as second geometry (source_id == 1)
     // (ring_map should be empty for source_id==1)
-
     Geometry empty;
-    assign_parents(geometry, empty, collection, ring_map, strategy,
-                   check_for_orientation, discard_double_negative);
+    assign_parents<OverlayType>(geometry, empty,
+            collection, ring_map, strategy);
 }
 
 
diff --git a/boost/geometry/algorithms/detail/overlay/cluster_info.hpp b/boost/geometry/algorithms/detail/overlay/cluster_info.hpp
index 5b460919f1..19343488f3 100644
--- a/boost/geometry/algorithms/detail/overlay/cluster_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/cluster_info.hpp
@@ -26,14 +26,11 @@ struct cluster_info
 {
     std::set<signed_size_type> turn_indices;
 
-    bool switch_source; // For clusters with a touch, conform turn_info uu
-
     //! Number of open spaces (e.g. 2 for touch)
     std::size_t open_count;
 
     inline cluster_info()
-        : switch_source(false)
-        , open_count(0)
+        : open_count(0)
     {}
 };
 
diff --git a/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp b/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp
index e25445651a..a35be052a0 100644
--- a/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp
+++ b/boost/geometry/algorithms/detail/overlay/enrich_intersection_points.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2007-2012 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -51,6 +52,21 @@ namespace boost { namespace geometry
 namespace detail { namespace overlay
 {
 
+template <typename Turns>
+struct discarded_turn
+{
+    discarded_turn(Turns const& turns)
+        : m_turns(turns)
+    {}
+
+    template <typename IndexedTurn>
+    inline bool operator()(IndexedTurn const& indexed) const
+    {
+        return m_turns[indexed.turn_index].discarded;
+    }
+
+    Turns const& m_turns;
+};
 
 // Sorts IP-s of this ring on segment-identifier, and if on same segment,
 //  on distance.
@@ -68,7 +84,6 @@ template
 >
 inline void enrich_sort(Operations& operations,
             Turns const& turns,
-            operation_type for_operation,
             Geometry1 const& geometry1,
             Geometry2 const& geometry2,
             RobustPolicy const& robust_policy,
@@ -84,12 +99,13 @@ inline void enrich_sort(Operations& operations,
                     RobustPolicy,
                     SideStrategy,
                     Reverse1, Reverse2
-                >(turns, for_operation, geometry1, geometry2, robust_policy, strategy));
+                >(turns, geometry1, geometry2, robust_policy, strategy));
 }
 
 
 template <typename Operations, typename Turns>
-inline void enrich_assign(Operations& operations, Turns& turns)
+inline void enrich_assign(Operations& operations, Turns& turns,
+                          bool check_turns)
 {
     typedef typename boost::range_value<Turns>::type turn_type;
     typedef typename turn_type::turn_operation_type op_type;
@@ -110,14 +126,17 @@ inline void enrich_assign(Operations& operations, Turns& turns)
             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)
+            if (check_turns && it->turn_index == next->turn_index)
             {
+                // Normal behaviour: next points at next turn, increase next.
+                // For dissolve this should not be done, turn_index is often
+                // the same for two consecutive operations
                 ++next;
             }
 
             // Cluster behaviour: next should point after cluster, unless
             // their seg_ids are not the same
+            // (For dissolve, this is still to be examined - TODO)
             while (turn.is_clustered()
                    && it->turn_index != next->turn_index
                    && turn.cluster_id == turns[next->turn_index].cluster_id
@@ -142,6 +161,11 @@ inline void enrich_assign(Operations& operations, Turns& turns)
                 // (this is one not circular therefore fraction is considered)
                 op.enriched.next_ip_index = static_cast<signed_size_type>(next->turn_index);
             }
+
+            if (! check_turns)
+            {
+                ++next;
+            }
         }
     }
 
@@ -176,6 +200,82 @@ inline void enrich_assign(Operations& operations, Turns& turns)
 
 }
 
+template <typename Operations, typename Turns>
+inline void enrich_adapt(Operations& operations, Turns& turns)
+{
+    typedef typename boost::range_value<Turns>::type turn_type;
+    typedef typename turn_type::turn_operation_type op_type;
+    typedef typename boost::range_value<Operations>::type indexed_turn_type;
+
+    if (operations.size() < 3)
+    {
+        // If it is empty, or contains one or two turns, it makes no sense
+        return;
+    }
+
+    // Operations is a vector of indexed_turn_operation<>
+
+    // Last index:
+    std::size_t const x = operations.size() - 1;
+    bool next_phase = false;
+
+    for (std::size_t i = 0; i < operations.size(); i++)
+    {
+        indexed_turn_type const& indexed = operations[i];
+
+        turn_type& turn = turns[indexed.turn_index];
+        op_type& op = turn.operations[indexed.operation_index];
+
+        // Previous/next index
+        std::size_t const p = i > 0 ? i - 1 : x;
+        std::size_t const n = i < x ? i + 1 : 0;
+
+        turn_type const& next_turn = turns[operations[n].turn_index];
+        op_type const& next_op = next_turn.operations[operations[n].operation_index];
+
+        if (op.seg_id.segment_index == next_op.seg_id.segment_index)
+        {
+            turn_type const& prev_turn = turns[operations[p].turn_index];
+            op_type const& prev_op = prev_turn.operations[operations[p].operation_index];
+            if (op.seg_id.segment_index == prev_op.seg_id.segment_index)
+            {
+                op.enriched.startable = false;
+                next_phase = true;
+            }
+        }
+    }
+
+    if (! next_phase)
+    {
+        return;
+    }
+
+    // Discard turns which are both non-startable
+    next_phase = false;
+    for (typename boost::range_iterator<Turns>::type
+            it = boost::begin(turns);
+         it != boost::end(turns);
+         ++it)
+    {
+        turn_type& turn = *it;
+        if (! turn.operations[0].enriched.startable
+            && ! turn.operations[1].enriched.startable)
+        {
+            turn.discarded = true;
+            next_phase = true;
+        }
+    }
+
+    if (! next_phase)
+    {
+        return;
+    }
+
+    // Remove discarded turns from operations to avoid having them as next turn
+    discarded_turn<Turns> const predicate(turns);
+    operations.erase(std::remove_if(boost::begin(operations),
+        boost::end(operations), predicate), boost::end(operations));
+}
 
 template <typename Turns, typename MappedVector>
 inline void create_map(Turns const& turns, MappedVector& mapped_vector)
@@ -255,8 +355,8 @@ inline void calculate_remaining_distance(Turns& turns)
             continue;
         }
 
-        int const to_index0 = op0.enriched.get_next_turn_index();
-        int const to_index1 = op1.enriched.get_next_turn_index();
+        signed_size_type const to_index0 = op0.enriched.get_next_turn_index();
+        signed_size_type const to_index1 = op1.enriched.get_next_turn_index();
         if (to_index0 >= 0
                 && to_index1 >= 0
                 && to_index0 != to_index1)
@@ -311,6 +411,7 @@ inline void enrich_intersection_points(Turns& turns,
             = target_operation == detail::overlay::operation_union
             ? detail::overlay::operation_intersection
             : detail::overlay::operation_union;
+    static const bool is_dissolve = OverlayType == overlay_dissolve;
 
     typedef typename boost::range_value<Turns>::type turn_type;
     typedef typename turn_type::turn_operation_type op_type;
@@ -338,31 +439,25 @@ inline void enrich_intersection_points(Turns& turns,
     {
         turn_type& turn = *it;
 
-        if (turn.both(detail::overlay::operation_none))
-        {
-            turn.discarded = true;
-            continue;
-        }
-
-        if (turn.both(opposite_operation))
+        if (turn.both(detail::overlay::operation_none)
+            || turn.both(opposite_operation)
+            || (detail::overlay::is_self_turn<OverlayType>(turn)
+                && ! turn.is_clustered()
+                && ! turn.both(target_operation)))
         {
             // For intersections, remove uu to avoid the need to travel
             // a union (during intersection) in uu/cc clusters (e.g. #31,#32,#33)
-            // Also, for union, discard ii
+
+            // Similarly, for union, discard ii
+
+            // Only keep self-uu-turns or self-ii-turns
+
+            // Blocked (or combination with blocked is still needed for difference)
             turn.discarded = true;
             turn.cluster_id = -1;
             continue;
         }
 
-        if (detail::overlay::is_self_turn<OverlayType>(turn)
-            && ! turn.is_clustered()
-            && ! turn.both(target_operation))
-        {
-            // Only keep self-uu-turns or self-ii-turns
-           turn.discarded = true;
-           continue;
-        }
-
         if (! turn.discarded
             && turn.both(detail::overlay::operation_continue))
         {
@@ -370,16 +465,19 @@ inline void enrich_intersection_points(Turns& turns,
         }
     }
 
-    detail::overlay::discard_closed_turns
-        <
-            OverlayType,
-            target_operation
-        >::apply(turns, clusters, geometry1, geometry2);
-    detail::overlay::discard_open_turns
-        <
-            OverlayType,
-            target_operation
-        >::apply(turns, clusters, geometry1, geometry2);
+    if (! is_dissolve)
+    {
+        detail::overlay::discard_closed_turns
+            <
+                OverlayType,
+                target_operation
+            >::apply(turns, clusters, geometry1, geometry2);
+        detail::overlay::discard_open_turns
+            <
+                OverlayType,
+                target_operation
+            >::apply(turns, clusters, geometry1, geometry2);
+    }
 
     // Create a map of vectors of indexed operation-types to be able
     // to sort intersection points PER RING
@@ -399,7 +497,7 @@ inline void enrich_intersection_points(Turns& turns,
         << mit->first << std::endl;
 #endif
         detail::overlay::enrich_sort<Reverse1, Reverse2>(
-                    mit->second, turns, target_operation,
+                    mit->second, turns,
                     geometry1, geometry2,
                     robust_policy, strategy);
     }
@@ -413,11 +511,13 @@ inline void enrich_intersection_points(Turns& turns,
     std::cout << "ENRICH-assign Ring "
         << mit->first << std::endl;
 #endif
-        detail::overlay::enrich_assign(mit->second, turns);
-    }
+        if (is_dissolve)
+        {
+            detail::overlay::enrich_adapt(mit->second, turns);
+        }
 
-    // Check some specific type of self-turns (after getting enriched info)
-    detail::overlay::discard_self_turns_which_loop<OverlayType>(turns);
+        detail::overlay::enrich_assign(mit->second, turns, ! is_dissolve);
+    }
 
     if (has_colocations)
     {
diff --git a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
index fdffd665e4..e01c13f749 100644
--- a/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/enrichment_info.hpp
@@ -35,6 +35,7 @@ struct enrichment_info
         , travels_to_ip_index(-1)
         , next_ip_index(-1)
         , startable(true)
+        , prefer_start(true)
         , count_left(0)
         , count_right(0)
         , rank(-1)
@@ -60,6 +61,7 @@ struct enrichment_info
     signed_size_type next_ip_index;
 
     bool startable; // Can be used to start in traverse
+    bool prefer_start; // Is preferred as starting point (if true)
 
     // Counts if polygons left/right of this operation
     std::size_t count_left;
diff --git a/boost/geometry/algorithms/detail/overlay/follow.hpp b/boost/geometry/algorithms/detail/overlay/follow.hpp
index 4a5993ea31..d948c4f670 100644
--- a/boost/geometry/algorithms/detail/overlay/follow.hpp
+++ b/boost/geometry/algorithms/detail/overlay/follow.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2007-2014 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2014, 2017.
 // Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
@@ -59,7 +60,7 @@ template
     typename Polygon,
     typename PtInPolyStrategy
 >
-static inline bool last_covered_by(Turn const& turn, Operation const& op,
+static inline bool last_covered_by(Turn const& /*turn*/, Operation const& op,
                 LineString const& linestring, Polygon const& polygon,
                 PtInPolyStrategy const& strategy)
 {
diff --git a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
index c634fc450f..6bb30fcce5 100644
--- a/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
+++ b/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -30,6 +31,7 @@
 #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>
+#include <boost/geometry/util/condition.hpp>
 
 #if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS)
 #  include <iostream>
@@ -312,17 +314,17 @@ inline void assign_cluster_to_turns(Turns& turns,
         {
             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())
+            typename ClusterPerSegment::const_iterator cit = cluster_per_segment.find(seg_frac);
+            if (cit != cluster_per_segment.end())
             {
 #if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS)
                 if (turn.is_clustered()
-                        && turn.cluster_id != it->second)
+                        && turn.cluster_id != cit->second)
                 {
                     std::cout << " CONFLICT " << std::endl;
                 }
 #endif
-                turn.cluster_id = it->second;
+                turn.cluster_id = cit->second;
                 clusters[turn.cluster_id].turn_indices.insert(turn_index);
             }
         }
@@ -439,7 +441,6 @@ inline bool is_ie_turn(segment_identifier const& ext_seg_0,
 template
 <
     bool Reverse0, bool Reverse1, // Reverse interpretation interior/exterior
-    overlay_type OverlayType,
     typename Turns,
     typename Clusters
 >
@@ -552,7 +553,7 @@ template
     typename Clusters
 >
 inline void check_colocation(bool& has_blocked,
-        int cluster_id, Turns const& turns, Clusters const& clusters)
+        signed_size_type cluster_id, Turns const& turns, Clusters const& clusters)
 {
     typedef typename boost::range_value<Turns>::type turn_type;
 
@@ -598,6 +599,8 @@ template
 inline bool handle_colocations(Turns& turns, Clusters& clusters,
         Geometry1 const& geometry1, Geometry2 const& geometry2)
 {
+    static const detail::overlay::operation_type target_operation
+            = detail::overlay::operation_from_overlay<OverlayType>::value;
     typedef std::map
         <
             segment_identifier,
@@ -610,7 +613,7 @@ inline bool handle_colocations(Turns& turns, Clusters& clusters,
     // that information can be used for the interior ring too
     map_type map;
 
-    int index = 0;
+    signed_size_type index = 0;
     for (typename boost::range_iterator<Turns>::type
             it = boost::begin(turns);
          it != boost::end(turns);
@@ -677,12 +680,15 @@ inline bool handle_colocations(Turns& turns, Clusters& clusters,
     // Get colocated information here and not later, to keep information
     // on turns which are discarded afterwards
     set_colocation<OverlayType>(turns, clusters);
-    discard_interior_exterior_turns
-        <
-            do_reverse<geometry::point_order<Geometry1>::value>::value != Reverse1,
-            do_reverse<geometry::point_order<Geometry2>::value>::value != Reverse2,
-            OverlayType
-        >(turns, clusters);
+
+    if (BOOST_GEOMETRY_CONDITION(target_operation == operation_intersection))
+    {
+        discard_interior_exterior_turns
+            <
+                do_reverse<geometry::point_order<Geometry1>::value>::value != Reverse1,
+                do_reverse<geometry::point_order<Geometry2>::value>::value != Reverse2
+            >(turns, clusters);
+    }
 
 #if defined(BOOST_GEOMETRY_DEBUG_HANDLE_COLOCATIONS)
     std::cout << "*** Colocations " << map.size() << std::endl;
@@ -794,9 +800,7 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns,
 
         cinfo.open_count = sbs.open_count(for_operation);
 
-        bool const set_startable
-                = OverlayType != overlay_dissolve_union
-                && OverlayType != overlay_dissolve_intersection;
+        bool const set_startable = OverlayType != overlay_dissolve;
 
         // Unset the startable flag for all 'closed' zones. This does not
         // apply for self-turns, because those counts are not from both
diff --git a/boost/geometry/algorithms/detail/overlay/handle_self_turns.hpp b/boost/geometry/algorithms/detail/overlay/handle_self_turns.hpp
index 9c4a3094e0..5ec2a10cf0 100644
--- a/boost/geometry/algorithms/detail/overlay/handle_self_turns.hpp
+++ b/boost/geometry/algorithms/detail/overlay/handle_self_turns.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2017 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // Use, modification and distribution is subject to the Boost Software License,
 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
@@ -14,6 +15,7 @@
 #include <boost/geometry/algorithms/detail/overlay/cluster_info.hpp>
 #include <boost/geometry/algorithms/detail/overlay/is_self_turn.hpp>
 #include <boost/geometry/algorithms/detail/overlay/overlay_type.hpp>
+#include <boost/geometry/algorithms/covered_by.hpp>
 #include <boost/geometry/algorithms/within.hpp>
 
 namespace boost { namespace geometry
@@ -23,6 +25,38 @@ namespace boost { namespace geometry
 namespace detail { namespace overlay
 {
 
+template <overlay_type OverlayType>
+struct check_within
+{
+    template <typename Turn, typename Geometry0, typename Geometry1>
+    static inline
+    bool apply(Turn const& turn, Geometry0 const& geometry0,
+                Geometry1 const& geometry1)
+    {
+        // Operations 0 and 1 have the same source index in self-turns
+        return turn.operations[0].seg_id.source_index == 0
+            ? geometry::within(turn.point, geometry1)
+            : geometry::within(turn.point, geometry0);
+    }
+
+};
+
+template <>
+struct check_within<overlay_difference>
+{
+    template <typename Turn, typename Geometry0, typename Geometry1>
+    static inline
+    bool apply(Turn const& turn, Geometry0 const& geometry0,
+                Geometry1 const& geometry1)
+    {
+        // difference = intersection(a, reverse(b))
+        // therefore we should reverse the meaning of within for geometry1
+        return turn.operations[0].seg_id.source_index == 0
+            ? ! geometry::covered_by(turn.point, geometry1)
+            : geometry::within(turn.point, geometry0);
+    }
+};
+
 struct discard_turns
 {
     template <typename Turns, typename Clusters, typename Geometry0, typename Geometry1>
@@ -41,7 +75,7 @@ struct discard_closed_turns<overlay_union, operation_union>
 
     template <typename Turns, typename Clusters, typename Geometry0, typename Geometry1>
     static inline
-    void apply(Turns& turns, Clusters const& clusters,
+    void apply(Turns& turns, Clusters const& /*clusters*/,
             Geometry0 const& geometry0, Geometry1 const& geometry1)
     {
         typedef typename boost::range_value<Turns>::type turn_type;
@@ -53,55 +87,25 @@ struct discard_closed_turns<overlay_union, operation_union>
         {
             turn_type& turn = *it;
 
-            if (turn.discarded || ! is_self_turn<overlay_union>(turn))
-            {
-                continue;
-            }
-
-            bool const within =
-                    turn.operations[0].seg_id.source_index == 0
-                    ? geometry::within(turn.point, geometry1)
-                    : geometry::within(turn.point, geometry0);
-
-            if (within)
+            if (! turn.discarded
+                && is_self_turn<overlay_union>(turn)
+                && check_within<overlay_union>::apply(turn,
+                        geometry0, geometry1))
             {
-                // It is in the interior of the other geometry
+                // Turn is in the interior of other geometry
                 turn.discarded = true;
             }
         }
     }
 };
 
+template <overlay_type OverlayType>
 struct discard_self_intersection_turns
 {
 private :
 
     template <typename Turns, typename Clusters>
     static inline
-    bool any_blocked(signed_size_type cluster_id,
-            const Turns& turns, Clusters const& clusters)
-    {
-        typename Clusters::const_iterator cit = clusters.find(cluster_id);
-        if (cit == clusters.end())
-        {
-            return false;
-        }
-        cluster_info const& cinfo = cit->second;
-        for (std::set<signed_size_type>::const_iterator it
-             = cinfo.turn_indices.begin();
-             it != cinfo.turn_indices.end(); ++it)
-        {
-            typename boost::range_value<Turns>::type const& turn = turns[*it];
-            if (turn.any_blocked())
-            {
-                return true;
-            }
-        }
-        return false;
-    }
-
-    template <typename Turns, typename Clusters>
-    static inline
     bool is_self_cluster(signed_size_type cluster_id,
             const Turns& turns, Clusters const& clusters)
     {
@@ -116,7 +120,7 @@ private :
              = cinfo.turn_indices.begin();
              it != cinfo.turn_indices.end(); ++it)
         {
-            if (! is_self_turn<overlay_intersection>(turns[*it]))
+            if (! is_self_turn<OverlayType>(turns[*it]))
             {
                 return false;
             }
@@ -125,16 +129,6 @@ private :
         return true;
     }
 
-    template <typename Turn, typename Geometry0, typename Geometry1>
-    static inline
-    bool within(Turn const& turn, Geometry0 const& geometry0,
-                Geometry1 const& geometry1)
-    {
-        return turn.operations[0].seg_id.source_index == 0
-            ? geometry::within(turn.point, geometry1)
-            : geometry::within(turn.point, geometry0);
-    }
-
     template <typename Turns, typename Clusters,
               typename Geometry0, typename Geometry1>
     static inline
@@ -144,17 +138,21 @@ private :
         for (typename Clusters::const_iterator cit = clusters.begin();
              cit != clusters.end(); ++cit)
         {
-            signed_size_type cluster_id = cit->first;
+            signed_size_type const cluster_id = cit->first;
 
             // If there are only self-turns in the cluster, the cluster should
             // be located within the other geometry, for intersection
-            if (is_self_cluster(cluster_id, turns, clusters))
+            if (! cit->second.turn_indices.empty()
+                && is_self_cluster(cluster_id, turns, clusters))
             {
                 cluster_info const& cinfo = cit->second;
-                if (! within(turns[*cinfo.turn_indices.begin()], geometry0, geometry1))
+                signed_size_type const index = *cinfo.turn_indices.begin();
+                if (! check_within<OverlayType>::apply(turns[index],
+                                                       geometry0, geometry1))
                 {
                     // Discard all turns in cluster
-                    for (std::set<signed_size_type>::const_iterator sit = cinfo.turn_indices.begin();
+                    for (std::set<signed_size_type>::const_iterator sit
+                         = cinfo.turn_indices.begin();
                          sit != cinfo.turn_indices.end(); ++sit)
                     {
                         turns[*sit].discarded = true;
@@ -183,109 +181,34 @@ public :
         {
             turn_type& turn = *it;
 
-            if (turn.discarded || ! is_self_turn<overlay_intersection>(turn))
-            {
-                continue;
-            }
-
-            segment_identifier const& id0 = turn.operations[0].seg_id;
-            segment_identifier const& id1 = turn.operations[1].seg_id;
-            if (id0.multi_index != id1.multi_index
-                    || (id0.ring_index == -1 && id1.ring_index == -1)
-                    || (id0.ring_index >= 0 && id1.ring_index >= 0))
-            {
-                // Not an ii ring (int/ext) on same ring
-                continue;
-            }
-
-            if (turn.is_clustered() && turn.has_colocated_both)
-            {
-                // Don't delete a self-ii-turn colocated with another ii-turn
-                // (for example #case_recursive_boxes_70)
-                // But for some cases (#case_58_iet) they should be deleted,
-                // there are many self-turns there and also blocked turns there
-                if (! any_blocked(turn.cluster_id, turns, clusters))
-                {
-                    continue;
-                }
-            }
-
             // It is a ii self-turn
             // Check if it is within the other geometry
-            // If not, it can be ignored
-            if (! within(turn, geometry0, geometry1))
+            if (! turn.discarded
+                && is_self_turn<overlay_intersection>(turn)
+                && ! check_within<OverlayType>::apply(turn, geometry0, geometry1))
             {
-                // It is not within another geometry, discard the turn
-                turn.discarded = true;
+                // It is not within another geometry, set it as non startable.
+                // It still might be traveled (#case_recursive_boxes_70)
+                turn.operations[0].enriched.startable = false;
+                turn.operations[1].enriched.startable = false;
             }
         }
     }
 };
 
+
 template <overlay_type OverlayType, operation_type OperationType>
 struct discard_open_turns : discard_turns {};
 
-// Handler it for intersection
+// Handler for intersection
 template <>
 struct discard_open_turns<overlay_intersection, operation_intersection>
-        : discard_self_intersection_turns {};
+        : discard_self_intersection_turns<overlay_intersection> {};
 
-// For difference, it should be done in a different way (TODO)
-
-
-template <overlay_type OverlayType, typename Turns>
-inline void discard_self_turns_which_loop(Turns& turns)
-{
-    if (operation_from_overlay<OverlayType>::value == operation_union)
-    {
-        // For union, self-turn i/u traveling to itself are allowed to form
-        // holes. #case_recursive_boxes_37
-        // TODO: this can be finetuned by inspecting the cluster too,
-        // and if there are non-self-turns the polygons on their sides can
-        // be checked
-        return;
-    }
-
-    typedef typename boost::range_value<Turns>::type turn_type;
-    typedef typename turn_type::turn_operation_type op_type;
-
-    signed_size_type turn_index = 0;
-    for (typename boost::range_iterator<Turns>::type
-            it = boost::begin(turns);
-         it != boost::end(turns);
-         ++it, ++turn_index)
-    {
-        turn_type& turn = *it;
-
-        if (! is_self_turn<OverlayType>(turn))
-        {
-            continue;
-        }
-        if (! turn.combination(operation_intersection, operation_union))
-        {
-            // ii may travel to itself
-            continue;
-        }
-
-        for (int i = 0; i < 2; i++)
-        {
-            op_type& op = turn.operations[i];
-
-            if (op.enriched.startable
-                && op.operation == operation_intersection
-                && op.enriched.get_next_turn_index() == turn_index)
-            {
-                // Self-turn i/u, i part traveling to itself. Discard it.
-                // (alternatively it might be made unstartable - but the
-                // intersection-operation may not be traveled anyway, and the
-                // union-operation is not traveled at all in intersections
-                // #case_recursive_boxes_77
-                turn.discarded = true;
-            }
-        }
-    }
-
-}
+// Handler for difference, with different meaning of 'within'
+template <>
+struct discard_open_turns<overlay_difference, operation_intersection>
+        : discard_self_intersection_turns<overlay_difference> {};
 
 }} // namespace detail::overlay
 #endif //DOXYGEN_NO_DETAIL
diff --git a/boost/geometry/algorithms/detail/overlay/is_self_turn.hpp b/boost/geometry/algorithms/detail/overlay/is_self_turn.hpp
index 9423a24b33..448c04404f 100644
--- a/boost/geometry/algorithms/detail/overlay/is_self_turn.hpp
+++ b/boost/geometry/algorithms/detail/overlay/is_self_turn.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2017-2017 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // Use, modification and distribution is subject to the Boost Software License,
 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
@@ -34,27 +35,17 @@ template <>
 struct is_self_turn_check<overlay_buffer>
 {
     template <typename Turn>
-    static inline bool apply(Turn const& turn)
+    static inline bool apply(Turn const& /*turn*/)
     {
         return false;
     }
 };
 
 template <>
-struct is_self_turn_check<overlay_dissolve_union>
+struct is_self_turn_check<overlay_dissolve>
 {
     template <typename Turn>
-    static inline bool apply(Turn const& turn)
-    {
-        return false;
-    }
-};
-
-template <>
-struct is_self_turn_check<overlay_dissolve_intersection>
-{
-    template <typename Turn>
-    static inline bool apply(Turn const& turn)
+    static inline bool apply(Turn const& /*turn*/)
     {
         return false;
     }
diff --git a/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp b/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp
index dd30635ee2..4b8752798a 100644
--- a/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp
+++ b/boost/geometry/algorithms/detail/overlay/less_by_segment_ratio.hpp
@@ -71,13 +71,11 @@ template
 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
             , SideStrategy const& strategy)
         : m_turns(turns)
-        , m_for_operation(for_operation)
         , m_geometry1(geometry1)
         , m_geometry2(geometry2)
         , m_robust_policy(robust_policy)
@@ -88,7 +86,6 @@ struct less_by_segment_ratio
 private :
 
     Turns const& m_turns;
-    operation_type m_for_operation;
     Geometry1 const& m_geometry1;
     Geometry2 const& m_geometry2;
     RobustPolicy const& m_robust_policy;
diff --git a/boost/geometry/algorithms/detail/overlay/overlay.hpp b/boost/geometry/algorithms/detail/overlay/overlay.hpp
index f24cde8b8f..5094c6c96c 100644
--- a/boost/geometry/algorithms/detail/overlay/overlay.hpp
+++ b/boost/geometry/algorithms/detail/overlay/overlay.hpp
@@ -1,7 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
-// Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland
+// Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland
 
 // This file was modified by Oracle on 2015, 2017.
 // Modifications copyright (c) 2015-2017, Oracle and/or its affiliates.
@@ -49,6 +49,7 @@
 
 #include <boost/geometry/policies/robustness/segment_ratio_type.hpp>
 
+#include <boost/geometry/util/condition.hpp>
 
 #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
 #  include <boost/geometry/io/dsv/write.hpp>
@@ -88,6 +89,10 @@ struct overlay_null_visitor
     template <typename Turns, typename Turn, typename Operation>
     void visit_traverse_reject(Turns const& , Turn const& , Operation const& , traverse_error_type )
     {}
+
+    template <typename Rings>
+    void visit_generated_rings(Rings const& )
+    {}
 };
 
 template
@@ -135,9 +140,9 @@ inline void get_ring_turn_info(TurnInfoMap& turn_info_map, Turns const& turns, C
 
             if (! is_self_turn<OverlayType>(turn)
                 && (
-                    (target_operation == operation_union
+                    (BOOST_GEOMETRY_CONDITION(target_operation == operation_union)
                       && op.enriched.count_left > 0)
-                  || (target_operation == operation_intersection
+                  || (BOOST_GEOMETRY_CONDITION(target_operation == operation_intersection)
                       && op.enriched.count_right <= 2)))
             {
                 // Avoid including untraversed rings which have polygons on
@@ -205,7 +210,7 @@ inline OutputIterator return_if_one_input_is_empty(Geometry1 const& geometry1,
             typename Strategy::template area_strategy
                 <
                     point_type1
-                >::type::return_type
+                >::type::template result_type<point_type1>::type
         > properties;
 
 // Silence warning C4127: conditional expression is constant
@@ -233,7 +238,7 @@ inline OutputIterator return_if_one_input_is_empty(Geometry1 const& geometry1,
 
     select_rings<OverlayType>(geometry1, geometry2, empty, all_of_one_of_them, strategy);
     ring_container_type rings;
-    assign_parents(geometry1, geometry2, rings, all_of_one_of_them, strategy);
+    assign_parents<OverlayType>(geometry1, geometry2, rings, all_of_one_of_them, strategy);
     return add_rings<GeometryOut>(all_of_one_of_them, geometry1, geometry2, rings, out,
                                   strategy.template get_area_strategy<point_type1>());
 }
@@ -306,7 +311,7 @@ std::cout << "get turns" << std::endl;
 
         visitor.visit_turns(1, turns);
 
-#ifdef BOOST_GEOMETRY_INCLUDE_SELF_TURNS
+#if ! defined(BOOST_GEOMETRY_NO_SELF_TURNS)
         if (needs_self_turns<Geometry1>::apply(geometry1))
         {
             self_get_turn_points::self_turns<Reverse1, assign_null_policy>(geometry1,
@@ -362,7 +367,7 @@ std::cout << "traverse" << std::endl;
         typedef ring_properties
             <
                 point_type,
-                typename area_strategy_type::return_type
+                typename area_strategy_type::template result_type<point_type>::type
             > properties;
 
         // Select all rings which are NOT touched by any intersection point
@@ -385,7 +390,8 @@ std::cout << "traverse" << std::endl;
             }
         }
 
-        assign_parents(geometry1, geometry2, rings, selected_ring_properties, strategy);
+        assign_parents<OverlayType>(geometry1, geometry2,
+            rings, selected_ring_properties, strategy);
 
         // NOTE: There is no need to check result area for union because
         // as long as the polygons in the input are valid the resulting
diff --git a/boost/geometry/algorithms/detail/overlay/overlay_type.hpp b/boost/geometry/algorithms/detail/overlay/overlay_type.hpp
index f3ec9eaa64..c76d440f91 100644
--- a/boost/geometry/algorithms/detail/overlay/overlay_type.hpp
+++ b/boost/geometry/algorithms/detail/overlay/overlay_type.hpp
@@ -21,8 +21,7 @@ enum overlay_type
     overlay_intersection,
     overlay_difference,
     overlay_buffer,
-    overlay_dissolve_union,
-    overlay_dissolve_intersection
+    overlay_dissolve
 };
 
 #ifndef DOXYGEN_NO_DETAIL
@@ -70,17 +69,11 @@ struct operation_from_overlay<overlay_difference>
 };
 
 template <>
-struct operation_from_overlay<overlay_dissolve_union>
+struct operation_from_overlay<overlay_dissolve>
 {
     static const operation_type value = operation_union;
 };
 
-template <>
-struct operation_from_overlay<overlay_dissolve_intersection>
-{
-    static const operation_type value = operation_intersection;
-};
-
 
 }} // namespace detail::overlay
 #endif //DOXYGEN_NO_DETAIL
diff --git a/boost/geometry/algorithms/detail/overlay/self_turn_points.hpp b/boost/geometry/algorithms/detail/overlay/self_turn_points.hpp
index a00606b087..3cc98d3b7e 100644
--- a/boost/geometry/algorithms/detail/overlay/self_turn_points.hpp
+++ b/boost/geometry/algorithms/detail/overlay/self_turn_points.hpp
@@ -77,7 +77,7 @@ struct self_section_visitor
     RobustPolicy const& m_rescale_policy;
     Turns& m_turns;
     InterruptPolicy& m_interrupt_policy;
-    std::size_t m_source_index;
+    int m_source_index;
     bool m_skip_adjacent;
 
     inline self_section_visitor(Geometry const& g,
@@ -85,7 +85,7 @@ struct self_section_visitor
                                 RobustPolicy const& rp,
                                 Turns& turns,
                                 InterruptPolicy& ip,
-                                std::size_t source_index,
+                                int source_index,
                                 bool skip_adjacent)
         : m_geometry(g)
         , m_intersection_strategy(is)
@@ -135,7 +135,7 @@ struct get_turns
             RobustPolicy const& robust_policy,
             Turns& turns,
             InterruptPolicy& interrupt_policy,
-            std::size_t source_index, bool skip_adjacent)
+            int source_index, bool skip_adjacent)
     {
         typedef model::box
             <
@@ -229,7 +229,7 @@ struct self_get_turn_points
             RobustPolicy const& ,
             Turns& ,
             InterruptPolicy& ,
-            std::size_t /*source_index*/,
+            int /*source_index*/,
             bool /*skip_adjacent*/)
     {
         return true;
@@ -292,7 +292,7 @@ inline void self_turns(Geometry const& geometry,
                        RobustPolicy const& robust_policy,
                        Turns& turns,
                        InterruptPolicy& interrupt_policy,
-                       std::size_t source_index = 0,
+                       int source_index = 0,
                        bool skip_adjacent = false)
 {
     concepts::check<Geometry const>();
@@ -339,7 +339,7 @@ inline void self_turns(Geometry const& geometry,
                        RobustPolicy const& robust_policy,
                        Turns& turns,
                        InterruptPolicy& interrupt_policy,
-                       std::size_t source_index = 0,
+                       int source_index = 0,
                        bool skip_adjacent = false)
 {
     concepts::check<Geometry const>();
diff --git a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
index fea5698ae1..6b929373b4 100644
--- a/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
+++ b/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp
@@ -1,6 +1,7 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
 // Copyright (c) 2015 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -239,7 +240,7 @@ public :
     {}
 
     template <typename Operation, typename Geometry1, typename Geometry2>
-    Point add(Operation const& op, signed_size_type turn_index, signed_size_type op_index,
+    Point add(Operation const& op, signed_size_type turn_index, int op_index,
             Geometry1 const& geometry1,
             Geometry2 const& geometry2,
             bool is_origin)
@@ -260,7 +261,7 @@ public :
     }
 
     template <typename Operation, typename Geometry1, typename Geometry2>
-    void add(Operation const& op, signed_size_type turn_index, signed_size_type op_index,
+    void add(Operation const& op, signed_size_type turn_index, int op_index,
             segment_identifier const& departure_seg_id,
             Geometry1 const& geometry1,
             Geometry2 const& geometry2,
@@ -277,7 +278,7 @@ public :
 
             if (is_origin)
             {
-                int const segment_distance = calculate_segment_distance(op, departure_seg_id, geometry1, geometry2);
+                signed_size_type const segment_distance = calculate_segment_distance(op, departure_seg_id, geometry1, geometry2);
                 if (m_origin_count == 0 ||
                         segment_distance < m_origin_segment_distance)
                 {
@@ -290,7 +291,7 @@ public :
     }
 
     template <typename Operation, typename Geometry1, typename Geometry2>
-    static int calculate_segment_distance(Operation const& op,
+    static signed_size_type calculate_segment_distance(Operation const& op,
             segment_identifier const& departure_seg_id,
             Geometry1 const& geometry1,
             Geometry2 const& geometry2)
@@ -303,7 +304,7 @@ public :
         // Suppose ring_count=10 (10 points, 9 segments), dep.seg_id=7, op.seg_id=2, then distance=10-9+2
         // Generic function (is this used somewhere else too?)
         ring_identifier const rid(op.seg_id.source_index, op.seg_id.multi_index, op.seg_id.ring_index);
-        int const segment_count
+        signed_size_type const segment_count
                     (op.seg_id.source_index == 0
                     ? geometry::num_points(detail::overlay::get_ring<typename geometry::tag<Geometry1>::type>::apply(rid, geometry1))
                     : geometry::num_points(detail::overlay::get_ring<typename geometry::tag<Geometry2>::type>::apply(rid, geometry2)));
@@ -434,7 +435,7 @@ public :
     container_type m_ranked_points;
     Point m_origin;
     std::size_t m_origin_count;
-    int m_origin_segment_distance;
+    signed_size_type m_origin_segment_distance;
     SideStrategy m_strategy;
 
 private :
diff --git a/boost/geometry/algorithms/detail/overlay/traversal.hpp b/boost/geometry/algorithms/detail/overlay/traversal.hpp
index 6a9b1def99..5c547c3278 100644
--- a/boost/geometry/algorithms/detail/overlay/traversal.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traversal.hpp
@@ -18,13 +18,12 @@
 
 #include <boost/range.hpp>
 
-#include <boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp>
 #include <boost/geometry/algorithms/detail/overlay/is_self_turn.hpp>
 #include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
-#include <boost/geometry/algorithms/detail/overlay/traversal_intersection_patterns.hpp>
 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
 #include <boost/geometry/core/access.hpp>
 #include <boost/geometry/core/assert.hpp>
+#include <boost/geometry/util/condition.hpp>
 
 #if defined(BOOST_GEOMETRY_DEBUG_INTERSECTION) \
     || defined(BOOST_GEOMETRY_OVERLAY_REPORT_WKT) \
@@ -232,50 +231,54 @@ struct traversal
     }
 
     template <signed_size_type segment_identifier::*Member>
-    inline bool select_source_generic(bool switch_source,
+    inline bool select_source_generic(turn_type const& turn,
             segment_identifier const& current,
             segment_identifier const& previous) const
     {
+        turn_operation_type const& op0 = turn.operations[0];
+        turn_operation_type const& op1 = turn.operations[1];
+
+        bool const switch_source = op0.enriched.region_id != -1
+                && op0.enriched.region_id == op1.enriched.region_id;
+
+#if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
+        if (switch_source)
+        {
+            std::cout << "Switch source at " << &turn << std::endl;
+        }
+        else
+        {
+            std::cout << "DON'T SWITCH SOURCES at " << &turn << std::endl;
+        }
+#endif
         return switch_source
                 ? current.*Member != previous.*Member
                 : current.*Member == previous.*Member;
     }
 
-    inline bool select_source(signed_size_type turn_index,
+    inline bool select_source(turn_type const& turn,
                               segment_identifier const& candidate_seg_id,
                               segment_identifier const& previous_seg_id) const
     {
         // For uu/ii, only switch sources if indicated
-        turn_type const& turn = m_turns[turn_index];
 
-#if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
-        if (turn.switch_source)
-        {
-            std::cout << "Switch source at " << turn_index << std::endl;
-        }
-        else
-        {
-            std::cout << "DON'T SWITCH SOURCES at " << turn_index << std::endl;
-        }
-#endif
-        if (OverlayType == overlay_buffer
-                || OverlayType == overlay_dissolve_union)
+        if (OverlayType == overlay_buffer)
         {
             // Buffer does not use source_index (always 0).
             return select_source_generic<&segment_identifier::multi_index>(
-                        turn.switch_source, candidate_seg_id, previous_seg_id);
+                        turn, candidate_seg_id, previous_seg_id);
         }
 
         if (is_self_turn<OverlayType>(turn))
         {
             // Also, if it is a self-turn, stay on same ring (multi/ring)
             return select_source_generic<&segment_identifier::multi_index>(
-                        turn.switch_source, candidate_seg_id, previous_seg_id);
+                        turn, candidate_seg_id, previous_seg_id);
         }
 
         // Use source_index
         return select_source_generic<&segment_identifier::source_index>(
-                    turn.switch_source, candidate_seg_id, previous_seg_id);
+                    turn, candidate_seg_id, previous_seg_id);
     }
 
     inline bool traverse_possible(signed_size_type turn_index) const
@@ -294,6 +297,51 @@ struct traversal
             || turn.has(operation_continue);
     }
 
+    inline std::size_t get_shortcut_level(turn_operation_type const& op,
+                             signed_size_type start_turn_index,
+                             signed_size_type origin_turn_index,
+                             std::size_t level = 1) const
+    {
+        signed_size_type next_turn_index = op.enriched.get_next_turn_index();
+        if (next_turn_index == -1)
+        {
+            return 0;
+        }
+        if (next_turn_index == start_turn_index)
+        {
+            // This operation finishes the ring
+            return 0;
+        }
+        if (next_turn_index == origin_turn_index)
+        {
+            // This operation travels to itself
+            return level;
+        }
+        if (level > 10)
+        {
+            // Avoid infinite recursion
+            return 0;
+        }
+
+        turn_type const& next_turn = m_turns[next_turn_index];
+        for (int i = 0; i < 2; i++)
+        {
+            turn_operation_type const& next_op = next_turn.operations[i];
+            if (next_op.operation == target_operation
+                && ! next_op.visited.finished()
+                && ! next_op.visited.visited())
+            {
+                // Recursively continue verifying
+                if (get_shortcut_level(next_op, start_turn_index,
+                                       origin_turn_index, level + 1))
+                {
+                    return level + 1;
+                }
+            }
+        }
+        return 0;
+    }
+
     inline
     bool select_cc_operation(turn_type const& turn,
                 signed_size_type start_turn_index,
@@ -342,7 +390,6 @@ struct traversal
 
     inline
     bool select_noncc_operation(turn_type const& turn,
-                signed_size_type turn_index,
                 segment_identifier const& previous_seg_id,
                 int& selected_op_index) const
     {
@@ -355,7 +402,7 @@ struct traversal
             if (op.operation == target_operation
                 && ! op.visited.finished()
                 && ! op.visited.visited()
-                && (! result || select_source(turn_index, op.seg_id, previous_seg_id)))
+                && (! result || select_source(turn, op.seg_id, previous_seg_id)))
             {
                 selected_op_index = i;
                 debug_traverse(turn, op, "Candidate");
@@ -367,6 +414,87 @@ struct traversal
     }
 
     inline
+    bool select_preferred_operation(turn_type const& turn,
+                signed_size_type turn_index,
+                signed_size_type start_turn_index,
+                int& selected_op_index) const
+    {
+        bool option[2] = {0};
+        bool finishing[2] = {0};
+        bool preferred[2] = {0};
+        std::size_t shortcut_level[2] = {0};
+        for (int i = 0; i < 2; i++)
+        {
+            turn_operation_type const& op = turn.operations[i];
+
+            if (op.operation == target_operation
+                && ! op.visited.finished()
+                && ! op.visited.visited())
+            {
+                option[i] = true;
+                if (op.enriched.get_next_turn_index() == start_turn_index)
+                {
+                    finishing[i] = true;
+                }
+                else
+                {
+                    shortcut_level[i] = get_shortcut_level(op, start_turn_index,
+                                                           turn_index);
+                }
+
+                if (op.enriched.prefer_start)
+                {
+                    preferred[i] = true;
+                }
+            }
+        }
+
+        if (option[0] != option[1])
+        {
+            // Only one operation is acceptable, take that one
+            selected_op_index = option[0] ? 0 : 1;
+            return true;
+        }
+
+        if (option[0] && option[1])
+        {
+            // Both operations are acceptable
+            if (finishing[0] != finishing[1])
+            {
+                // Prefer operation finishing the ring
+                selected_op_index = finishing[0] ? 0 : 1;
+                return true;
+            }
+
+            if (shortcut_level[0] != shortcut_level[1])
+            {
+                // If a turn can travel to itself again (without closing the
+                // ring), take the shortest one
+                selected_op_index = shortcut_level[0] < shortcut_level[1] ? 0 : 1;
+                return true;
+            }
+
+            if (preferred[0] != preferred[1])
+            {
+                // Only one operation is preferred (== was not intersection)
+                selected_op_index = preferred[0] ? 0 : 1;
+                return true;
+            }
+        }
+
+        for (int i = 0; i < 2; i++)
+        {
+            if (option[i])
+            {
+                selected_op_index = 0;
+                return true;
+            }
+        }
+
+        return false;
+    }
+
+    inline
     bool select_operation(const turn_type& turn,
                 signed_size_type turn_index,
                 signed_size_type start_turn_index,
@@ -380,10 +508,15 @@ struct traversal
             result = select_cc_operation(turn, start_turn_index,
                                          selected_op_index);
         }
+        else if (OverlayType == overlay_dissolve)
+        {
+            result = select_preferred_operation(turn, turn_index,
+                start_turn_index, selected_op_index);
+        }
         else
         {
-            result = select_noncc_operation(turn, turn_index,
-                                            previous_seg_id, selected_op_index);
+            result = select_noncc_operation(turn, previous_seg_id,
+                selected_op_index);
         }
         if (result)
         {
@@ -417,6 +550,13 @@ struct traversal
         return true;
     }
 
+
+    template <typename RankedPoint>
+    inline turn_operation_type const& operation_from_rank(RankedPoint const& rp) const
+    {
+        return m_turns[rp.turn_index].operations[rp.operation_index];
+    }
+
     inline int select_turn_in_cluster_union(std::size_t selected_rank,
             typename sbs_type::rp const& ranked_point,
             signed_size_type start_turn_index, int start_op_index) const
@@ -431,8 +571,7 @@ struct traversal
             return 0;
         }
 
-        turn_type const& turn = m_turns[ranked_point.turn_index];
-        turn_operation_type const& op = turn.operations[ranked_point.operation_index];
+        turn_operation_type const& op = operation_from_rank(ranked_point);
 
         // Check finalized: TODO: this should be finetuned, it is not necessary
         if (op.visited.finalized())
@@ -440,7 +579,7 @@ struct traversal
             return 0;
         }
 
-        if (OverlayType != overlay_dissolve_union
+        if (OverlayType != overlay_dissolve
             && (op.enriched.count_left != 0 || op.enriched.count_right == 0))
         {
             // Check counts: in some cases interior rings might be generated with
@@ -455,27 +594,45 @@ struct traversal
             ;
     }
 
-    inline bool select_from_cluster_union(signed_size_type& turn_index,
-        int& op_index, sbs_type& sbs,
-        signed_size_type start_turn_index, int start_op_index) const
+    inline signed_size_type select_rank(sbs_type const& sbs,
+                                        bool skip_isolated) const
     {
-        std::vector<sort_by_side::rank_with_rings> aggregation;
-        sort_by_side::aggregate_operations(sbs, aggregation, m_turns, operation_union);
-
-        sort_by_side::rank_with_rings const& incoming = aggregation.front();
+        // Take the first outgoing rank corresponding to incoming region,
+        // or take another region if it is not isolated
+        turn_operation_type const& incoming_op
+                = operation_from_rank(sbs.m_ranked_points.front());
 
-        // Take the first one outgoing for the incoming region
-        std::size_t selected_rank = 0;
-        for (std::size_t i = 1; i < aggregation.size(); i++)
+        for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
         {
-            sort_by_side::rank_with_rings const& rwr = aggregation[i];
-            if (rwr.all_to()
-                    && rwr.region_id() == incoming.region_id())
+            typename sbs_type::rp const& rp = sbs.m_ranked_points[i];
+            if (rp.rank == 0 || rp.direction == sort_by_side::dir_from)
             {
-                selected_rank = rwr.rank;
-                break;
+                continue;
+            }
+            turn_operation_type const& op = operation_from_rank(rp);
+
+            if (op.operation != target_operation
+                && op.operation != operation_continue)
+            {
+                continue;
+            }
+
+            if (op.enriched.region_id == incoming_op.enriched.region_id
+                || (skip_isolated && ! op.enriched.isolated))
+            {
+                // Region corresponds to incoming region, or (for intersection)
+                // there is a non-isolated other region which should be taken
+                return rp.rank;
             }
         }
+        return -1;
+    }
+
+    inline bool select_from_cluster_union(signed_size_type& turn_index,
+        int& op_index, sbs_type const& sbs,
+        signed_size_type start_turn_index, int start_op_index) const
+    {
+        std::size_t const selected_rank = select_rank(sbs, false);
 
         int best_code = 0;
         bool result = false;
@@ -508,87 +665,7 @@ struct traversal
     inline bool analyze_cluster_intersection(signed_size_type& turn_index,
                 int& op_index, sbs_type const& sbs) const
     {
-        std::vector<sort_by_side::rank_with_rings> aggregation;
-        sort_by_side::aggregate_operations(sbs, aggregation, m_turns, operation_intersection);
-
-        std::size_t selected_rank = 0;
-
-
-        // Detect specific pattern(s)
-        bool const detected
-            = intersection_pattern_common_interior1(selected_rank, aggregation)
-            || intersection_pattern_common_interior2(selected_rank, aggregation)
-            || intersection_pattern_common_interior3(selected_rank, aggregation)
-            || intersection_pattern_common_interior4(selected_rank, aggregation)
-            || intersection_pattern_common_interior5(selected_rank, aggregation)
-            || intersection_pattern_common_interior6(selected_rank, aggregation)
-                ;
-
-        if (! detected)
-        {
-            signed_size_type incoming_region_id = 0;
-            std::set<signed_size_type> outgoing_region_ids;
-
-            for (std::size_t i = 0; i < aggregation.size(); i++)
-            {
-                sort_by_side::rank_with_rings const& rwr = aggregation[i];
-
-                if (rwr.all_to()
-                        && rwr.traversable(m_turns)
-                        && selected_rank == 0)
-                {
-                    // Take the first (= right) where segments leave,
-                    // having the polygon on the right side
-                    selected_rank = rwr.rank;
-                }
-
-                if (rwr.all_from()
-                        && selected_rank > 0
-                        && outgoing_region_ids.empty())
-                {
-                    // Incoming
-                    break;
-                }
-
-                if (incoming_region_id == 0)
-                {
-                    sort_by_side::ring_with_direction const& rwd = *rwr.rings.begin();
-                    turn_type const& turn = m_turns[rwd.turn_index];
-                    incoming_region_id = turn.operations[rwd.operation_index].enriched.region_id;
-                }
-                else
-                {
-                    if (rwr.rings.size() == 1)
-                    {
-                        sort_by_side::ring_with_direction const& rwd = *rwr.rings.begin();
-                        turn_type const& turn = m_turns[rwd.turn_index];
-                        if (rwd.direction == sort_by_side::dir_to
-                                && turn.both(operation_intersection))
-                        {
-
-                            turn_operation_type const& op = turn.operations[rwd.operation_index];
-                            if (op.enriched.region_id != incoming_region_id
-                                    && op.enriched.isolated)
-                            {
-                                outgoing_region_ids.insert(op.enriched.region_id);
-                            }
-                        }
-                        else if (! outgoing_region_ids.empty())
-                        {
-                            for (int i = 0; i < 2; i++)
-                            {
-                                signed_size_type const region_id = turn.operations[i].enriched.region_id;
-                                if (outgoing_region_ids.count(region_id) == 1)
-                                {
-                                    selected_rank = 0;
-                                    outgoing_region_ids.erase(region_id);
-                                }
-                            }
-                        }
-                    }
-                }
-            }
-        }
+        std::size_t const selected_rank = select_rank(sbs, true);
 
         if (selected_rank > 0)
         {
@@ -602,10 +679,9 @@ struct traversal
 
                 if (ranked_point.rank == selected_rank)
                 {
-                    turn_type const& ranked_turn = m_turns[ranked_point.turn_index];
-                    turn_operation_type const& ranked_op = ranked_turn.operations[ranked_point.operation_index];
+                    turn_operation_type const& op = operation_from_rank(ranked_point);
 
-                    if (ranked_op.visited.finalized())
+                    if (op.visited.finalized())
                     {
                         // This direction is already traveled before, the same
                         // cannot be traveled again
@@ -614,10 +690,10 @@ struct traversal
 
                     // Take turn with the smallest remaining distance
                     if (selected_index == sbs.m_ranked_points.size()
-                            || ranked_op.remaining_distance < min_remaining_distance)
+                            || op.remaining_distance < min_remaining_distance)
                     {
                         selected_index = i;
-                        min_remaining_distance = ranked_op.remaining_distance;
+                        min_remaining_distance = op.remaining_distance;
                     }
                 }
             }
@@ -725,9 +801,7 @@ struct traversal
                 turn_operation_type const& start_op,
                 int start_op_index) const
     {
-        if (OverlayType != overlay_buffer
-                && OverlayType != overlay_dissolve_union
-                && OverlayType != overlay_dissolve_intersection)
+        if (OverlayType != overlay_buffer && OverlayType != overlay_dissolve)
         {
             return;
         }
@@ -830,7 +904,7 @@ struct traversal
     {
         turn_type const& current_turn = m_turns[turn_index];
 
-        if (target_operation == operation_intersection)
+        if (BOOST_GEOMETRY_CONDITION(target_operation == operation_intersection))
         {
             bool const back_at_start_cluster
                     = current_turn.is_clustered()
diff --git a/boost/geometry/algorithms/detail/overlay/traversal_intersection_patterns.hpp b/boost/geometry/algorithms/detail/overlay/traversal_intersection_patterns.hpp
deleted file mode 100644
index a8abea230e..0000000000
--- a/boost/geometry/algorithms/detail/overlay/traversal_intersection_patterns.hpp
+++ /dev/null
@@ -1,451 +0,0 @@
-// Boost.Geometry (aka GGL, Generic Geometry Library)
-
-// Copyright (c) 2007-2017 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_TRAVERSAL_INTERSECTION_PATTERNS_HPP
-#define BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_TRAVERSAL_INTERSECTION_PATTERNS_HPP
-
-#include <cstddef>
-#include <vector>
-
-#include <boost/geometry/algorithms/detail/overlay/aggregate_operations.hpp>
-#include <boost/geometry/algorithms/detail/overlay/sort_by_side.hpp>
-
-namespace boost { namespace geometry
-{
-
-#ifndef DOXYGEN_NO_DETAIL
-namespace detail { namespace overlay
-{
-
-inline bool check_pairs(std::vector<sort_by_side::rank_with_rings> const& aggregation,
-                        signed_size_type incoming_region_id,
-                        std::size_t first, std::size_t last)
-{
-    // Check if pairs 1,2 (and possibly 3,4 and 5,6 etc) satisfy
-
-    for (std::size_t i = first; i <= last; i += 2)
-    {
-        sort_by_side::rank_with_rings const& curr = aggregation[i];
-        sort_by_side::rank_with_rings const& next = aggregation[i + 1];
-        signed_size_type const curr_id = curr.region_id();
-        signed_size_type const next_id = next.region_id();
-
-        bool const possible =
-                curr.rings.size() == 2
-                && curr.is_isolated()
-                && curr.has_unique_region_id()
-                && next.rings.size() == 2
-                && next.is_isolated()
-                && next.has_unique_region_id()
-                && curr_id == next_id
-                && curr_id != incoming_region_id;
-
-        if (! possible)
-        {
-            return false;
-        }
-    }
-
-    return true;
-}
-
-inline bool intersection_pattern_common_interior1(std::size_t& selected_rank,
-           std::vector<sort_by_side::rank_with_rings> const& aggregation)
-{
-    // Pattern: coming from exterior ring, encountering an isolated
-    // parallel interior ring, which should be skipped, and the first
-    // left (normally intersection takes first right) should be taken.
-    // Solves cases #case_133_multi
-    // and #case_recursive_boxes_49
-
-    std::size_t const n = aggregation.size();
-    if (n < 4)
-    {
-        return false;
-    }
-
-    sort_by_side::rank_with_rings const& incoming = aggregation.front();
-    sort_by_side::rank_with_rings const& outgoing = aggregation.back();
-
-    bool const incoming_ok =
-        incoming.all_from()
-        && incoming.rings.size() == 1
-        && incoming.has_only(operation_intersection);
-
-    if (! incoming_ok)
-    {
-        return false;
-    }
-
-    bool const outgoing_ok =
-        outgoing.all_to()
-        && outgoing.rings.size() == 1
-        && outgoing.has_only(operation_intersection)
-        && outgoing.region_id() == incoming.region_id();
-
-    if (! outgoing_ok)
-    {
-        return false;
-    }
-
-    if (check_pairs(aggregation, incoming.region_id(), 1, n - 2))
-    {
-        selected_rank = n - 1;
-        return true;
-    }
-    return false;
-}
-
-inline bool intersection_pattern_common_interior2(std::size_t& selected_rank,
-           std::vector<sort_by_side::rank_with_rings> const& aggregation)
-{
-    // Pattern: coming from two exterior rings, encountering two isolated
-    // equal interior rings
-
-    // See (for example, for ii) #case_recursive_boxes_53:
-
-    // INCOMING:
-    // Rank 0  {11[0] (s:0, m:0) i F rgn: 1 ISO}             {13[1] (s:1, m:0) i F rgn: 1 ISO}
-
-    // PAIR:
-    // Rank 1  {13[0] (s:0, r:1, m:0) i T rgn: 3 ISO ->16}   {11[1] (s:1, r:5, m:0) i T rgn: 3 ISO ->16}
-    // Rank 2  {13[0] (s:0, r:1, m:0) i F rgn: 3 ISO}        {11[1] (s:1, r:5, m:0) i F rgn: 3 ISO}
-
-    // LEAVING (in the same direction, take last one)
-    // Rank 3  {11[0] (s:0, m:0) i T rgn: 1 ISO ->10}        {13[1] (s:1, m:0) i T rgn: 1 ISO ->10}
-
-
-    std::size_t const n = aggregation.size();
-    if (n < 4)
-    {
-        return false;
-    }
-
-    sort_by_side::rank_with_rings const& incoming = aggregation.front();
-    sort_by_side::rank_with_rings const& outgoing = aggregation.back();
-
-    bool const incoming_ok =
-        incoming.all_from()
-        && incoming.rings.size() == 2
-        && incoming.has_unique_region_id();
-
-    if (! incoming_ok)
-    {
-        return false;
-    }
-
-    bool const outgoing_ok =
-        outgoing.all_to()
-        && outgoing.rings.size() == 2
-        && outgoing.has_unique_region_id()
-        && outgoing.region_id() == incoming.region_id();
-
-    if (! outgoing_ok)
-    {
-        return false;
-    }
-
-    bool const operation_ok =
-            (incoming.has_only(operation_continue) && outgoing.has_only(operation_continue))
-             || (incoming.has_only(operation_intersection) && outgoing.has_only(operation_intersection));
-
-    if (! operation_ok)
-    {
-        return false;
-    }
-
-    // Check if pairs 1,2 (and possibly 3,4 and 5,6 etc) satisfy
-    if (check_pairs(aggregation, incoming.region_id(), 1, n - 2))
-    {
-        selected_rank = n - 1;
-        return true;
-    }
-    return false;
-}
-
-inline bool intersection_pattern_common_interior3(std::size_t& selected_rank,
-           std::vector<sort_by_side::rank_with_rings> const& aggregation)
-{
-    // Pattern: approaches colocated turn (exterior+interior) from two
-    // different directions, and both leaves in the same direction
-
-    // See #case_136_multi:
-    // INCOMING:
-    //Rank 0  {10[0] (s:0, m:0) c F rgn: 1 ISO}
-
-    // PAIR:
-    //Rank 1  {14[0] (s:0, r:0, m:0) i T rgn: 2 ISO ->16} {11[1] (s:1, r:1, m:0) i T rgn: 2 ISO ->16}
-    //Rank 2  {14[0] (s:0, r:0, m:0) i F rgn: 2 ISO}      {11[1] (s:1, r:1, m:0) i F rgn: 2 ISO}
-
-    // LEAVING (select this one):
-    //Rank 3  {10[0] (s:0, m:0) c T rgn: 1 ISO ->12}      {10[1] (s:1, m:0) c T rgn: 1 ISO ->12}
-
-    // ADDITIONALLY: (other polygon coming in)
-    //Rank 4  {10[1] (s:1, m:0) c F rgn: 1 ISO}
-
-    std::size_t const n = aggregation.size();
-    if (n < 4)
-    {
-        return false;
-    }
-
-    sort_by_side::rank_with_rings const& incoming = aggregation.front();
-    sort_by_side::rank_with_rings const& outgoing = aggregation[n - 2];
-    sort_by_side::rank_with_rings const& last = aggregation.back();
-
-    bool const incoming_ok =
-        incoming.all_from()
-        && incoming.rings.size() == 1
-        && incoming.has_only(operation_continue);
-
-    if (! incoming_ok)
-    {
-        return false;
-    }
-
-    bool const outgoing_ok =
-        outgoing.all_to()
-        && outgoing.rings.size() == 2
-        && outgoing.has_only(operation_continue)
-        && outgoing.has_unique_region_id()
-        && outgoing.region_id() == incoming.region_id()
-        && last.all_from()
-        && last.rings.size() == 1
-        && last.region_id() == incoming.region_id()
-        && last.all_from();
-
-    if (! outgoing_ok)
-    {
-        return false;
-    }
-
-    // Check if pairs 1,2 (and possibly 3,4 and 5,6 etc) satisfy
-    if (check_pairs(aggregation, incoming.region_id(), 1, n - 3))
-    {
-        selected_rank = n - 2;
-        return true;
-    }
-    return false;
-}
-
-
-inline bool intersection_pattern_common_interior4(std::size_t& selected_rank,
-           std::vector<sort_by_side::rank_with_rings> const& aggregation)
-{
-    // Pattern: approaches colocated turn (exterior+interior) from same
-    // direction, but leaves in two different directions
-
-    // See #case_137_multi:
-
-    // INCOMING:
-    //Rank 0  {11[0] (s:0, m:0) i F rgn: 1 ISO}            {10[1] (s:1, m:0) i F rgn: 1 ISO}
-
-    // PAIR:
-    //Rank 1  {13[0] (s:0, r:0, m:0) i T rgn: 2 ISO ->15}  {11[1] (s:1, r:1, m:0) i T rgn: 2 ISO ->15}
-    //Rank 2  {13[0] (s:0, r:0, m:0) i F rgn: 2 ISO}       {11[1] (s:1, r:1, m:0) i F rgn: 2 ISO}
-
-    // LEAVING (in two different directions, take penultimate one)
-    //Rank 3  {10[1] (s:1, m:0) i T rgn: 1 ISO ->0}
-    //Rank 4  {11[0] (s:0, m:0) i T rgn: 1 ISO ->12}
-
-    std::size_t const n = aggregation.size();
-    if (n < 4)
-    {
-        return false;
-    }
-
-    sort_by_side::rank_with_rings const& incoming = aggregation.front();
-    sort_by_side::rank_with_rings const& extra = aggregation[n - 2];
-    sort_by_side::rank_with_rings const& outgoing = aggregation.back();
-
-    bool const incoming_ok =
-        incoming.all_from()
-        && incoming.rings.size() == 2
-        && incoming.has_unique_region_id()
-        && incoming.has_only(operation_intersection);
-
-    if (! incoming_ok)
-    {
-        return false;
-    }
-
-    bool const outgoing_ok =
-        outgoing.all_to()
-        && outgoing.rings.size() == 1
-        && outgoing.has_only(operation_intersection)
-        && outgoing.region_id() == incoming.region_id()
-        && extra.all_to()
-        && extra.rings.size() == 1
-        && extra.has_only(operation_intersection)
-        && extra.region_id() == incoming.region_id();
-
-    if (! outgoing_ok)
-    {
-        return false;
-    }
-
-    // Check if pairs 1,2 (and possibly 3,4 and 5,6 etc) satisfy
-    if (check_pairs(aggregation, incoming.region_id(), 1, n - 3))
-    {
-        selected_rank = n - 2;
-        return true;
-    }
-    return false;
-}
-
-inline bool intersection_pattern_common_interior5(std::size_t& selected_rank,
-           std::vector<sort_by_side::rank_with_rings> const& aggregation)
-{
-    // Pattern: isolated regions
-
-    // See #case_recursive_boxes_65
-
-    // INCOMING:
-    //    Rank 0  {19[0] (s:0, r:2, m:0) i F rgn: 4 ISO}
-
-    //    Rank 1  {19[1] (s:1, m:0) i T rgn: 1 ISO FIN ->18 (2.5)}
-    //    Rank 2  {21[1] (s:1, m:2) i F rgn: 1 ISO}
-    //    Rank 3  {21[1] (s:1, m:2) i T rgn: 1 ISO ->17 (2)}
-    //    Rank 4  {19[1] (s:1, m:0) i F rgn: 1 ISO FIN}
-
-    // LEAVING (take this one):
-    //    Rank 5  {19[0] (s:0, r:2, m:0) i T rgn: 4 ISO ->22 (1)}
-
-    std::size_t const n = aggregation.size();
-    if (n < 3)
-    {
-        return false;
-    }
-
-    sort_by_side::rank_with_rings const& incoming = aggregation.front();
-    sort_by_side::rank_with_rings const& outgoing = aggregation.back();
-
-    bool const incoming_ok =
-        incoming.all_from()
-        && incoming.has_unique_region_id()
-        && incoming.is_isolated();
-
-    if (! incoming_ok)
-    {
-        return false;
-    }
-
-    signed_size_type const incoming_region_id = incoming.region_id();
-
-    bool const outgoing_ok =
-        outgoing.all_to()
-        && outgoing.has_unique_region_id()
-        && outgoing.is_isolated()
-        && outgoing.region_id() == incoming_region_id;
-
-    if (! outgoing_ok)
-    {
-        return false;
-    }
-
-    selected_rank = n - 1;
-    bool other_region = true;
-
-    // Assumed is that other regions go (T) and come back (F)
-    for (std::size_t i = 1; i < n - 1; i++)
-    {
-        sort_by_side::rank_with_rings const& rwr = aggregation[i];
-        if (! rwr.has_unique_region_id() || ! rwr.is_isolated())
-        {
-            return false;
-        }
-        signed_size_type const region_id = rwr.region_id();
-        if (other_region && region_id != incoming_region_id)
-        {
-            // OK
-        }
-        else if (other_region && region_id == incoming_region_id)
-        {
-            // OK, next phase (same region as incoming region)
-            selected_rank = i;
-            other_region = false;
-        }
-        else if (! other_region && region_id != incoming_region_id)
-        {
-            // After that the region is the same is incoming, it should
-            // stay like that
-            return false;
-        }
-    }
-
-    return true;
-}
-
-inline bool intersection_pattern_common_interior6(std::size_t& selected_rank,
-           std::vector<sort_by_side::rank_with_rings> const& aggregation)
-{
-    // Pattern: isolated regions in between
-
-    // See #case_recursive_boxes_75
-
-    // Incoming: one region
-    // In between: several rings having isolated region, all the same
-    // Outging == incoming
-
-    std::size_t const n = aggregation.size();
-    if (n < 3)
-    {
-        return false;
-    }
-
-    sort_by_side::rank_with_rings const& incoming = aggregation.front();
-    sort_by_side::rank_with_rings const& outgoing = aggregation.back();
-    sort_by_side::rank_with_rings const& first_isolated = aggregation[2];
-
-    bool const incoming_ok =
-        incoming.all_from()
-        && incoming.has_unique_region_id()
-        && ! incoming.is_isolated();
-
-    if (! incoming_ok)
-    {
-        return false;
-    }
-
-    signed_size_type const incoming_region_id = incoming.region_id();
-
-    bool const outgoing_ok =
-        outgoing.all_to()
-        && outgoing.has_unique_region_id()
-        && ! outgoing.is_isolated()
-        && outgoing.region_id() == incoming_region_id;
-
-    if (! outgoing_ok)
-    {
-        return false;
-    }
-
-    const signed_size_type isolated_region_id = first_isolated.region_id();
-
-    for (std::size_t i = 1; i < n - 1; i++)
-    {
-        sort_by_side::rank_with_rings const& rwr = aggregation[i];
-        if (! rwr.has_unique_region_id()
-                || ! rwr.is_isolated()
-                || rwr.region_id() != isolated_region_id)
-        {
-            return false;
-        }
-    }
-
-    selected_rank = n - 1;
-
-    return true;
-}
-
-}} // namespace detail::overlay
-#endif // DOXYGEN_NO_DETAIL
-
-}} // namespace boost::geometry
-
-#endif // BOOST_GEOMETRY_ALGORITHMS_DETAIL_OVERLAY_TRAVERSAL_INTERSECTION_PATTERNS_HPP
diff --git a/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp b/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp
index 4df3f6e7ac..7f80c8313a 100644
--- a/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traversal_ring_creator.hpp
@@ -162,7 +162,7 @@ struct traversal_ring_creator
         // Update registration and append point
         turn_type& current_turn = m_turns[turn_index];
         turn_operation_type& op = current_turn.operations[op_index];
-        detail::overlay::append_no_dups_or_spikes(current_ring, current_turn.point,
+        detail::overlay::append_no_collinear(current_ring, current_turn.point,
             m_intersection_strategy.get_side_strategy(),
             m_robust_policy);
 
@@ -180,7 +180,7 @@ struct traversal_ring_creator
         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,
+        detail::overlay::append_no_collinear(ring, start_turn.point,
             m_intersection_strategy.get_side_strategy(),
             m_robust_policy);
 
@@ -344,6 +344,60 @@ struct traversal_ring_creator
         }
     }
 
+    template <typename Rings>
+    void iterate_with_preference(std::size_t phase,
+                 Rings& rings, std::size_t& finalized_ring_size,
+                 typename Backtrack::state_type& state)
+    {
+        for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index)
+        {
+            turn_type const& turn = m_turns[turn_index];
+
+            if (turn.discarded || turn.blocked())
+            {
+                // Skip discarded and blocked turns
+                continue;
+            }
+
+            turn_operation_type const& op0 = turn.operations[0];
+            turn_operation_type const& op1 = turn.operations[1];
+
+            if (phase == 0)
+            {
+                if (! op0.enriched.prefer_start && ! op1.enriched.prefer_start)
+                {
+                    // Not preferred, take next one
+                    continue;
+                }
+            }
+
+            if (turn.both(operation_continue))
+            {
+                // Traverse only one turn, the one with the SMALLEST remaining distance
+                // to avoid skipping a turn in between, which can happen in rare cases
+                // (e.g. #130)
+                int const op_index
+                        = op0.remaining_distance <= op1.remaining_distance ? 0 : 1;
+
+                traverse_with_operation(turn, turn_index, op_index,
+                        rings, finalized_ring_size, state);
+            }
+            else
+            {
+                bool const forward = op0.enriched.prefer_start;
+
+                int op_index = forward ? 0 : 1;
+                int const increment = forward ? 1 : -1;
+
+                for (int i = 0; i < 2; i++, op_index += increment)
+                {
+                    traverse_with_operation(turn, turn_index, op_index,
+                            rings, finalized_ring_size, state);
+                }
+            }
+        }
+    }
+
 private:
     traversal_type m_trav;
 
diff --git a/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp b/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp
index 72f9c4f12a..8bdb03df5d 100644
--- a/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp
+++ b/boost/geometry/algorithms/detail/overlay/traversal_switch_detector.hpp
@@ -20,6 +20,7 @@
 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
 #include <boost/geometry/core/access.hpp>
 #include <boost/geometry/core/assert.hpp>
+#include <boost/geometry/util/condition.hpp>
 
 namespace boost { namespace geometry
 {
@@ -48,6 +49,13 @@ template
 >
 struct traversal_switch_detector
 {
+    static const operation_type target_operation
+            = operation_from_overlay<OverlayType>::value;
+    static const operation_type opposite_operation
+            = target_operation == operation_union
+            ? operation_intersection
+            : operation_union;
+
     enum isolation_type
     {
         isolation_unknown = -1,
@@ -119,33 +127,6 @@ struct traversal_switch_detector
     {
     }
 
-    bool inspect_difference(set_type& turn_id_difference,
-            set_type const& turn_ids,
-            set_type const& other_turn_ids) const
-    {
-        // TODO: consider if std::set_difference can be used in the final version
-        int const turn_count = turn_ids.size();
-        int const other_turn_count = other_turn_ids.size();
-
-        // First quick check on size (TODO: implement multiple-multiple connections)
-        if (turn_count - other_turn_count > 1)
-        {
-            return false;
-        }
-
-        // Check if all turns are also present in the connection.
-        // The difference is returned
-        for (set_iterator it = turn_ids.begin(); it != turn_ids.end(); ++it)
-        {
-            signed_size_type const& id = *it;
-            if (other_turn_ids.count(id) == 0)
-            {
-                turn_id_difference.insert(id);
-            }
-        }
-        return true;
-    }
-
     bool one_connection_to_another_region(region_properties const& region) const
     {
         // For example:
@@ -224,12 +205,83 @@ struct traversal_switch_detector
         return true;
     }
 
+    bool ii_turn_connects_two_regions(region_properties const& region,
+            region_properties const& connected_region,
+            signed_size_type turn_index) const
+    {
+        turn_type const& turn = m_turns[turn_index];
+        if (! turn.both(operation_intersection))
+        {
+            return false;
+        }
+
+        signed_size_type const id0 = turn.operations[0].enriched.region_id;
+        signed_size_type const id1 = turn.operations[1].enriched.region_id;
+
+        return (id0 == region.region_id && id1 == connected_region.region_id)
+            || (id1 == region.region_id && id0 == connected_region.region_id);
+    }
+
+
+    bool isolated_multiple_connection(region_properties const& region,
+            region_properties const& connected_region) const
+    {
+        if (connected_region.isolated != isolation_multiple)
+        {
+            return false;
+        }
+
+        // First step: compare turns of regions with turns of connected region
+        set_type turn_ids = region.unique_turn_ids;
+        for (set_iterator sit = connected_region.unique_turn_ids.begin();
+             sit != connected_region.unique_turn_ids.end(); ++sit)
+        {
+            turn_ids.erase(*sit);
+        }
+
+        // There should be one connection (turn or cluster) left
+        if (turn_ids.size() != 1)
+        {
+            return false;
+        }
+
+        for (set_iterator sit = connected_region.unique_turn_ids.begin();
+             sit != connected_region.unique_turn_ids.end(); ++sit)
+        {
+            signed_size_type const id_or_index = *sit;
+            if (id_or_index >= 0)
+            {
+                if (! ii_turn_connects_two_regions(region, connected_region, id_or_index))
+                {
+                    return false;
+                }
+            }
+            else
+            {
+                signed_size_type const cluster_id = -id_or_index;
+                typename Clusters::const_iterator it = m_clusters.find(cluster_id);
+                if (it != m_clusters.end())
+                {
+                    cluster_info const& cinfo = it->second;
+                    for (set_iterator cit = cinfo.turn_indices.begin();
+                         cit != cinfo.turn_indices.end(); ++cit)
+                    {
+                        if (! ii_turn_connects_two_regions(region, connected_region, *cit))
+                        {
+                            return false;
+                        }
+                    }
+                }
+            }
+        }
+
+        return true;
+    }
+
     bool has_only_isolated_children(region_properties const& region) const
     {
         bool first_with_turn = true;
-        bool first_with_multiple = true;
         signed_size_type first_turn_id = 0;
-        signed_size_type first_multiple_region_id = 0;
 
         for (typename connection_map::const_iterator it = region.connected_region_counts.begin();
              it != region.connected_region_counts.end(); ++it)
@@ -246,43 +298,17 @@ struct traversal_switch_detector
 
             region_properties const& connected_region = mit->second;
 
-            bool const multiple = connected_region.isolated == isolation_multiple;
-
             if (cprop.count != 1)
             {
-                if (! multiple)
+                // If there are more connections, check their isolation
+                if (! isolated_multiple_connection(region, connected_region))
                 {
                     return false;
                 }
-
-                // It connects multiple times to an isolated region.
-                // This is allowed as long as it happens only once
-                if (first_with_multiple)
-                {
-                    first_multiple_region_id = connected_region.region_id;
-                    first_with_multiple = false;
-                }
-                else if (first_multiple_region_id != connected_region.region_id)
-                {
-                    return false;
-                }
-
-                // Turns in region should be either present in the connection,
-                // of form part of the connection with the other region
-                set_type diff;
-                if (! inspect_difference(diff, region.unique_turn_ids,
-                                  connected_region.unique_turn_ids))
-                {
-                    return false;
-                }
-                if (diff.size() > 1)
-                {
-                    // For now:
-                    return false;
-                }
             }
 
-            if (connected_region.isolated != isolation_yes && ! multiple)
+            if (connected_region.isolated != isolation_yes
+                && connected_region.isolated != isolation_multiple)
             {
                 signed_size_type const unique_turn_id = *cprop.unique_turn_ids.begin();
                 if (first_with_turn)
@@ -296,6 +322,7 @@ struct traversal_switch_detector
                 }
             }
         }
+
         // If there is only one connection (with a 'parent'), and all other
         // connections are itself isolated, it is isolated
         return true;
@@ -381,6 +408,12 @@ struct traversal_switch_detector
             {
                 turn_type& turn = m_turns[*sit];
 
+                if (! acceptable(turn))
+                {
+                    // No assignment necessary
+                    continue;
+                }
+
                 for (int i = 0; i < 2; i++)
                 {
                     turn_operation_type& op = turn.operations[i];
@@ -441,12 +474,19 @@ struct traversal_switch_detector
         }
     }
 
+    inline bool acceptable(turn_type const& turn) const
+    {
+        // Discarded turns don't connect rings to the same region
+        // Also xx are not relevant
+        // (otherwise discarded colocated uu turn could make a connection)
+        return ! turn.discarded
+            && ! turn.both(operation_blocked);
+    }
+
     inline bool connects_same_region(turn_type const& turn) const
     {
-        if (turn.discarded)
+        if (! acceptable(turn))
         {
-            // Discarded turns don't connect same region (otherwise discarded colocated uu turn
-            // could make a connection)
             return false;
         }
 
@@ -456,7 +496,7 @@ struct traversal_switch_detector
             return ! (turn.both(operation_union) || turn.both(operation_intersection));
         }
 
-        if (operation_from_overlay<OverlayType>::value == operation_union)
+        if (BOOST_GEOMETRY_CONDITION(target_operation == operation_union))
         {
             // It is a cluster, check zones
             // (assigned by sort_by_side/handle colocations) of both operations
@@ -540,7 +580,7 @@ struct traversal_switch_detector
     void iterate()
     {
 #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
-        std::cout << "SWITCH BEGIN ITERATION" << std::endl;
+        std::cout << "BEGIN ITERATION GETTING REGION_IDS" << std::endl;
 #endif
 
         // Collect turns per ring
@@ -552,7 +592,7 @@ struct traversal_switch_detector
             turn_type const& turn = m_turns[turn_index];
 
             if (turn.discarded
-                    && operation_from_overlay<OverlayType>::value == operation_intersection)
+                && BOOST_GEOMETRY_CONDITION(target_operation == operation_intersection))
             {
                 // Discarded turn (union currently still needs it to determine regions)
                 continue;
@@ -580,64 +620,8 @@ struct traversal_switch_detector
             assign_isolation();
         }
 
-        // Now that all regions are filled, assign switch_source property
-        // Iterate through all clusters
-        for (typename Clusters::iterator it = m_clusters.begin(); it != m_clusters.end(); ++it)
-        {
-            cluster_info& cinfo = it->second;
-            if (cinfo.open_count <= 1)
-            {
-                // Not a touching cluster
-                continue;
-            }
-
-            // A touching cluster, gather regions
-            set_type regions;
-            set_type const& ids = cinfo.turn_indices;
-
-#if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
-                std::cout << "SWITCH EXAMINE CLUSTER " << it->first << std::endl;
-#endif
-
-            for (set_iterator sit = ids.begin(); sit != ids.end(); ++sit)
-            {
-                signed_size_type turn_index = *sit;
-                turn_type const& turn = m_turns[turn_index];
-                for (int oi = 0; oi < 2; oi++)
-                {
-                    signed_size_type const region_id = get_region_id(turn.operations[oi]);
-                    regions.insert(region_id);
-                }
-            }
-            // Switch source if this cluster connects the same region
-            cinfo.switch_source = regions.size() <= 1;
-        }
-
-        // Iterate through all uu/ii turns (non-clustered)
-        for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index)
-        {
-            turn_type& turn = m_turns[turn_index];
-
-            if (turn.discarded
-                    || turn.blocked()
-                    || turn.is_clustered()
-                    || ! (turn.both(operation_union) || turn.both(operation_intersection)))
-            {
-                // Skip discarded, blocked, non-uu/ii and clustered turns
-                continue;
-            }
-
-
-            signed_size_type const region0 = get_region_id(turn.operations[0]);
-            signed_size_type const region1 = get_region_id(turn.operations[1]);
-
-            // Switch sources for same region
-            turn.switch_source = region0 == region1;
-        }
-
-
 #if defined(BOOST_GEOMETRY_DEBUG_TRAVERSAL_SWITCH_DETECTOR)
-        std::cout << "SWITCH END ITERATION" << std::endl;
+        std::cout << "END ITERATION GETTIN REGION_IDS" << std::endl;
 
         for (std::size_t turn_index = 0; turn_index < m_turns.size(); ++turn_index)
         {
@@ -646,25 +630,19 @@ struct traversal_switch_detector
             if ((turn.both(operation_union) || turn.both(operation_intersection))
                  && ! turn.is_clustered())
             {
-                std::cout << "UU/II SWITCH RESULT "
+                std::cout << "UU/II RESULT "
                              << turn_index << " -> "
-                          << turn.switch_source << std::endl;
+                          << turn.operations[0].enriched.region_id
+                          << " " << turn.operations[1].enriched.region_id
+                          << std::endl;
             }
         }
 
         for (typename Clusters::const_iterator it = m_clusters.begin(); it != m_clusters.end(); ++it)
         {
             cluster_info const& cinfo = it->second;
-            if (cinfo.open_count > 1)
-            {
-                std::cout << "CL SWITCH RESULT " << it->first
-                             << " -> " << cinfo.switch_source << std::endl;
-            }
-            else
-            {
-                std::cout << "CL SWITCH RESULT " << it->first
-                          << " is not registered as open" << std::endl;
-            }
+            std::cout << "CL RESULT " << it->first
+                         << " -> " << cinfo.open_count << std::endl;
         }
 #endif
 
diff --git a/boost/geometry/algorithms/detail/overlay/turn_info.hpp b/boost/geometry/algorithms/detail/overlay/turn_info.hpp
index 3ed75cf09e..545b5e902c 100644
--- a/boost/geometry/algorithms/detail/overlay/turn_info.hpp
+++ b/boost/geometry/algorithms/detail/overlay/turn_info.hpp
@@ -94,7 +94,6 @@ struct turn_info
     bool discarded;
 
     bool has_colocated_both; // Colocated with a uu turn (for union) or ii (other)
-    bool switch_source; // For u/u turns which can either switch or not
 
     Container operations;
 
@@ -104,7 +103,6 @@ struct turn_info
         , cluster_id(-1)
         , discarded(false)
         , has_colocated_both(false)
-        , switch_source(false)
     {}
 
     inline bool both(operation_type type) const
diff --git a/boost/geometry/algorithms/detail/partition.hpp b/boost/geometry/algorithms/detail/partition.hpp
index db134d548d..8eb0c7b1fd 100644
--- a/boost/geometry/algorithms/detail/partition.hpp
+++ b/boost/geometry/algorithms/detail/partition.hpp
@@ -765,7 +765,7 @@ public:
                              std::size_t min_elements)
     {
         return apply(forward_range1, forward_range2, visitor,
-                     expand_policy1, overlaps_policy1, expand_policy2, overlaps_policy1,
+                     expand_policy1, overlaps_policy1, expand_policy2, overlaps_policy2,
                      min_elements, detail::partition::visit_no_policy());
     }
 
diff --git a/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp b/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp
index b8ea5e30e6..ccd3af92d5 100644
--- a/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp
+++ b/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp
@@ -37,8 +37,8 @@ namespace detail
 
 template <typename Point1, typename Point2, typename Point3>
 inline bool collinear_point_is_spike_or_equal(Point1 const& last_point,
-                                                Point2 const& segment_a,
-                                                Point3 const& segment_b)
+                                              Point2 const& segment_a,
+                                              Point3 const& segment_b)
 {
     // Check if segment is equal
     int const sgn_x1 = sign_of_difference<0>(last_point, segment_b);
@@ -70,10 +70,10 @@ template
     typename Point1, typename Point2, typename Point3,
     typename SideStrategy
 >
-static inline bool point_is_spike_or_equal(Point1 const& last_point, // prev | back
-                                           Point2 const& segment_a,  // next | back - 2
-                                           Point3 const& segment_b,  // curr | back - 1 | spike's vertex
-                                           SideStrategy const& strategy)
+inline bool point_is_spike_or_equal(Point1 const& last_point, // prev | back
+                                    Point2 const& segment_a,  // next | back - 2
+                                    Point3 const& segment_b,  // curr | back - 1 | spike's vertex
+                                    SideStrategy const& strategy)
 {
     int const side = strategy.apply(segment_a, segment_b, last_point);
     if (side == 0)
@@ -100,7 +100,7 @@ template
     typename SideStrategy,
     typename RobustPolicy
 >
-static inline bool point_is_spike_or_equal(Point1 const& last_point,
+inline bool point_is_spike_or_equal(Point1 const& last_point,
             Point2 const& segment_a,
             Point3 const& segment_b,
             SideStrategy const& strategy,
@@ -137,6 +137,48 @@ static inline bool point_is_spike_or_equal(Point1 const& last_point,
         );
 }
 
+template
+<
+    typename Point1,
+    typename Point2,
+    typename Point3,
+    typename SideStrategy,
+    typename RobustPolicy
+>
+inline bool point_is_collinear(Point1 const& last_point,
+            Point2 const& segment_a,
+            Point3 const& segment_b,
+            SideStrategy const& strategy,
+            RobustPolicy const& robust_policy)
+{
+    int const side = strategy.apply(segment_a, segment_b, last_point);
+    if (side == 0)
+    {
+        return true;
+    }
+
+    // This part (or whole method, because it is then trivial)
+    // will be removed after rescaling
+    if (BOOST_GEOMETRY_CONDITION(! RobustPolicy::enabled))
+    {
+        return false;
+    }
+
+    // Redo, using specified robust policy
+    typedef typename geometry::robust_point_type
+    <
+        Point1,
+        RobustPolicy
+    >::type robust_point_type;
+
+    robust_point_type last_point_rob, segment_a_rob, segment_b_rob;
+    geometry::recalculate(last_point_rob, last_point, robust_policy);
+    geometry::recalculate(segment_a_rob, segment_a, robust_policy);
+    geometry::recalculate(segment_b_rob, segment_b, robust_policy);
+
+    int const side_rob = strategy.apply(segment_a_rob, segment_b_rob, last_point_rob);
+    return side_rob == 0;
+}
 
 } // namespace detail
 #endif
diff --git a/boost/geometry/algorithms/length.hpp b/boost/geometry/algorithms/length.hpp
index 39a567a26b..cfc46a07bd 100644
--- a/boost/geometry/algorithms/length.hpp
+++ b/boost/geometry/algorithms/length.hpp
@@ -51,6 +51,7 @@
 #include <boost/geometry/algorithms/detail/multi_sum.hpp>
 // #include <boost/geometry/algorithms/detail/throw_on_empty_input.hpp>
 #include <boost/geometry/views/closeable_view.hpp>
+#include <boost/geometry/strategies/default_strategy.hpp>
 #include <boost/geometry/strategies/distance.hpp>
 #include <boost/geometry/strategies/default_length_result.hpp>
 
@@ -184,6 +185,33 @@ struct length<MultiLinestring, multi_linestring_tag> : detail::multi_sum
 #endif // DOXYGEN_NO_DISPATCH
 
 
+namespace resolve_strategy {
+
+struct length
+{
+    template <typename Geometry, typename Strategy>
+    static inline typename default_length_result<Geometry>::type
+    apply(Geometry const& geometry, Strategy const& strategy)
+    {
+        return dispatch::length<Geometry>::apply(geometry, strategy);
+    }
+
+    template <typename Geometry>
+    static inline typename default_length_result<Geometry>::type
+    apply(Geometry const& geometry, default_strategy)
+    {
+        typedef typename strategy::distance::services::default_strategy
+            <
+                point_tag, point_tag, typename point_type<Geometry>::type
+            >::type strategy_type;
+
+        return dispatch::length<Geometry>::apply(geometry, strategy_type());
+    }
+};
+
+} // namespace resolve_strategy
+
+
 namespace resolve_variant {
 
 template <typename Geometry>
@@ -193,7 +221,7 @@ struct length
     static inline typename default_length_result<Geometry>::type
     apply(Geometry const& geometry, Strategy const& strategy)
     {
-        return dispatch::length<Geometry>::apply(geometry, strategy);
+        return resolve_strategy::length::apply(geometry, strategy);
     }
 };
 
@@ -255,13 +283,7 @@ length(Geometry const& geometry)
 
     // detail::throw_on_empty_input(geometry);
 
-    // TODO put this into a resolve_strategy stage
-    typedef typename strategy::distance::services::default_strategy
-        <
-            point_tag, point_tag, typename point_type<Geometry>::type
-        >::type strategy_type;
-
-    return resolve_variant::length<Geometry>::apply(geometry, strategy_type());
+    return resolve_variant::length<Geometry>::apply(geometry, default_strategy());
 }
 
 
diff --git a/boost/geometry/algorithms/point_on_surface.hpp b/boost/geometry/algorithms/point_on_surface.hpp
index 3f4d0f4afe..3bd18d83dd 100644
--- a/boost/geometry/algorithms/point_on_surface.hpp
+++ b/boost/geometry/algorithms/point_on_surface.hpp
@@ -3,7 +3,7 @@
 // Copyright (c) 2007-2013 Barend Gehrels, Amsterdam, the Netherlands.
 // Copyright (c) 2008-2013 Bruno Lalande, Paris, France.
 // Copyright (c) 2009-2013 Mateusz Loskot, London, UK.
-// Copyright (c) 2013 Adam Wulkiewicz, Lodz, Poland.
+// Copyright (c) 2013-2017 Adam Wulkiewicz, Lodz, Poland.
 
 // This file was modified by Oracle on 2014, 2017.
 // Modifications copyright (c) 2014-2017 Oracle and/or its affiliates.
@@ -31,6 +31,7 @@
 #include <boost/geometry/geometries/concepts/check.hpp>
 
 #include <boost/geometry/algorithms/detail/extreme_points.hpp>
+#include <boost/geometry/algorithms/detail/signed_size_type.hpp>
 
 #include <boost/geometry/strategies/cartesian/centroid_bashein_detmer.hpp>
 #include <boost/geometry/strategies/side.hpp>
@@ -154,7 +155,6 @@ inline void calculate_average(Point& point, std::vector<P> const& points)
 {
     typedef typename geometry::coordinate_type<Point>::type coordinate_type;
     typedef typename std::vector<P>::const_iterator iterator_type;
-    typedef typename std::vector<P>::size_type size_type;
 
     coordinate_type x = 0;
     coordinate_type y = 0;
@@ -166,7 +166,7 @@ inline void calculate_average(Point& point, std::vector<P> const& points)
         y += geometry::get<1>(*it);
     }
 
-    size_type const count = points.size();
+    signed_size_type const count = points.size();
     geometry::set<0>(point, x / count);
     geometry::set<1>(point, y / count);
 }
diff --git a/boost/geometry/algorithms/simplify.hpp b/boost/geometry/algorithms/simplify.hpp
index cfeb542220..298f9c640a 100644
--- a/boost/geometry/algorithms/simplify.hpp
+++ b/boost/geometry/algorithms/simplify.hpp
@@ -36,9 +36,13 @@
 #include <boost/geometry/strategies/default_strategy.hpp>
 #include <boost/geometry/strategies/distance.hpp>
 
+#include <boost/geometry/algorithms/area.hpp>
 #include <boost/geometry/algorithms/clear.hpp>
 #include <boost/geometry/algorithms/convert.hpp>
+#include <boost/geometry/algorithms/detail/equals/point_point.hpp>
 #include <boost/geometry/algorithms/not_implemented.hpp>
+#include <boost/geometry/algorithms/is_empty.hpp>
+#include <boost/geometry/algorithms/perimeter.hpp>
 
 #include <boost/geometry/algorithms/detail/distance/default_strategies.hpp>
 
@@ -49,6 +53,14 @@ namespace boost { namespace geometry
 namespace detail { namespace simplify
 {
 
+template <typename Range>
+inline bool is_degenerate(Range const& range)
+{
+    return boost::size(range) == 2
+        && detail::equals::equals_point_point(geometry::range::front(range),
+                                              geometry::range::back(range));
+}
+
 struct simplify_range_insert
 {
     template<typename Range, typename Strategy, typename OutputIterator, typename Distance>
@@ -57,7 +69,11 @@ struct simplify_range_insert
     {
         boost::ignore_unused(strategy);
 
-        if (boost::size(range) <= 2 || max_distance < 0)
+        if (is_degenerate(range))
+        {
+            std::copy(boost::begin(range), boost::begin(range) + 1, out);
+        }
+        else if (boost::size(range) <= 2 || max_distance < 0)
         {
             std::copy(boost::begin(range), boost::end(range), out);
         }
@@ -71,41 +87,32 @@ struct simplify_range_insert
 
 struct simplify_copy
 {
-    template <typename Range, typename Strategy, typename Distance>
-    static inline void apply(Range const& range, Range& out,
+    template <typename RangeIn, typename RangeOut, typename Strategy, typename Distance>
+    static inline void apply(RangeIn const& range, RangeOut& out,
                              Distance const& , Strategy const& )
     {
         std::copy
             (
-                boost::begin(range), boost::end(range), geometry::range::back_inserter(out)
+                boost::begin(range), boost::end(range),
+                    geometry::range::back_inserter(out)
             );
     }
 };
 
 
-template<std::size_t Minimum>
+template <std::size_t MinimumToUseStrategy>
 struct simplify_range
 {
-    template <typename Range, typename Strategy, typename Distance>
-    static inline void apply(Range const& range, Range& out,
+    template <typename RangeIn, typename RangeOut, typename Strategy, typename Distance>
+    static inline void apply(RangeIn const& range, RangeOut& out,
                     Distance const& max_distance, Strategy const& strategy)
     {
-        // Call do_container for a linestring / ring
-
-        /* For a RING:
-            The first/last point (the closing point of the ring) should maybe
-            be excluded because it lies on a line with second/one but last.
-            Here it is never excluded.
+        // For a RING:
+        // Note that, especially if max_distance is too large,
+        // the output ring might be self intersecting while the input ring is
+        // not, although chances are low in normal polygons
 
-            Note also that, especially if max_distance is too large,
-            the output ring might be self intersecting while the input ring is
-            not, although chances are low in normal polygons
-
-            Finally the inputring might have 3 (open) or 4 (closed) points (=correct),
-                the output < 3 or 4(=wrong)
-        */
-
-        if (boost::size(range) <= int(Minimum) || max_distance < 0.0)
+        if (boost::size(range) <= MinimumToUseStrategy || max_distance < 0)
         {
             simplify_copy::apply(range, out, max_distance, strategy);
         }
@@ -116,34 +123,172 @@ struct simplify_range
                     range, geometry::range::back_inserter(out), max_distance, strategy
                 );
         }
+
+        // Verify the two remaining points are equal. If so, remove one of them.
+        // This can cause the output being under the minimum size
+        if (is_degenerate(out))
+        {
+            range::resize(out, 1);
+        }
+    }
+};
+
+struct simplify_ring
+{
+private :
+    template <typename Area>
+    static inline int area_sign(Area const& area)
+    {
+        return area > 0 ? 1 : area < 0 ? -1 : 0;
+    }
+
+    template <typename Strategy, typename Ring>
+    static std::size_t get_opposite(std::size_t index, Ring const& ring)
+    {
+        typename Strategy::distance_strategy_type distance_strategy;
+
+        // Verify if it is NOT the case that all points are less than the
+        // simplifying distance. If so, output is empty.
+        typename Strategy::distance_type max_distance(-1);
+
+        typename geometry::point_type<Ring>::type point = range::at(ring, index);
+        std::size_t i = 0;
+        for (typename boost::range_iterator<Ring const>::type
+                it = boost::begin(ring); it != boost::end(ring); ++it, ++i)
+        {
+            // This actually is point-segment distance but will result
+            // in point-point distance
+            typename Strategy::distance_type dist = distance_strategy.apply(*it, point, point);
+            if (dist > max_distance)
+            {
+                max_distance = dist;
+                index = i;
+            }
+        }
+        return index;
+    }
+
+public :
+    template <typename Ring, typename Strategy, typename Distance>
+    static inline void apply(Ring const& ring, Ring& out,
+                    Distance const& max_distance, Strategy const& strategy)
+    {
+        std::size_t const size = boost::size(ring);
+        if (size == 0)
+        {
+            return;
+        }
+
+        int const input_sign = area_sign(geometry::area(ring));
+
+        std::set<std::size_t> visited_indexes;
+
+        // Rotate it into a copied vector
+        // (vector, because source type might not support rotation)
+        // (duplicate end point will be simplified away)
+        typedef typename geometry::point_type<Ring>::type point_type;
+
+        std::vector<point_type> rotated(size);
+
+        // Closing point (but it will not start here)
+        std::size_t index = 0;
+
+        // Iterate (usually one iteration is enough)
+        for (std::size_t iteration = 0; iteration < 4u; iteration++)
+        {
+            // Always take the opposite. Opposite guarantees that no point
+            // "halfway" is chosen, creating an artefact (very narrow triangle)
+            // Iteration 0: opposite to closing point (1/2, = on convex hull)
+            //              (this will start simplification with that point
+            //               and its opposite ~0)
+            // Iteration 1: move a quarter on that ring, then opposite to 1/4
+            //              (with its opposite 3/4)
+            // Iteration 2: move an eight on that ring, then opposite (1/8)
+            // Iteration 3: again move a quarter, then opposite (7/8)
+            // So finally 8 "sides" of the ring have been examined (if it were
+            // a semi-circle). Most probably, there are only 0 or 1 iterations.
+            switch (iteration)
+            {
+                case 1 : index = (index + size / 4) % size; break;
+                case 2 : index = (index + size / 8) % size; break;
+                case 3 : index = (index + size / 4) % size; break;
+            }
+            index = get_opposite<Strategy>(index, ring);
+
+            if (visited_indexes.count(index) > 0)
+            {
+                // Avoid trying the same starting point more than once
+                continue;
+            }
+
+            std::rotate_copy(boost::begin(ring), range::pos(ring, index),
+                             boost::end(ring), rotated.begin());
+
+            // Close the rotated copy
+            rotated.push_back(range::at(ring, index));
+
+            simplify_range<0>::apply(rotated, out, max_distance, strategy);
+
+            // Verify that what was positive, stays positive (or goes to 0)
+            // and what was negative stays negative (or goes to 0)
+            int const output_sign = area_sign(geometry::area(out));
+            if (output_sign == input_sign)
+            {
+                // Result is considered as satisfactory (usually this is the
+                // first iteration - only for small rings, having a scale
+                // similar to simplify_distance, next iterations are tried
+                return;
+            }
+
+            // Original is simplified away. Possibly there is a solution
+            // when another starting point is used
+            geometry::clear(out);
+
+            if (iteration == 0
+                && geometry::perimeter(ring) < 3 * max_distance)
+            {
+                // Check if it is useful to iterate. A minimal triangle has a
+                // perimeter of a bit more than 3 times the simplify distance
+                return;
+            }
+
+            // Prepare next try
+            visited_indexes.insert(index);
+            rotated.resize(size);
+        }
     }
 };
 
+
 struct simplify_polygon
 {
 private:
 
     template
     <
-        std::size_t Minimum,
         typename IteratorIn,
-        typename IteratorOut,
+        typename InteriorRingsOut,
         typename Distance,
         typename Strategy
     >
     static inline void iterate(IteratorIn begin, IteratorIn end,
-                    IteratorOut it_out,
+                    InteriorRingsOut& interior_rings_out,
                     Distance const& max_distance, Strategy const& strategy)
     {
-        for (IteratorIn it_in = begin; it_in != end;  ++it_in, ++it_out)
+        typedef typename boost::range_value<InteriorRingsOut>::type single_type;
+        for (IteratorIn it = begin; it != end; ++it)
         {
-            simplify_range<Minimum>::apply(*it_in, *it_out, max_distance, strategy);
+            single_type out;
+            simplify_ring::apply(*it, out, max_distance, strategy);
+            if (! geometry::is_empty(out))
+            {
+                range::push_back(interior_rings_out, out);
+            }
         }
     }
 
     template
     <
-        std::size_t Minimum,
         typename InteriorRingsIn,
         typename InteriorRingsOut,
         typename Distance,
@@ -154,12 +299,11 @@ private:
                     InteriorRingsOut& interior_rings_out,
                     Distance const& max_distance, Strategy const& strategy)
     {
-        traits::resize<InteriorRingsOut>::apply(interior_rings_out,
-            boost::size(interior_rings_in));
+        range::clear(interior_rings_out);
 
-        iterate<Minimum>(
+        iterate(
             boost::begin(interior_rings_in), boost::end(interior_rings_in),
-            boost::begin(interior_rings_out),
+            interior_rings_out,
             max_distance, strategy);
     }
 
@@ -168,21 +312,14 @@ public:
     static inline void apply(Polygon const& poly_in, Polygon& poly_out,
                     Distance const& max_distance, Strategy const& strategy)
     {
-        std::size_t const minimum = core_detail::closure::minimum_ring_size
-            <
-                geometry::closure<Polygon>::value
-            >::value;
-
         // Note that if there are inner rings, and distance is too large,
         // they might intersect with the outer ring in the output,
         // while it didn't in the input.
-        simplify_range<minimum>::apply(exterior_ring(poly_in),
-                                       exterior_ring(poly_out),
-                                       max_distance, strategy);
+        simplify_ring::apply(exterior_ring(poly_in), exterior_ring(poly_out),
+            max_distance, strategy);
 
-        apply_interior_rings<minimum>(interior_rings(poly_in),
-                                      interior_rings(poly_out),
-                                      max_distance, strategy);
+        apply_interior_rings(interior_rings(poly_in),
+            interior_rings(poly_out), max_distance, strategy);
     }
 };
 
@@ -194,16 +331,19 @@ struct simplify_multi
     static inline void apply(MultiGeometry const& multi, MultiGeometry& out,
                     Distance const& max_distance, Strategy const& strategy)
     {
-        traits::resize<MultiGeometry>::apply(out, boost::size(multi));
+        range::clear(out);
+
+        typedef typename boost::range_value<MultiGeometry>::type single_type;
 
-        typename boost::range_iterator<MultiGeometry>::type it_out
-                = boost::begin(out);
         for (typename boost::range_iterator<MultiGeometry const>::type
-                it_in = boost::begin(multi);
-             it_in != boost::end(multi);
-             ++it_in, ++it_out)
+                it = boost::begin(multi); it != boost::end(multi); ++it)
         {
-            Policy::apply(*it_in, *it_out, max_distance, strategy);
+            single_type single_out;
+            Policy::apply(*it, single_out, max_distance, strategy);
+            if (! geometry::is_empty(single_out))
+            {
+                range::push_back(out, single_out);
+            }
         }
     }
 };
@@ -236,7 +376,7 @@ struct simplify<Point, point_tag>
     }
 };
 
-
+// Linestring, keep 2 points (unless those points are the same)
 template <typename Linestring>
 struct simplify<Linestring, linestring_tag>
     : detail::simplify::simplify_range<2>
@@ -244,13 +384,7 @@ struct simplify<Linestring, linestring_tag>
 
 template <typename Ring>
 struct simplify<Ring, ring_tag>
-    : detail::simplify::simplify_range
-            <
-                core_detail::closure::minimum_ring_size
-                    <
-                        geometry::closure<Ring>::value
-                    >::value
-            >
+    : detail::simplify::simplify_ring
 {};
 
 template <typename Polygon>