[get_turn_info] handle collinear as equal if both segments arrive at intersection point

Author: barendgehrels

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

@@ -186,7 +186,8 @@ struct base_turn_handler
     {
         ti.method = method;
 
-        // For touch/touch interior always take the intersection point 0 (there is only one).
+        // For touch/touch interior always take the intersection point 0
+        // (usually there is only one - but if collinear is handled as touch, both could be taken).
         static int const index = 0;
 
         geometry::convert(info.intersections[index], ti.point);
@@ -965,8 +966,57 @@ template
 >
 struct collinear : public base_turn_handler
 {
+    template
+    <
+        typename IntersectionInfo,
+        typename UniqueSubRange1,
+        typename UniqueSubRange2,
+        typename DirInfo
+    >
+    static bool handle_as_equal(IntersectionInfo const& info,
+                                UniqueSubRange1 const& range_p,
+                                UniqueSubRange2 const& range_q,
+                                DirInfo const& dir_info)
+    {
+#if defined(BOOST_GEOMETRY_USE_RESCALING)
+        return false;
+#endif
+        int const arrival_p = dir_info.arrival[0];
+        int const arrival_q = dir_info.arrival[1];
+        if (arrival_p * arrival_q != -1 || info.count != 2)
+        {
+            // Code below assumes that either p or q arrives in the other segment
+            return false;
+        }
+
+        auto const location = distance_measure(info.intersections[1],
+                arrival_p == 1 ? range_q.at(1) : range_p.at(1));
+        decltype(location) const zero = 0;
+        return math::equals(location, zero);
+    }
+
     /*
-        arrival P   pk//p1  qk//q1   product*  case    result
+        Either P arrives within Q (arrival_p == -1) or Q arrives within P.
+
+        Typical situation:
+              ^q2
+             /
+            ^q1
+           /         ____ ip[1]
+          //|p1  } this section of p/q is colllinear
+       q0// |    }   ____ ip[0]
+         /  |
+        /   v
+       p0   p2
+
+       P arrives (at p1) in segment Q (between q0 and q1).
+       Therefore arrival_p == 1
+       P (p2) goes to the right (-1). Follow P for intersection, or follow Q for union.
+       Therefore if (arrival_p==1) and side_p==-1, result = iu
+
+       Complete table:
+
+        arrival P   pk//p1  qk//q1   product   case    result
          1           1                1        CLL1    ui
         -1                   1       -1        CLL2    iu
          1           1                1        CLR1    ui
@@ -980,7 +1030,9 @@ struct collinear : public base_turn_handler
          1           0                0        CC1     cc
         -1                   0        0        CC2     cc
 
-         *product = arrival * (pk//p1 or qk//q1)
+         Resulting in the rule:
+         The arrival-info multiplied by the relevant side delivers the result.
+         product = arrival * (pk//p1 or qk//q1)
 
          Stated otherwise:
          - if P arrives: look at turn P
@@ -989,13 +1041,6 @@ struct collinear : public base_turn_handler
          - if P arrives and P turns right: intersection for P
          - if Q arrives and Q turns left: union for Q (=intersection for P)
          - if Q arrives and Q turns right: intersection for Q (=union for P)
-
-         ROBUSTNESS: p and q are collinear, so you would expect
-         that side qk//p1 == pk//q1. But that is not always the case
-         in near-epsilon ranges. Then decision logic is different.
-         If p arrives, q is further, so the angle qk//p1 is (normally)
-         more precise than pk//p1
-
     */
     template
     <
@@ -1016,29 +1061,25 @@ struct collinear : public base_turn_handler
         // Copy the intersection point in TO direction
         assign_point(ti, method_collinear, info, non_opposite_to_index(info));
 
-        int const arrival = dir_info.arrival[0];
+        int const arrival_p = dir_info.arrival[0];
         // Should not be 0, this is checked before
-        BOOST_GEOMETRY_ASSERT(arrival != 0);
+        BOOST_GEOMETRY_ASSERT(arrival_p != 0);
 
         bool const has_pk = ! range_p.is_last_segment();
         bool const has_qk = ! range_q.is_last_segment();
         int const side_p = has_pk ? side.pk_wrt_p1() : 0;
         int const side_q = has_qk ? side.qk_wrt_q1() : 0;
 
         // If p arrives, use p, else use q
-        int const side_p_or_q = arrival == 1
+        int const side_p_or_q = arrival_p == 1
             ? side_p
             : side_q
             ;
 
-        // See comments above,
-        // resulting in a strange sort of mathematic rule here:
-        // The arrival-info multiplied by the relevant side
-        // delivers a consistent result.
+        // Calculate product according to comments above.
+        int const product = arrival_p * side_p_or_q;
 
-        int const product = arrival * side_p_or_q;
-
-        if(product == 0)
+        if (product == 0)
         {
             both(ti, operation_continue);
         }
@@ -1186,11 +1227,11 @@ private :
     {
         TurnInfo tp = tp_model;
 
-        int const p_arrival = info.d_info().arrival[0];
-        int const q_arrival = info.d_info().arrival[1];
+        int const arrival_p = info.d_info().arrival[0];
+        int const arrival_q = info.d_info().arrival[1];
 
         // If P arrives within Q, there is a turn dependent on P
-        if ( p_arrival == 1
+        if ( arrival_p == 1
           && ! range_p.is_last_segment()
           && set_tp<0>(side.pk_wrt_p1(), tp, info.i_info()) )
         {
@@ -1200,7 +1241,7 @@ private :
         }
 
         // If Q arrives within P, there is a turn dependent on Q
-        if ( q_arrival == 1
+        if ( arrival_q == 1
           && ! range_q.is_last_segment()
           && set_tp<1>(side.qk_wrt_q1(), tp, info.i_info()) )
         {
@@ -1212,8 +1253,8 @@ private :
         if (AssignPolicy::include_opposite)
         {
             // Handle cases not yet handled above
-            if ((q_arrival == -1 && p_arrival == 0)
-                || (p_arrival == -1 && q_arrival == 0))
+            if ((arrival_q == -1 && arrival_p == 0)
+                || (arrival_p == -1 && arrival_q == 0))
             {
                 for (unsigned int i = 0; i < 2; i++)
                 {
@@ -1420,9 +1461,10 @@ struct get_turn_info
             // Collinear
             if ( ! inters.d_info().opposite )
             {
-
-                if ( inters.d_info().arrival[0] == 0 )
+                if (inters.d_info().arrival[0] == 0
+                    || collinear<TurnInfo>::handle_as_equal(inters.i_info(), range_p, range_q, inters.d_info()))
                 {
+                    // Both segments arrive at the second intersection point
                     handle_as_equal = true;
                 }
                 else