From 132fdd87ee50d211382b32276328e1dc0dac9e86 Mon Sep 17 00:00:00 2001
From: Barend Gehrels <barend@xs4all.nl>
Date: Tue, 9 Feb 2010 09:24:02 +0000
Subject: [PATCH] Bugfix in assemble for multi Bugfix in get_turns in one
 constellation Other preparations for dissolve/buffer

[SVN r59593]
---
 .../detail/overlay/get_turn_info.hpp          |  19 +-
 .../detail/overlay/ring_properties.hpp        | 274 ++++++-------
 .../algorithms/detail/overlay/turn_info.hpp   |   4 +
 .../geometry/algorithms/intersection.hpp      |   2 +-
 .../geometry/algorithms/overlay/assemble.hpp  | 380 +++++++++---------
 .../geometry/algorithms/overlay/traverse.hpp  | 330 ++++++++-------
 include/boost/geometry/algorithms/union.hpp   |   2 +-
 .../multi/algorithms/overlay/assemble.hpp     |   4 +-
 8 files changed, 536 insertions(+), 479 deletions(-)

diff --git a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp
index 161506fd3..c61a53880 100644
--- a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp
+++ b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp
@@ -47,6 +47,10 @@ struct base_turn_handler
     {
         ti.operations[0].operation = op;
         ti.operations[1].operation = op;
+        if (op == operation_union)
+        {
+            ti.ignore = true;
+        }
     }
 
     // If condition, first union/second intersection, else vice versa
@@ -63,10 +67,7 @@ struct base_turn_handler
     template <typename TurnInfo>
     static inline void uu_else_ii(bool condition, TurnInfo& ti)
     {
-        ti.operations[0].operation = condition
-                    ? operation_union : operation_intersection;
-        ti.operations[1].operation = condition
-                    ? operation_union : operation_intersection;
+        both(ti, condition ? operation_union : operation_intersection);
     }
 };
 
@@ -241,7 +242,7 @@ struct touch : public base_turn_handler
             {
                 // Collinear -> lines join, continue
                 // (#BRL2)
-                if (side_pk_q2 == 0)
+                if (side_pk_q2 == 0 && ! block_q)
                 {
                     both(ti, operation_continue);
                     return;
@@ -292,7 +293,7 @@ struct touch : public base_turn_handler
             else
             {
                 // Pk at other side than Qi/Pk
-                int const side_qk_q  = SideStrategy::apply(qi, qj, qk);
+                int const side_qk_q = SideStrategy::apply(qi, qj, qk);
                 bool const q_turns_left = side_qk_q == 1;
 
                 ti.operations[0].operation = q_turns_left
@@ -304,6 +305,12 @@ struct touch : public base_turn_handler
                             ? operation_union
                             : operation_intersection;
 
+                if (ti.operations[0].operation == operation_union
+                    && ti.operations[1].operation == operation_union)
+                {
+                    ti.ignore = true;
+                }
+
                 return;
             }
         }
diff --git a/include/boost/geometry/algorithms/detail/overlay/ring_properties.hpp b/include/boost/geometry/algorithms/detail/overlay/ring_properties.hpp
index e4a2457d3..abcf69ff7 100644
--- a/include/boost/geometry/algorithms/detail/overlay/ring_properties.hpp
+++ b/include/boost/geometry/algorithms/detail/overlay/ring_properties.hpp
@@ -8,12 +8,6 @@
 #ifndef BOOST_GEOMETRY_ALG_DET_OV_RING_PROPERTIES_HPP
 #define BOOST_GEOMETRY_ALG_DET_OV_RING_PROPERTIES_HPP
 
-#define BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-#include <deque>
-#else
-#include <boost/array.hpp>
-#endif
 
 #include <boost/geometry/algorithms/area.hpp>
 #include <boost/geometry/algorithms/envelope.hpp>
@@ -31,89 +25,132 @@ namespace boost { namespace geometry
 template <typename Point>
 struct ring_properties
 {
-    struct parent
-    {
-        int c;
-        ring_identifier id;
+    typedef Point point_type;
+    typedef geometry::box<Point> box_type;
 
-        inline parent()
-            : c(0)
-        {}
+    ring_identifier ring_id;
+    typename area_result<Point>::type area;
+    int signum;
 
-        inline parent(int c_, ring_identifier i)
-            : c(c_), id(i)
-        {}
-    };
+    bool intersects;
+    bool produced;
 
-#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-    inline void push(parent const& p)
+    // "Stack"/counter of non-intersecting parent rings.
+    // This denotes if it a negative ring should be included,
+    int parent_count;
+
+    // ID of the parent
+    ring_identifier parent_ring_id;
+
+    box_type box;
+
+    point_type point;
+    bool has_point;
+
+    // Default constructor (necessary for vector, but not called)
+    inline ring_properties()
+        : intersects(false)
+        , produced(false)
+        , parent_count(0)
+        , has_point(false)
     {
-        parents.push_back(p);
+        parent_ring_id.source_index = -1;
     }
-    inline void pop()
+
+    template <typename Geometry>
+    inline ring_properties(ring_identifier id, Geometry const& geometry,
+            bool i, bool p = false)
+        : ring_id(id)
+        , area(geometry::area(geometry))
+        , intersects(i)
+        , produced(p)
+        , parent_count(0)
+        , box(geometry::make_envelope<box_type>(geometry))
     {
-        parents.pop_back();
+        has_point = geometry::point_on_border(point, geometry, true);
+        typedef typename coordinate_type<Geometry>::type coordinate_type;
+        coordinate_type zero = coordinate_type();
+        signum = area > zero ? 1 : area < zero ? -1 : 0;
+        parent_ring_id.source_index = -1;
     }
-    inline bool parent_stack_empty() const
+
+    inline bool positive() const { return signum == 1; }
+    inline bool negative() const { return signum == -1; }
+
+    inline bool operator<(ring_properties<Point> const& other) const
     {
-        return parents.empty();
+        // Normal sorting: in reverse order
+        return abs(area) > abs(other.area);
     }
-#else
-    inline void push(parent const& p)
+
+    inline ring_identifier const& id(int direction) const
     {
-        parents[array_size++] = p;
+        // Return the id of ifself, or of the parent
+        return positive() || parent_ring_id.source_index < 0
+            ? ring_id
+            : parent_ring_id;
     }
-    inline void pop()
+
+
+    inline void push(ring_properties<Point> const& r,
+                int direction, bool dissolve)
     {
-        array_size--;
-    }
-    inline bool parent_stack_empty() const
-    {
-        return array_size == 0;
-    }
+        if (//(r.produced || r.untouched()) &&
+            r.included(direction, dissolve))
+        {
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << " id.push " << r.ring_id;
 #endif
+            parent_ring_id = r.ring_id;
+        }
+        if (! r.produced || dissolve)
+        {
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << " or.push " << r.ring_id;
+#endif
+            parent_count++;
+        }
+    }
+
+
+    inline void pop(ring_properties<Point> const& r)
+    {
+        if (! r.produced && parent_count > 0)
+        {
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << " or.pop";
+#endif
+
+            parent_count--;
+        }
+    }
 
 
     inline bool interior_included(int direction) const
     {
-        if (area < 0 && ! parent_stack_empty())
+        if (negative())
         {
-            // Inner rings are included if the last encountered parent
-            // matches the operation
-#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-            int d = parents.back().c;
-            if (d == 0
-                && parents.size() == 2
-                && parents.front().c == 0)
-#else
-            int d = parents[array_size - 1].c;
-            if (d == 0
-                && array_size == 2
-                && parents[0].c == 0)
-#endif
+            // Original inner rings are included if there
+            // are two untouched parents (union) or one (intersection);
+            // Produced are ones are included if there is a parent found
+            if (produced)
             {
-                // It is contained by two outer rings, both of them
-                // are not included in the other. So they must be
-                // equal. Then always included, both in union and
-                // in intersection.
-                return direction == -1;
+                return parent_count > 0;
             }
-
-            if (d == 0)
-            {
-                d = 1;
-            }
-            return d * direction == 1;
+            return
+                (direction == 1  && parent_count == 1)
+                || (direction == -1 && parent_count > 1);
         }
         return false;
     }
 
-    inline bool included(int direction) const
+    inline bool included(int direction, bool dissolve) const
     {
-        if (produced)
+        if (produced && ! dissolve)
         {
             // Traversed rings are included in all operations,
             // because traversal was direction-dependant.
+            // On dissolve, this is not the case.
             return true;
         }
         if (intersects)
@@ -123,13 +160,19 @@ struct ring_properties
             return false;
         }
 
-        if (area > 0)
+        if (positive())
         {
             // Outer rings are included if they don't have parents
             // (union) or have parents (intersection)
-            return (parent_stack_empty() ? 1 : -1) * direction == 1;
+            if (produced)
+            {
+                return parent_count == 0;
+            }
+            return
+                (direction == 1  && parent_count == 0)
+                || (direction == -1 && parent_count > 0);
         }
-        else if (area < 0 && ! parent_stack_empty())
+        else if (negative())
         {
             // Inner rings are included if the last encountered parent
             // matches the operation
@@ -138,86 +181,28 @@ struct ring_properties
         return false;
     }
 
-    typedef Point point_type;
-    typedef geometry::box<Point> box_type;
-
-    ring_identifier ring_id;
-    typename area_result<Point>::type area;
-
-    int c;
-    bool intersects;
-    bool produced;
-
-    // Stack of non-intersecting parent rings.
-    // This denotes if it a negative ring should be included,
-    // and which is the parent.
-#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-    std::deque<parent> parents;
-#else
-    parent parents[3];
-    int array_size;
-#endif
-
-    // If the parent is an intersecting ring, take
-    // this as the parent.
-    ring_identifier parent_id;
-    bool has_parent_id;
-
-
-    box_type box;
-
-    point_type point;
-    bool has_point;
-
-    inline ring_properties()
-        : c(0)
-        , intersects(false)
-        , produced(false)
-        , has_point(false)
-        , has_parent_id(false)
-#ifndef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-        , array_size(0)
-#endif
-    {}
-
-    template <typename Geometry>
-    inline ring_properties(ring_identifier id, Geometry const& geometry,
-            bool i, bool p = false)
-        : ring_id(id)
-        , area(geometry::area(geometry))
-        , c(0)
-        , intersects(i)
-        , produced(p)
-        , has_parent_id(false)
-        , box(geometry::make_envelope<box_type>(geometry))
-#ifndef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-        , array_size(0)
-#endif
+    inline bool untouched() const
     {
-        has_point = geometry::point_on_border(point, geometry);
+        // It should be in comparisons on parent/child if:
+        // it is produced
+        // it is not produced, and not intersecting
+        return ! produced && ! intersects;
     }
 
-    inline bool operator<(ring_properties<Point> const& other) const
-    {
-        // Normal sorting: in reverse order
-        return std::abs(area) > std::abs(other.area);
-    }
 
-    inline ring_identifier const& id() const
+#if defined(BOOST_GEOMETRY_DEBUG_IDENTIFIER)
+    friend std::ostream& operator<<(std::ostream &os, ring_properties<Point> const& prop)
     {
-        if (has_parent_id)
-        {
-            return parent_id;
-        }
-        return parent_stack_empty() || area > 0
-            ? ring_id
-#ifdef BOOST_GEOMETRY_ASSEMBLE_PARENT_DEQUE
-            : parents.back().id;
-#else
-            : parents[array_size - 1].id;
+        os << "prop: " << prop.ring_id << " " << prop.area;
+        os << " count: " << prop.parent_count;
+        std::cout << " parent: " << prop.parent_ring_id;
+        if (prop.produced) std::cout << " produced";
+        if (prop.intersects) std::cout << " intersects";
+        if (prop.included(1, false)) std::cout << " union";
+        if (prop.included(-1, false)) std::cout << " intersection";
+        return os;
+    }
 #endif
-    }
-
 
 };
 
@@ -225,14 +210,21 @@ struct ring_properties
 template<typename Prop>
 struct sort_on_id_or_parent_id
 {
+private :
+    int m_direction;
+public :
+    sort_on_id_or_parent_id(int direction)
+        : m_direction(direction)
+    {}
+
     inline bool operator()(Prop const& left, Prop const& right) const
     {
-        ring_identifier const& left_id = left.id();
-        ring_identifier const& right_id = right.id();
+        ring_identifier const& left_id = left.id(m_direction);
+        ring_identifier const& right_id = right.id(m_direction);
 
         // If it is the same, sort on size descending
         return left_id == right_id
-            ? std::abs(left.area) > std::abs(right.area)
+            ? abs(left.area) > abs(right.area)
             : left_id < right_id;
     }
 };
diff --git a/include/boost/geometry/algorithms/detail/overlay/turn_info.hpp b/include/boost/geometry/algorithms/detail/overlay/turn_info.hpp
index 3255fc9b7..6aa8b2de9 100644
--- a/include/boost/geometry/algorithms/detail/overlay/turn_info.hpp
+++ b/include/boost/geometry/algorithms/detail/overlay/turn_info.hpp
@@ -82,11 +82,15 @@ struct turn_info
 
     Point point;
     method_type method;
+    bool ignore;
+    bool rejected;
 
     Container operations;
 
     turn_info()
         : method(method_none)
+        , ignore(false)
+        , rejected(false)
     {}
 
 };
diff --git a/include/boost/geometry/algorithms/intersection.hpp b/include/boost/geometry/algorithms/intersection.hpp
index 41575bfe3..57ca814a5 100644
--- a/include/boost/geometry/algorithms/intersection.hpp
+++ b/include/boost/geometry/algorithms/intersection.hpp
@@ -90,7 +90,7 @@ template
     typename Strategy
 >
 struct intersection_inserter
-    : detail::overlay::overlay_and_assemble
+    : detail::overlay::overlay
         <G1, G2, OutputIterator, GeometryOut, -1, Strategy>
 {};
 
diff --git a/include/boost/geometry/algorithms/overlay/assemble.hpp b/include/boost/geometry/algorithms/overlay/assemble.hpp
index da28380a9..ce3569bdc 100644
--- a/include/boost/geometry/algorithms/overlay/assemble.hpp
+++ b/include/boost/geometry/algorithms/overlay/assemble.hpp
@@ -21,7 +21,6 @@
 
 #include <boost/geometry/algorithms/overlay/get_turns.hpp>
 #include <boost/geometry/algorithms/overlay/enrich_intersection_points.hpp>
-//#include <boost/geometry/algorithms/detail/overlay/enrich.hpp>
 #include <boost/geometry/algorithms/overlay/traverse.hpp>
 
 
@@ -37,8 +36,6 @@
 
 #include <boost/geometry/strategies/intersection.hpp>
 
-//#define BOOST_GEOMETRY_DEBUG_ASSEMBLE
-
 #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
 #  include <boost/geometry/util/write_dsv.hpp>
 #endif
@@ -165,10 +162,16 @@ struct add_to_containment<ring_tag, Ring>
 {
     template <typename ContainmentContainer, typename Map>
     static inline void apply(ContainmentContainer& container,
-            ring_identifier const& id, Ring const& ring, Map const& map)
+            ring_identifier const& id, Ring const& ring, Map const& map,
+            bool dissolve)
     {
         typedef typename range_value<ContainmentContainer>::type prop;
-        container.push_back(prop(id, ring, map.find(id) != map.end()));
+        bool found = map.find(id) != map.end();
+        if (! dissolve || ! found)
+        {
+            // For dissolve, do not add intersected rings
+            container.push_back(prop(id, ring, found));
+        }
     }
 };
 
@@ -177,10 +180,15 @@ struct add_to_containment<box_tag, Box>
 {
     template <typename ContainmentContainer, typename Map>
     static inline void apply(ContainmentContainer& container,
-            ring_identifier const& id, Box const& box, Map const& map)
+            ring_identifier const& id, Box const& box, Map const& map,
+            bool dissolve)
     {
         typedef typename range_value<ContainmentContainer>::type prop;
-        container.push_back(prop(id, box, map.find(id) != map.end()));
+        bool found = map.find(id) != map.end();
+        if (! dissolve || ! found)
+        {
+            container.push_back(prop(id, box, found));
+        }
     }
 };
 
@@ -190,7 +198,8 @@ struct add_to_containment<polygon_tag, Polygon>
 {
     template <typename ContainmentContainer, typename Map>
     static inline void apply(ContainmentContainer& container,
-            ring_identifier const& id, Polygon const& polygon, Map const& map)
+            ring_identifier const& id, Polygon const& polygon, Map const& map,
+            bool dissolve)
     {
         // Add exterior ring and interior rings
         ring_identifier copy = id;
@@ -201,7 +210,7 @@ struct add_to_containment<polygon_tag, Polygon>
                 typename ring_type<Polygon>::type
             > policy;
 
-        policy::apply(container, copy, exterior_ring(polygon), map);
+        policy::apply(container, copy, exterior_ring(polygon), map, dissolve);
         copy.ring_index = 0;
         for (typename boost::range_iterator
                 <
@@ -210,7 +219,7 @@ struct add_to_containment<polygon_tag, Polygon>
             it != boost::end(interior_rings(polygon));
             ++it, ++copy.ring_index)
         {
-            policy::apply(container, copy, *it, map);
+            policy::apply(container, copy, *it, map, dissolve);
         }
     }
 };
@@ -229,19 +238,22 @@ inline void map_turns(Map& map, TurnPoints const& turn_points)
          it != boost::end(turn_points);
          ++it, ++index)
     {
-        int op_index = 0;
-        for (typename boost::range_iterator<container_type const>::type
-                op_it = boost::begin(it->operations);
-            op_it != boost::end(it->operations);
-            ++op_it, ++op_index)
+        if (! it->ignore)
         {
-            ring_identifier ring_id
-                (
-                    op_it->seg_id.source_index,
-                    op_it->seg_id.multi_index,
-                    op_it->seg_id.ring_index
-                );
-            map[ring_id]++;
+            int op_index = 0;
+            for (typename boost::range_iterator<container_type const>::type
+                    op_it = boost::begin(it->operations);
+                op_it != boost::end(it->operations);
+                ++op_it, ++op_index)
+            {
+                ring_identifier ring_id
+                    (
+                        op_it->seg_id.source_index,
+                        op_it->seg_id.multi_index,
+                        op_it->seg_id.ring_index
+                    );
+                map[ring_id]++;
+            }
         }
     }
 }
@@ -259,46 +271,27 @@ struct enrich_containment
 {
 
     typedef typename boost::range_value<Container>::type item_type;
-    typedef typename item_type::parent parent_type;
     typedef typename geometry::tag<Geometry1>::type tag1;
     typedef typename geometry::tag<Geometry2>::type tag2;
     typedef void tag3; // For the ring-container
 
 
-    static inline void assign(item_type& item1, item_type& item2)
+    static inline void assign(item_type& larger, item_type& smaller, int direction, bool dissolve)
     {
-        int or1 = item1.area > 0 ? 1 : item1.area < 0 ? -1 : 0;
-        int or2 = item2.area > 0 ? 1 : item2.area < 0 ? -1 : 0;
-        if (or1 != 0 && or2 != 0)
-        {
-            if (or1 == 1 && or2 == 1)
-            {
-                item1.c = 1;
-                item2.c = -1;
-            }
+        typedef typename geometry::coordinate_type
+            <
+                Geometry1
+            >::type coordinate_type;
 
-            if (! item1.produced)
+        if (larger.signum != 0 && smaller.signum != 0)
+        {
+            if (larger.signum == 1)
             {
-                // This is an original ring, used to set info about
-                // the parents of interiors
-                if (or1 == 1)
-                {
-                    item2.push(parent_type(item1.c, item1.ring_id));
-                }
-                else if (or1 == -1 && ! item1.intersects)
-                {
-                    item2.pop();
-                }
+                smaller.push(larger, direction, dissolve);
             }
-            else
+            else if (larger.signum == -1)
             {
-                if (or1 == 1 && or2 == -1 && ! item2.intersects)
-                {
-                    // Assign this produced parent as the parent
-                    // to this hole
-                    item2.has_parent_id = true;
-                    item2.parent_id = item1.ring_id;
-                }
+                smaller.pop(/*from*/larger);
             }
         }
     }
@@ -348,7 +341,7 @@ struct enrich_containment
     static inline void enrich(Selection& selection, Map& map,
             Geometry1 const& geometry1, Geometry2 const& geometry2,
             RingCollection const& collection,
-            int direction)
+            int direction, bool dissolve)
     {
         typedef typename boost::range_iterator<Selection>::type iterator;
 
@@ -370,12 +363,12 @@ struct enrich_containment
 
                     item_type& item2 = **it2;
                     if (geometry::within(item2.point, item1.box)
-                        && std::abs(item2.area) < std::abs(item1.area)
+                        && abs(item2.area) < abs(item1.area)
                         && contains(item1, item2, geometry1, geometry2,
                                     collection)
                         )
                     {
-                        assign(item1, item2);
+                        assign(item1, item2, direction, dissolve);
                     }
                 }
             }
@@ -386,7 +379,7 @@ struct enrich_containment
     static inline void divide_and_conquer(Selection& selection, Map& map,
             Geometry1 const& geometry1, Geometry2 const& geometry2,
             RingCollection const& collection,
-            int direction, Box const& box,
+            int direction, bool dissolve, Box const& box,
             std::size_t iteration = 0, std::size_t previous_count = 0)
     {
         std::size_t n = boost::size(selection);
@@ -400,14 +393,15 @@ std::cout << "spatial divide n="
 
         if (iteration > 3)
         {
-            enrich(selection, map, geometry1, geometry2, collection, direction);
+            enrich(selection, map, geometry1, geometry2, collection, direction, dissolve);
             return;
         }
 
         // Divide the complete box in two (alternating) halves
         Box lower = box, upper = box;
-        typename geometry::coordinate_type<Box>::type two = 2.0;
-        typename geometry::coordinate_type<Box>::type mid
+        typedef typename geometry::coordinate_type<Box>::type coordinate_type;
+        coordinate_type two = 2.0;
+        coordinate_type mid
             = (geometry::get<min_corner, Dimension>(box)
                 + geometry::get<max_corner, Dimension>(box)) / two;
 
@@ -420,15 +414,15 @@ std::cout << "spatial divide n="
         select_by_box(upper_sel, selection, upper);
 
         divide_and_conquer<1 - Dimension>(lower_sel, map, geometry1, geometry2,
-                    collection, direction, lower, iteration + 1, n);
+                    collection, direction, dissolve, lower, iteration + 1, n);
         divide_and_conquer<1 - Dimension>(upper_sel, map, geometry1, geometry2,
-                    collection, direction, upper, iteration + 1, n);
+                    collection, direction, dissolve, upper, iteration + 1, n);
     }
 
     static inline void enrich(Container& container,
             Geometry1 const& geometry1, Geometry2 const& geometry2,
             RingCollection const& collection,
-            int direction)
+            int direction, bool dissolve)
     {
         typedef typename boost::range_iterator<Container>::type iterator;
 
@@ -439,17 +433,26 @@ std::cout << "spatial divide n="
             ++it1)
         {
             item_type& item1 = *it1;
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << item1.ring_id << " area: " << item1.area << std::endl;
+#endif
             iterator it2 = it1;
             for (it2++; it2 != boost::end(container); ++it2)
             {
                 item_type& item2 = *it2;
                 if (geometry::within(item2.point, item1.box)
-                    && std::abs(item2.area) < std::abs(item1.area)
+                    && abs(item2.area) < abs(item1.area)
                     && contains(item1, item2, geometry1, geometry2, collection)
                     )
                 {
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << " -> contains " << item2.ring_id;
+#endif
                     n++;
-                    assign(item1, item2);
+                    assign(item1, item2, direction, dissolve);
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << std::endl;
+#endif
                 }
             }
         }
@@ -461,13 +464,23 @@ std::cout << "spatial divide n="
     static inline void apply(Container& container,
             Geometry1 const& geometry1, Geometry2 const& geometry2,
             RingCollection const& collection,
-            int direction, Box const& box)
+            int direction, bool dissolve, Box const& box)
     {
         if (boost::size(container) == 0)
         {
             return;
         }
-        enrich(container, geometry1, geometry2, collection, direction);
+        enrich(container, geometry1, geometry2, collection, direction, dissolve);
+
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+        for (typename boost::range_iterator<Container const>::type
+                it = boost::begin(container);
+            it != boost::end(container);
+            ++it)
+        {
+            std::cout << *it << std::endl;
+        }
+#endif
         return;
 
         // Method using divide and conquer (does NOT work corretly!)
@@ -484,7 +497,7 @@ std::cout << "spatial divide n="
 
         std::map<std::pair<item_type*, item_type*>, bool> map;
         divide_and_conquer<1>(selection, map, geometry1, geometry2, collection,
-                    direction, box);
+                    direction, dissolve, box);
     }
 };
 
@@ -501,7 +514,7 @@ template
 inline OutputIterator add_all_rings(Container& container,
             Geometry1 const& geometry1, Geometry2 const& geometry2,
             RingCollection const& collection,
-            int direction,
+            int direction, bool dissolve,
             OutputIterator out)
 {
     typedef typename boost::range_iterator<Container>::type iterator;
@@ -520,7 +533,7 @@ inline OutputIterator add_all_rings(Container& container,
         it != boost::end(container);
         ++it)
     {
-        if (it->area > 0)
+        if (it->positive())
         {
             if (result_filled)
             {
@@ -531,7 +544,7 @@ inline OutputIterator add_all_rings(Container& container,
 
             // If it is an outer ring, it is included if there are no parents
             // (union) or if there are parents (intersection)
-            if (it->included(direction))
+            if (it->included(direction, dissolve))
             {
                 geometry::clear(result);
                 previous_id = it->ring_id;
@@ -558,8 +571,8 @@ inline OutputIterator add_all_rings(Container& container,
             // If it is an interior ring, it is included if
             // it's parent-id matches the id of the outputted exterior ring
             if (result_filled
-                && it->id() == previous_id
-                && it->included(direction))
+                && it->id(direction) == previous_id
+                && it->included(direction, dissolve))
             {
                 if (it->ring_id.source_index == 0)
                 {
@@ -589,140 +602,58 @@ inline OutputIterator add_all_rings(Container& container,
     return out;
 }
 
-template <typename Container>
-inline void sort_on_parent_id(Container& container, int direction)
-{
-    typedef typename boost::range_iterator<Container>::type iterator;
-    typedef typename boost::range_value<Container>::type item_type;
-
-    // For negative rings produced by traversal,
-    // find the parent of the selected operation
-    // (by popping from the stack if necessary)
-    for (iterator it = boost::begin(container);
-        it != boost::end(container);
-        ++it)
-    {
-        if (it->produced && it->area < 0)
-        {
-            while (! it->interior_included(direction)
-                && ! it->parent_stack_empty())
-            {
-                it->pop();
-            }
-        }
-    }
-
-    // Sort container on parent-id
-#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
-std::cout << "assemble.properties sort on parent-id "
-    << boost::size(container) << std::endl;
-#endif
-    std::sort(boost::begin(container), boost::end(container),
-            sort_on_id_or_parent_id<item_type>());
-}
 
 template
 <
+    typename GeometryOut,
+    typename Rings, typename Turns,
     typename Geometry1, typename Geometry2,
-    typename OutputIterator, typename GeometryOut,
-    int Direction,
-    typename Strategy
+    typename OutputIterator
 >
-struct overlay_and_assemble
+inline OutputIterator assemble(Rings const& rings, Turns& turn_points,
+            Geometry1 const& geometry1, 
+            Geometry2 const& geometry2,
+            int direction, bool dissolve,
+            OutputIterator out)
 {
-    typedef typename geometry::tag<Geometry1>::type tag1;
-    typedef typename geometry::tag<Geometry2>::type tag2;
-    typedef typename geometry::tag<GeometryOut>::type tag_out;
-
-    static inline OutputIterator apply(Geometry1 const& geometry1,
-                Geometry2 const& geometry2, OutputIterator out,
-                Strategy const& strategy)
-    {
-        if (geometry::num_points(geometry1) == 0
-            || geometry::num_points(geometry2) == 0)
-        {
-            return out;
-        }
+        typedef typename geometry::tag<Geometry1>::type tag1;
+        typedef typename geometry::tag<Geometry2>::type tag2;
+        typedef typename geometry::tag<GeometryOut>::type tag_out;
 
         typedef typename geometry::point_type<GeometryOut>::type point_type;
+        typedef typename boost::range_value<Rings>::type ring_type;
 
-
-        typedef detail::overlay::traversal_turn_info<point_type> turn_info;
-        typedef std::deque<turn_info> container_type;
-
-        // "Abuse" rangetype for ringtype: for polygon, it is the type of the
-        // exterior ring. For ring, it is the ring itself. That is what is
-        // wished here as well.
-        typedef typename geometry::range_type<GeometryOut>::type ring_type;
-
-        container_type turn_points;
-
-#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
-std::cout << "get turns" << std::endl;
-#endif
-        boost::geometry::get_turns
-            <
-                detail::overlay::CalculateDistancePolicy
-            >(geometry1, geometry2, turn_points);
-
-#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
-std::cout << "enrich" << std::endl;
-#endif
-        typename Strategy::side_strategy_type side_strategy;
-        geometry::enrich_intersection_points(turn_points, geometry1, geometry2,
-                    side_strategy);
-
-        std::vector<ring_type> rings;
-#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
-std::cout << "traverse" << std::endl;
-#endif
-        geometry::traverse<typename Strategy::side_strategy_type, ring_type>
-            (
-                geometry1,
-                geometry2,
-                Direction == -1
-                    ? boost::geometry::detail::overlay::operation_intersection
-                    : boost::geometry::detail::overlay::operation_union
-                    ,
-                turn_points,
-                std::back_inserter(rings)
-            );
 #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
 std::cout << "assemble" << std::endl;
 #endif
 
-
         // Map intersection-points per ring-identifier to <count>
         std::map<ring_identifier, int> map;
         map_turns(map, turn_points);
 
-// TODO: translate/implement this (sorry for the dutch)
-//1. loop ook door originele ringen heen, map[ring] = 0
-//2. verander dan map.find=.end in .find -> it=0
-//3. we weten nu dat er "untouched"  ringen zijn, als dat zo is,
-// dan "enrich assemble" etc.kl
-
         typedef std::vector
             <
                 ring_properties<point_type>
             > ring_properties_container_type;
         ring_properties_container_type ring_properties_container;
 
-
         add_to_containment
             <
                 tag1,
                 Geometry1
             >::apply(ring_properties_container,
-                        ring_identifier(0,-1,-1), geometry1, map);
-        add_to_containment
-            <
-                tag2,
-                Geometry2
-            >::apply(ring_properties_container,
-                        ring_identifier(1,-1,-1), geometry2, map);
-
-
+                        ring_identifier(0, -1,-1), geometry1, 
+                        map, dissolve);
+        if (! dissolve)
+        {
+            add_to_containment
+                <
+                    tag2,
+                    Geometry2
+                >::apply(ring_properties_container,
+                            ring_identifier(1, -1,-1), geometry2, 
+                            map, dissolve);
+        }
 
         // Add all produced rings using source index 2
         {
@@ -745,7 +676,7 @@ std::cout << "assemble" << std::endl;
                 it != boost::end(ring_properties_container);
                 ++it)
             {
-                if (it->area>0)
+                if (it->positive())
                 {
                     pos++;
                 }
@@ -788,15 +719,100 @@ std::cout << "assemble.enrich containment" << std::endl;
                     std::vector<ring_type>,
                     geometry::box<point_type>
                 >::apply(ring_properties_container,
-                        geometry1, geometry2, rings, Direction, total);
+                        geometry1, geometry2, rings, direction, dissolve, total);
 
-            sort_on_parent_id(ring_properties_container, Direction);
+    // Sort container on parent-id
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << "assemble.properties sort on parent-id "
+    << boost::size(ring_properties_container) << std::endl;
+#endif
+            std::sort(boost::begin(ring_properties_container), 
+                    boost::end(ring_properties_container),
+                    sort_on_id_or_parent_id
+                        <
+                            ring_properties<point_type> 
+                        >(direction));
         }
 #ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
 std::cout << "assemble.add rings" << std::endl;
 #endif
         return add_all_rings<GeometryOut>(ring_properties_container,
-                    geometry1, geometry2, rings, Direction, out);
+                    geometry1, geometry2, rings, direction, dissolve, out);
+}
+
+
+template
+<
+    typename Geometry1, typename Geometry2,
+    typename OutputIterator, typename GeometryOut,
+    int Direction,
+    typename Strategy
+>
+struct overlay
+{
+    static inline OutputIterator apply(
+                Geometry1 const& geometry1, Geometry2 const& geometry2, 
+                OutputIterator out,
+                Strategy const& strategy)
+    {
+        if (geometry::num_points(geometry1) == 0 && geometry::num_points(geometry2) == 0)
+        {
+            return out;
+        }
+
+        typedef typename geometry::point_type<GeometryOut>::type point_type;
+        typedef detail::overlay::traversal_turn_info<point_type> turn_info;
+        typedef std::deque<turn_info> container_type;
+
+        // "Use" rangetype for ringtype: 
+        // for polygon, it is the type of the exterior ring. 
+        // for ring, it is the ring itself. That is what is
+        // for multi-polygon, it is also the type of the ring.
+        typedef typename geometry::range_type<GeometryOut>::type ring_type;
+
+        container_type turn_points;
+        std::vector<ring_type> rings;
+
+        // If one input is empty, output the other one for a union. 
+        // For an intersection, the intersection is empty.
+        if (geometry::num_points(geometry1) == 0 
+            || geometry::num_points(geometry2) == 0)
+        {
+            if (Direction == 1)
+            {
+                return assemble<GeometryOut>(rings, turn_points, 
+                                geometry1, geometry2, Direction, false, out);
+            }
+            return out;
+        }
+
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << "get turns" << std::endl;
+#endif
+        boost::geometry::get_turns
+            <
+                detail::overlay::CalculateDistancePolicy
+            >(geometry1, geometry2, turn_points);
+
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << "enrich" << std::endl;
+#endif
+        typename Strategy::side_strategy_type side_strategy;
+        geometry::enrich_intersection_points(turn_points, geometry1, geometry2,
+                    side_strategy);
+
+#ifdef BOOST_GEOMETRY_DEBUG_ASSEMBLE
+std::cout << "traverse" << std::endl;
+#endif
+        geometry::traverse(geometry1, geometry2,
+                Direction == -1
+                    ? boost::geometry::detail::overlay::operation_intersection
+                    : boost::geometry::detail::overlay::operation_union
+                    ,
+                turn_points, rings);
+
+        return assemble<GeometryOut>(rings, turn_points, 
+                        geometry1, geometry2, Direction, false, out);
     }
 };
 
diff --git a/include/boost/geometry/algorithms/overlay/traverse.hpp b/include/boost/geometry/algorithms/overlay/traverse.hpp
index 928e3f09d..dbaf4d5af 100644
--- a/include/boost/geometry/algorithms/overlay/traverse.hpp
+++ b/include/boost/geometry/algorithms/overlay/traverse.hpp
@@ -18,29 +18,51 @@
 
 #include <boost/geometry/geometries/concepts/check.hpp>
 
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-#include <boost/geometry/extensions/gis/io/wkt/write_wkt.hpp>
+#if defined(BOOST_GEOMETRY_DEBUG_INTERSECTION) || defined(BOOST_GEOMETRY_OVERLAY_REPORT_WKT)
+#include <boost/geometry/extensions/gis/io/wkt/wkt.hpp>
 #endif
 
 
 #include <boost/geometry/algorithms/detail/overlay/turn_info.hpp>
 
 
-
-
 namespace boost { namespace geometry
 {
 
 
+
+
 #ifndef DOXYGEN_NO_DETAIL
 namespace detail { namespace overlay {
 
 
+template <typename Turns>
+inline void clear_visit_info(Turns& turns)
+{
+    typedef typename boost::range_value<Turns>::type tp_type;
+
+    for (typename boost::range_iterator<Turns>::type
+        it = boost::begin(turns);
+        it != boost::end(turns);
+        ++it)
+    {
+        for (typename boost::range_iterator
+            <
+                typename tp_type::container_type
+            >::type op_it = boost::begin(it->operations);
+            op_it != boost::end(it->operations);
+            ++op_it)
+        {
+            op_it->visited.set_none();
+        }
+    }
+}
+
+
 template <typename Info, typename Turn>
 inline void set_visited_for_contine(Info& info, Turn const& turn)
 {
     // On "continue", set "visited" for ALL directions
-
     if (turn.operation == detail::overlay::operation_continue)
     {
         for (typename boost::range_iterator
@@ -64,19 +86,16 @@ template
     typename GeometryOut,
     typename G1,
     typename G2,
-    typename TurnPoints,
+    typename Turns,
     typename IntersectionInfo
 >
 inline void assign_next_ip(G1 const& g1, G2 const& g2,
-            TurnPoints& turn_points,
-            typename boost::range_iterator<TurnPoints>::type & ip,
+            Turns& turns,
+            typename boost::range_iterator<Turns>::type & ip,
             GeometryOut& current_output,
             IntersectionInfo & info,
             segment_identifier& seg_id)
 {
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-    //std::cout << "traverse: take: " << info << std::endl;
-#endif
     info.visited.set_visited();
     set_visited_for_contine(*ip, info);
 
@@ -95,12 +114,12 @@ inline void assign_next_ip(G1 const& g1, G2 const& g2,
                     info.enriched.travels_to_vertex_index,
                     current_output);
         }
-        ip = boost::begin(turn_points) + info.enriched.travels_to_ip_index;
+        ip = boost::begin(turns) + info.enriched.travels_to_ip_index;
         seg_id = info.seg_id;
     }
     else
     {
-        ip = boost::begin(turn_points) + info.enriched.next_ip_index;
+        ip = boost::begin(turns) + info.enriched.next_ip_index;
         seg_id = info.seg_id;
     }
     *(std::back_inserter(current_output)++) = ip->point;
@@ -109,9 +128,7 @@ inline void assign_next_ip(G1 const& g1, G2 const& g2,
 
 inline bool select_source(operation_type operation, int source1, int source2)
 {
-    return
-        (operation == operation_intersection && source1 != source2)
-        || (operation == operation_union && source1 == source2);
+    return operation == operation_intersection && source1 != source2;
 }
 
 
@@ -134,7 +151,7 @@ inline bool select_next_ip(operation_type operation,
     {
         // In some cases there are two alternatives.
         // For "ii", take the other one (alternate)
-        // For "uu", take the same one
+        // For "uu", take the same one (see above);
         // For "cc", take either one, but if there is a starting one,
         //           take that one.
         if (   (it->operation == operation_continue
@@ -157,34 +174,50 @@ inline bool select_next_ip(operation_type operation,
         if (it->visited.started())
         {
             selected = it;
-            has_tp = true;
             return true;
         }
-
     }
 
     return has_tp;
 }
 
 
-template <typename Container>
-inline void stop_gracefully(Container& container, bool& stop,
-            std::string const& reason)
+// This backtracking approach is necessary for invalid (== self-interescting)
+// geometries and the approach will be enhanced. TODO
+template
+<
+    typename Rings,
+    typename Turns,
+    typename Geometry1,
+    typename Geometry2
+>
+inline void backtrack(std::size_t size_at_start, bool& fail,
+            Rings& rings, typename boost::range_value<Rings>::type& ring,
+            Turns& turns, typename boost::range_value<Turns>::type& turn,
+            std::string const& reason,
+            Geometry1 const& geometry1,
+            Geometry2 const& geometry2)
 {
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-    std::cout << "STOPPING: " << reason << std::endl;
+    fail = true;
+
+    // Make bad output clean
+    rings.resize(size_at_start);
+    ring.clear();
+
+    // Reject this as a starting point
+    turn.rejected = true;
+
+    // And clear all visit info
+    clear_visit_info(turns);
+
+
+#ifdef BOOST_GEOMETRY_OVERLAY_REPORT_WKT
+    std::cout << " BT (" << reason << " )";
+    std::cout
+        << geometry::wkt(geometry1) << std::endl
+        << geometry::wkt(geometry2) << std::endl;
 #endif
 
-    stop = true;
-    if (container.size() > 0)
-    {
-        // Two-step adding first value, without assignment
-        // and without (erroneous) push_back of reference
-        typename boost::range_value<Container>::type p;
-        geometry::copy_coordinates(container.front(), p);
-        *(std::back_inserter(container)++) = p;
-
-    }
 }
 
 }} // namespace detail::overlay
@@ -197,146 +230,151 @@ inline void stop_gracefully(Container& container, bool& stop,
  */
 template
 <
-    typename SideStrategy,
-    typename GeometryOut,
     typename Geometry1,
     typename Geometry2,
-    typename TurnPoints,
-    typename OutputIterator
+    typename Turns,
+    typename Rings
 >
 inline void traverse(Geometry1 const& geometry1,
             Geometry2 const& geometry2,
             detail::overlay::operation_type operation,
-            TurnPoints& turn_points,
-            OutputIterator out)
+            Turns& turns, Rings& rings)
 {
-    concept::check<const Geometry1>();
-    concept::check<const Geometry2>();
-
-    typedef typename boost::range_iterator<TurnPoints>::type turn_iterator;
-    typedef typename boost::range_value<TurnPoints>::type turn_type;
+    typedef typename boost::range_iterator<Turns>::type turn_iterator;
+    typedef typename boost::range_value<Turns>::type turn_type;
     typedef typename boost::range_iterator
         <
             typename turn_type::container_type
         >::type turn_operation_iterator_type;
 
+    std::size_t size_at_start = boost::size(rings);
 
-    // Iterate through all unvisited points
-    for (turn_iterator it = boost::begin(turn_points);
-        it != boost::end(turn_points);
-        ++it)
+    bool fail = false;
+    do
     {
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-        //std::cout << "traversal: try: " << *it;
-#endif
-
-        for (turn_operation_iterator_type iit = boost::begin(it->operations);
-            iit != boost::end(it->operations);
-            ++iit)
+        fail = false;
+        // Iterate through all unvisited points
+        for (turn_iterator it = boost::begin(turns);
+            ! fail && it != boost::end(turns);
+            ++it)
         {
-            if (iit->visited.none()
-                && (iit->operation == operation
-                    || iit->operation == detail::overlay::operation_continue)
-                )
+            // Skip the ones marked ignore (these are: "uu", so self-tangent)
+            if (! it->ignore)
             {
-                set_visited_for_contine(*it, *iit);
-
-                GeometryOut current_output;
-                *(std::back_inserter(current_output)++) = it->point;
-
-                turn_iterator current = it;
-                turn_operation_iterator_type current_iit = iit;
-                segment_identifier current_seg_id;
-
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-                //std::cout << "traversal: start: " << *current;
-#endif
-
-                detail::overlay::assign_next_ip(geometry1, geometry2,
-                            turn_points,
-                            current, current_output,
-                            *iit, current_seg_id);
-                detail::overlay::select_next_ip(
-                                operation,
-                                *current,
-                                current_seg_id,
-                                boost::end(it->operations),
-                                current_iit);
-
-                iit->visited.set_started();
-
-
-                unsigned int i = 0;
-                bool stop = false;
-
-                while (current_iit != iit && ! stop)
+                for (turn_operation_iterator_type iit = boost::begin(it->operations);
+                    ! fail && ! it->rejected && iit != boost::end(it->operations);
+                    ++iit)
                 {
-
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-                    //std::cout << "traverse: " << *current;
-#endif
-
-                    if (current_iit->visited.visited())
+                    if (iit->visited.none()
+                        && (iit->operation == operation
+                            || iit->operation == detail::overlay::operation_continue)
+                        )
                     {
-                        // It visits a visited node again, without passing the start node.
-                        // This makes it suspicious for endless loops
-                        // Check if it is really same node
-                        detail::overlay::stop_gracefully(
-                            current_output, stop, "Visit again");
-                    }
+                        set_visited_for_contine(*it, *iit);
+
+                        typename boost::range_value<Rings>::type current_output;
+                        *(std::back_inserter(current_output)++) = it->point;
+
+                        turn_iterator current = it;
+                        turn_operation_iterator_type current_iit = iit;
+                        segment_identifier current_seg_id;
+
+                        detail::overlay::assign_next_ip(geometry1, geometry2,
+                                    turns,
+                                    current, current_output,
+                                    *iit, current_seg_id);
+
+                        if (! detail::overlay::select_next_ip(
+                                        operation,
+                                        *current,
+                                        current_seg_id,
+                                        boost::end(it->operations),
+                                        current_iit))
+                        {
+                            detail::overlay::backtrack(
+                                size_at_start, fail,
+                                rings, current_output, turns, *it,
+                                "Dead end at start",
+                                geometry1, geometry2);
+                        }
+                        else
+                        {
+
+                            iit->visited.set_started();
+
+                            unsigned int i = 0;
+
+                            while (current_iit != iit && ! fail)
+                            {
+                                if (current_iit->visited.visited())
+                                {
+                                    // It visits a visited node again, without passing the start node.
+                                    // This makes it suspicious for endless loops
+                                    detail::overlay::backtrack(
+                                        size_at_start, fail,
+                                        rings,  current_output, turns, *it,
+                                        "Visit again",
+                                        geometry1, geometry2);
+                                }
+                                else
+                                {
 
 
-                    // We assume clockwise polygons only, non self-intersecting, closed.
-                    // However, the input might be different, and checking validity
-                    // is up to the library user.
+                                    // We assume clockwise polygons only, non self-intersecting, closed.
+                                    // However, the input might be different, and checking validity
+                                    // is up to the library user.
 
-                    // Therefore we make here some sanity checks. If the input
-                    // violates the assumptions, the output polygon will not be correct
-                    // but the routine will stop and output the current polygon, and
-                    // will continue with the next one.
+                                    // Therefore we make here some sanity checks. If the input
+                                    // violates the assumptions, the output polygon will not be correct
+                                    // but the routine will stop and output the current polygon, and
+                                    // will continue with the next one.
 
-                    // Below three reasons to stop.
-                    assign_next_ip(geometry1, geometry2,
-                        turn_points, current, current_output,
-                        *current_iit, current_seg_id);
+                                    // Below three reasons to stop.
+                                    assign_next_ip(geometry1, geometry2,
+                                        turns, current, current_output,
+                                        *current_iit, current_seg_id);
 
-                    if (! detail::overlay::select_next_ip(
-                                operation,
-                                *current,
-                                current_seg_id,
-                                iit,
-                                current_iit))
-                    {
-                        // Should not occur in valid (non-self-intersecting) polygons
-                        // Should not occur in self-intersecting polygons without spikes
-                        // Might occur in polygons with spikes
-                        detail::overlay::stop_gracefully(
-                            current_output, stop, "Dead end");
-                    }
+                                    if (! detail::overlay::select_next_ip(
+                                                operation,
+                                                *current,
+                                                current_seg_id,
+                                                iit,
+                                                current_iit))
+                                    {
+                                        // Should not occur in valid (non-self-intersecting) polygons
+                                        // Should not occur in self-intersecting polygons without spikes
+                                        // Might occur in polygons with spikes
+                                        detail::overlay::backtrack(
+                                            size_at_start, fail,
+                                            rings,  current_output, turns, *it,
+                                            "Dead end",
+                                            geometry1, geometry2);
+                                    }
 
-                    if (i++ > turn_points.size())
-                    {
-                        // Sanity check: there may be never more loops
-                        // than intersection points.
-                        detail::overlay::stop_gracefully(
-                            current_output, stop, "Endless loop");
+                                    if (i++ > turns.size())
+                                    {
+                                        // Sanity check: there may be never more loops
+                                        // than intersection points.
+                                        detail::overlay::backtrack(
+                                            size_at_start, fail,
+                                            rings,  current_output, turns, *it,
+                                            "Endless loop",
+                                            geometry1, geometry2);
+                                    }
+                                }
+                            }
+
+                            if (! fail)
+                            {
+                                iit->visited.set_finished();
+                                rings.push_back(current_output);
+                            }
+                        }
                     }
                 }
-
-                iit->visited.set_finished();
-
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-                //std::cout << "traversal: finish: " << *current;
-                std::cout << geometry::wkt(current_output) << std::endl;
-#endif
-
-                *out++ = current_output;
             }
         }
-    }
-#ifdef BOOST_GEOMETRY_DEBUG_INTERSECTION
-    std::cout << "traversal: all IP's handled" << std::endl;
-#endif
+    } while (fail);
 }
 
 
diff --git a/include/boost/geometry/algorithms/union.hpp b/include/boost/geometry/algorithms/union.hpp
index 9fca810a9..5f9f88c44 100644
--- a/include/boost/geometry/algorithms/union.hpp
+++ b/include/boost/geometry/algorithms/union.hpp
@@ -60,7 +60,7 @@ template
     typename Strategy
 >
 struct union_inserter
-    : detail::overlay::overlay_and_assemble
+    : detail::overlay::overlay
         <G1, G2, OutputIterator, GeometryOut, 1, Strategy>
 {
 };
diff --git a/include/boost/geometry/multi/algorithms/overlay/assemble.hpp b/include/boost/geometry/multi/algorithms/overlay/assemble.hpp
index ca014b8d2..10c712292 100644
--- a/include/boost/geometry/multi/algorithms/overlay/assemble.hpp
+++ b/include/boost/geometry/multi/algorithms/overlay/assemble.hpp
@@ -43,7 +43,7 @@ struct add_to_containment<multi_polygon_tag, MultiPolygon>
     template <typename ContainmentContainer, typename Map>
     static inline void apply(ContainmentContainer& container,
             ring_identifier const& id, MultiPolygon const& multi_polygon,
-            Map const& map)
+            Map const& map, bool dissolve)
     {
         ring_identifier copy = id;
         copy.multi_index = 0;
@@ -58,7 +58,7 @@ struct add_to_containment<multi_polygon_tag, MultiPolygon>
             <
                 polygon_tag,
                 typename boost::range_value<MultiPolygon>::type
-            >::apply(container, copy, *it, map);
+            >::apply(container, copy, *it, map, dissolve);
         }
     }
 };