[overlay] Fix possible instabilities in sort_by_side

plus some minor code refactoring

Author: barendgehrels

include/boost/geometry/algorithms/detail/overlay/handle_colocations.hpp

@@ -752,6 +752,46 @@ struct is_turn_index
     signed_size_type m_index;
 };
 
+template
+<
+    typename Sbs,
+    typename Point,
+    typename Turns,
+    typename Geometry1,
+    typename Geometry2
+>
+inline bool fill_sbs(Sbs& sbs, Point& turn_point,
+                     cluster_info const& cinfo,
+                     Turns const& turns,
+                     Geometry1 const& geometry1, Geometry2 const& geometry2)
+{
+    typedef typename boost::range_value<Turns>::type turn_type;
+
+    std::set<signed_size_type> const& ids = cinfo.turn_indices;
+
+    if (ids.empty())
+    {
+        return false;
+    }
+
+    bool first = true;
+    for (std::set<signed_size_type>::const_iterator sit = ids.begin();
+         sit != ids.end(); ++sit)
+    {
+        signed_size_type turn_index = *sit;
+        turn_type const& turn = turns[turn_index];
+        if (first )
+        {
+            turn_point = turn.point;
+        }
+        for (int i = 0; i < 2; i++)
+        {
+            sbs.add(turn.operations[i], turn_index, i, geometry1, geometry2, first);
+            first = false;
+        }
+    }
+    return true;
+}
 
 template
 <
@@ -783,32 +823,14 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns,
          mit != clusters.end(); ++mit)
     {
         cluster_info& cinfo = mit->second;
-        std::set<signed_size_type> const& ids = cinfo.turn_indices;
-        if (ids.empty())
-        {
-            continue;
-        }
 
         sbs_type sbs(strategy);
         point_type turn_point; // should be all the same for all turns in cluster
-
-        bool first = true;
-        for (std::set<signed_size_type>::const_iterator sit = ids.begin();
-             sit != ids.end(); ++sit)
+        if (! fill_sbs(sbs, turn_point, cinfo, turns, geometry1, geometry2))
         {
-            signed_size_type turn_index = *sit;
-            turn_type const& turn = turns[turn_index];
-            if (first)
-            {
-                turn_point = turn.point;
-            }
-            for (int i = 0; i < 2; i++)
-            {
-                turn_operation_type const& op = turn.operations[i];
-                sbs.add(op, turn_index, i, geometry1, geometry2, first);
-                first = false;
-            }
+            continue;
         }
+
         sbs.apply(turn_point);
 
         sbs.find_open();
@@ -823,7 +845,7 @@ inline void gather_cluster_properties(Clusters& clusters, Turns& turns,
         // polygons
         for (std::size_t i = 0; i < sbs.m_ranked_points.size(); i++)
         {
-            const typename sbs_type::rp& ranked = sbs.m_ranked_points[i];
+            typename sbs_type::rp const& ranked = sbs.m_ranked_points[i];
             turn_type& turn = turns[ranked.turn_index];
             turn_operation_type& op = turn.operations[ranked.operation_index];
 

include/boost/geometry/algorithms/detail/overlay/sort_by_side.hpp

@@ -186,7 +186,14 @@ struct less_by_side
             return on_same(first, second);
         }
 
-        int const side_first_wrt_second = -side_second_wrt_first;
+        int const side_first_wrt_second = m_strategy.apply(m_turn_point, second.point, first.point);
+        if (side_second_wrt_first != -side_first_wrt_second)
+        {
+            // (FP) accuracy error in side calculation, the sides are not opposite.
+            // In that case they can be handled as collinear.
+            // If not, then the sort-order might not be stable.
+            return on_same(first, second);
+        }
 
         // Both are on same side, and not collinear
         // Union: return true if second is right w.r.t. first, so -1,
@@ -367,9 +374,11 @@ public :
         }
     }
 
-    template <signed_size_type segment_identifier::*Member, typename Map>
-    void find_open_generic(Map& handled, bool check)
+    void find_open_by_piece_index()
     {
+        // For buffers, use piece index
+        std::set<signed_size_type> handled;
+
         for (std::size_t i = 0; i < m_ranked_points.size(); i++)
         {
             const rp& ranked = m_ranked_points[i];
@@ -378,39 +387,47 @@ public :
                 continue;
             }
 
-            signed_size_type const& index = ranked.seg_id.*Member;
-            if (check && (index < 0 || index > 1))
+            signed_size_type const& index = ranked.seg_id.piece_index;
+            if (handled.count(index) > 0)
             {
-                // Should not occur
                 continue;
             }
-            if (! handled[index])
+            find_polygons_for_source<&segment_identifier::piece_index>(index, i);
+            handled.insert(index);
+        }
+    }
+
+    void find_open_by_source_index()
+    {
+        // Check for source index 0 and 1
+        bool handled[2] = {false, false};
+        for (std::size_t i = 0; i < m_ranked_points.size(); i++)
+        {
+            const rp& ranked = m_ranked_points[i];
+            if (ranked.direction != dir_from)
+            {
+                continue;
+            }
+
+            signed_size_type const& index = ranked.seg_id.source_index;
+            if (index < 0 || index > 1 || handled[index])
             {
-                find_polygons_for_source<Member>(index, i);
-                handled[index] = true;
+                continue;
             }
+            find_polygons_for_source<&segment_identifier::source_index>(index, i);
+            handled[index] = true;
         }
     }
 
     void find_open()
     {
         if (BOOST_GEOMETRY_CONDITION(OverlayType == overlay_buffer))
         {
-            // For buffers, use piece index
-            std::map<signed_size_type, bool> handled;
-            find_open_generic
-                <
-                    &segment_identifier::piece_index
-                >(handled, false);
+            find_open_by_piece_index();
         }
         else
         {
-            // For other operations, by source (there should only source 0,1)
-            bool handled[2] = {false, false};
-            find_open_generic
-                <
-                    &segment_identifier::source_index
-                >(handled, true);
+            find_open_by_source_index();
         }
     }