Merge pull request #736 from barendgehrels/fix/imperfect-touch

barendgehrels · web-flow · commit 9759656d2bb7 · 2020-07-22T11:57:35.000+02:00
Fix errors in some floating point cases
diff --git a/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp b/include/boost/geometry/algorithms/detail/overlay/get_turn_info.hpp
@@ -111,6 +111,60 @@ struct base_turn_handler
         both(ti, condition ? operation_union : operation_intersection);
     }
 
+
+#if ! defined(BOOST_GEOMETRY_USE_RESCALING)
+    template
+    <
+        typename UniqueSubRange1,
+        typename UniqueSubRange2
+    >
+    static inline int side_with_distance_measure(UniqueSubRange1 const& range_p,
+            UniqueSubRange2 const& range_q,
+            int range_index, int point_index)
+    {
+        if (range_index >= 1 && range_p.is_last_segment())
+        {
+            return 0;
+        }
+        if (point_index >= 2 && range_q.is_last_segment())
+        {
+            return 0;
+        }
+
+        typedef typename select_coordinate_type
+            <
+                typename UniqueSubRange1::point_type,
+                typename UniqueSubRange2::point_type
+            >::type coordinate_type;
+
+        typedef detail::distance_measure<coordinate_type> dm_type;
+
+        dm_type const dm = get_distance_measure(range_p.at(range_index), range_p.at(range_index + 1), range_q.at(point_index));
+        return dm.measure == 0 ? 0 : dm.measure > 0 ? 1 : -1;
+    }
+
+    template
+    <
+        typename UniqueSubRange1,
+        typename UniqueSubRange2
+    >
+    static inline int verified_side(int side,
+                                    UniqueSubRange1 const& range_p,
+                                    UniqueSubRange2 const& range_q,
+                                    int range_index,
+                                    int point_index)
+    {
+        return side == 0 ? side_with_distance_measure(range_p, range_q, range_index, point_index) : side;
+    }
+#else
+    template <typename T1, typename T2>
+    static inline int verified_side(int side, T1 const& , T2 const& , int , int)
+    {
+        return side;
+    }
+#endif
+
+
     template <typename TurnInfo, typename IntersectionInfo>
     static inline void assign_point(TurnInfo& ti,
                 method_type method,
@@ -178,7 +232,7 @@ struct base_turn_handler
                 <
                     typename UniqueSubRange1::point_type,
                     typename UniqueSubRange2::point_type
-                    >::type
+                >::type
             > dm_type;
 
         const bool p_closer =
@@ -369,14 +423,80 @@ struct touch : public base_turn_handler
         return side1 == side2 && ! opposite(side1, turn);
     }
 
-    /*static inline void block_second(bool block, TurnInfo& ti)
+#if ! defined(BOOST_GEOMETRY_USE_RESCALING)
+    template
+    <
+        typename UniqueSubRange1,
+        typename UniqueSubRange2
+    >
+    static inline bool handle_imperfect_touch(UniqueSubRange1 const& range_p,
+            UniqueSubRange2 const& range_q, TurnInfo& ti)
     {
-        if (block)
+        //  Q
+        //  ^
+        // ||
+        // ||
+        // |^----
+        // >----->P
+        // *            * they touch here (P/Q are (nearly) on top)
+        //
+        // Q continues from where P comes.
+        // P continues from where Q comes
+        // This is often a blocking situation,
+        // unless there are FP issues: there might be a distance
+        // between Pj and Qj, in that case handle it as a union.
+        //
+        // Exaggerated:
+        //  Q
+        //  ^           Q is nearly vertical
+        //   \          but not completely - and still ends above P
+        // |  \qj       In this case it should block P and
+        // |  ^------   set Q to Union
+        // >----->P     qj is LEFT of P1 and pi is LEFT of Q2
+        //              (the other way round is also possible)
+
+        typedef typename select_coordinate_type
+            <
+                typename UniqueSubRange1::point_type,
+                typename UniqueSubRange2::point_type
+            >::type coordinate_type;
+
+        typedef detail::distance_measure<coordinate_type> dm_type;
+
+        dm_type const dm_qj_p1 = get_distance_measure(range_p.at(0), range_p.at(1), range_q.at(1));
+        dm_type const dm_pi_q2 = get_distance_measure(range_q.at(1), range_q.at(2), range_p.at(0));
+
+        if (dm_qj_p1.measure > 0 && dm_pi_q2.measure > 0)
         {
-            ti.operations[1].operation = operation_blocked;
+            // Even though there is a touch, Q(j) is left of P1
+            // and P(i) is still left from Q2.
+            // It can continue.
+            ti.operations[0].operation = operation_blocked;
+            // Q turns right -> union (both independent),
+            // Q turns left -> intersection
+            ti.operations[1].operation = operation_union;
+            ti.touch_only = true;
+            return true;
         }
-    }*/
 
+        dm_type const dm_pj_q1 = get_distance_measure(range_q.at(0), range_q.at(1), range_p.at(1));
+        dm_type const dm_qi_p2 = get_distance_measure(range_p.at(1), range_p.at(2), range_q.at(0));
+
+        if (dm_pj_q1.measure > 0 && dm_qi_p2.measure > 0)
+        {
+            // Even though there is a touch, Q(j) is left of P1
+            // and P(i) is still left from Q2.
+            // It can continue.
+            ti.operations[0].operation = operation_union;
+            // Q turns right -> union (both independent),
+            // Q turns left -> intersection
+            ti.operations[1].operation = operation_blocked;
+            ti.touch_only = true;
+            return true;
+        }
+        return false;
+    }
+#endif
 
     template
     <
@@ -400,9 +520,10 @@ struct touch : public base_turn_handler
         bool const has_pk = ! range_p.is_last_segment();
         bool const has_qk = ! range_q.is_last_segment();
 
-        int const side_qi_p1 = dir_info.sides.template get<1, 0>();
-        int const side_qk_p1 = has_qk ? side.qk_wrt_p1() : 0;
+        int const side_pk_q1 = has_pk ? side.pk_wrt_q1() : 0;
 
+        int const side_qi_p1 = verified_side(dir_info.sides.template get<1, 0>(), range_p, range_q, 0, 0);
+        int const side_qk_p1 = has_qk ? verified_side(side.qk_wrt_p1(), range_p, range_q, 0, 2) : 0;
 
         // If Qi and Qk are both at same side of Pi-Pj,
         // or collinear (so: not opposite sides)
@@ -413,6 +534,7 @@ struct touch : public base_turn_handler
             int const side_qk_q  = has_qk ? side.qk_wrt_q1() : 0;
 
             bool const q_turns_left = side_qk_q == 1;
+
             bool const block_q = side_qk_p1 == 0
                         && ! same(side_qi_p1, side_qk_q)
                         ;
@@ -422,9 +544,18 @@ struct touch : public base_turn_handler
             // or Q is fully collinear && P turns not to left
             if (side_pk_p == side_qi_p1
                 || side_pk_p == side_qk_p1
-                || (side_qi_p1 == 0 && side_qk_p1 == 0 && side_pk_p != -1)
-                )
+                || (side_qi_p1 == 0 && side_qk_p1 == 0 && side_pk_p != -1))
             {
+#if ! defined(BOOST_GEOMETRY_USE_RESCALING)
+                if (side_qk_p1 == 0 && side_pk_q1 == 0
+                    && has_qk && has_qk
+                    && handle_imperfect_touch(range_p, range_q, ti))
+                {
+                    // If q continues collinearly (opposite) with p, it should be blocked
+                    // but (FP) not if there is just a tiny space in between
+                    return;
+                }
+#endif
                 // Collinear -> lines join, continue
                 // (#BRL2)
                 if (side_pk_q2 == 0 && ! block_q)
@@ -433,8 +564,6 @@ struct touch : public base_turn_handler
                     return;
                 }
 
-                int const side_pk_q1 = has_pk && has_qk ? side.pk_wrt_q1() : 0;
-
                 // Collinear opposite case -> block P
                 // (#BRL4, #BLR8)
                 if (side_pk_q1 == 0)
@@ -457,7 +586,6 @@ struct touch : public base_turn_handler
                     {
                         ti.operations[1].operation = operation_blocked;
                     }
-                    //block_second(block_q, ti);
                     return;
                 }
 
@@ -483,7 +611,6 @@ struct touch : public base_turn_handler
                     {
                         ti.touch_only = true;
                     }
-                    //block_second(block_q, ti);
                     return;
                 }
             }
@@ -508,8 +635,9 @@ struct touch : public base_turn_handler
         }
         else
         {
+            // The qi/qk are opposite to each other, w.r.t. p1
             // From left to right or from right to left
-            int const side_pk_p = has_pk ? side.pk_wrt_p1() : 0;
+            int const side_pk_p = has_pk ? verified_side(side.pk_wrt_p1(), range_p, range_p, 0, 2) : 0;
             bool const right_to_left = side_qk_p1 == 1;
 
             // If p turns into direction of qi (1,2)
@@ -601,23 +729,23 @@ struct equal : public base_turn_handler
             // Without rescaling, to check for union/intersection,
             // try to check side values (without any thresholds)
             typedef typename select_coordinate_type
-                    <
-                        typename UniqueSubRange1::point_type,
-                        typename UniqueSubRange2::point_type
-                        >::type coordinate_type;
+                <
+                    typename UniqueSubRange1::point_type,
+                    typename UniqueSubRange2::point_type
+                >::type coordinate_type;
 
             typedef detail::distance_measure<coordinate_type> dm_type;
 
-            dm_type const dm_qk_p
+            dm_type const dm_pk_q2
                = get_distance_measure(range_q.at(1), range_q.at(2), range_p.at(2));
-            dm_type const dm_pk_q
+            dm_type const dm_qk_p2
                = get_distance_measure(range_p.at(1), range_p.at(2), range_q.at(2));
 
-            if (dm_pk_q.measure != dm_qk_p.measure)
+            if (dm_qk_p2.measure != dm_pk_q2.measure)
             {
                 // A (possibly very small) difference is detected, which
                 // can be used to distinguish between union/intersection
-                ui_else_iu(dm_pk_q.measure < dm_qk_p.measure, ti);
+                ui_else_iu(dm_qk_p2.measure < dm_pk_q2.measure, ti);
                 return;
             }
         }
diff --git a/include/boost/geometry/strategies/cartesian/intersection.hpp b/include/boost/geometry/strategies/cartesian/intersection.hpp
@@ -371,6 +371,38 @@ struct cartesian_segments
         return unified<ratio_type>(sinfo, p, q, policy, modelled_range_p, modelled_range_q);
     }
 
+    //! Returns true if two segments do not overlap.
+    //! If not, then no further calculations need to be done.
+    template
+    <
+        std::size_t Dimension,
+        typename CoordinateType,
+        typename PointP,
+        typename PointQ
+    >
+    static inline bool disjoint_by_range(PointP const& p1, PointP const& p2,
+                                         PointQ const& q1, PointQ const& q2)
+    {
+        CoordinateType minp = get<Dimension>(p1);
+        CoordinateType maxp = get<Dimension>(p2);
+        CoordinateType minq = get<Dimension>(q1);
+        CoordinateType maxq = get<Dimension>(q2);
+        if (minp > maxp)
+        {
+            std::swap(minp, maxp);
+        }
+        if (minq > maxq)
+        {
+            std::swap(minq, maxq);
+        }
+
+        // In this case, max(p) < min(q)
+        //     P         Q
+        // <-------> <------->
+        // (and the space in between is not extremely small)
+        return math::smaller(maxp, minq) || math::smaller(maxq, minp);
+    }
+
     // Implementation for either rescaled or non rescaled versions.
     template
     <
@@ -413,6 +445,12 @@ struct cartesian_segments
                 ;
         }
 
+        if (disjoint_by_range<0, coordinate_type>(p1, p2, q1, q2)
+         || disjoint_by_range<1, coordinate_type>(p1, p2, q1, q2))
+        {
+            return Policy::disjoint();
+        }
+
         side_info sides;
         sides.set<0>(side_strategy_type::apply(q1, q2, p1),
                      side_strategy_type::apply(q1, q2, p2));
@@ -467,8 +505,8 @@ struct cartesian_segments
             {
                 // If this is the case, no rescaling is done for FP precision.
                 // We set it to collinear, but it indicates a robustness issue.
-                sides.set<0>(0,0);
-                sides.set<1>(0,0);
+                sides.set<0>(0, 0);
+                sides.set<1>(0, 0);
                 collinear = true;
             }
             else
diff --git a/test/algorithms/buffer/buffer_linestring_aimes.cpp b/test/algorithms/buffer/buffer_linestring_aimes.cpp
@@ -483,8 +483,8 @@ int test_main(int, char* [])
     test_aimes<bg::model::point<default_test_type, 2, bg::cs::cartesian> >();
 
 #if defined(BOOST_GEOMETRY_TEST_FAILURES)
-    // Non-rescaled reports 22 failures (failures in validity only)
-    BoostGeometryWriteExpectedFailures(BG_NO_FAILURES, 22);
+    // Non-rescaled reports failures, but in validity only
+    BoostGeometryWriteExpectedFailures(BG_NO_FAILURES, 4);
 #endif
 
     return 0;
diff --git a/test/algorithms/buffer/buffer_multi_polygon.cpp b/test/algorithms/buffer/buffer_multi_polygon.cpp
@@ -307,13 +307,6 @@ static std::string const mysql_report_2015_07_05_1
 static std::string const mysql_report_2015_07_05_2
     = "MULTIPOLYGON(((19777 -21893,3.22595e+307 6.86823e+307,-40 -13,19777 -21893)),((-1322 4851,8.49998e+307 3.94481e+307,75 -69,8.64636e+307 3.94909e+307,-1.15292e+18 7.20576e+16,-1322 4851)))";
 
-#if ! defined(BOOST_GEOMETRY_USE_RESCALING) \
-    && defined(BOOST_GEOMETRY_USE_KRAMER_RULE) \
-    && ! defined(BOOST_GEOMETRY_TEST_FAILURES)
-// These testcases are failing for non-rescaled Kramer rule
-#define BOOST_GEOMETRY_EXCLUDE
-#endif
-
 template <bool Clockwise, typename P>
 void test_all()
 {
@@ -453,18 +446,14 @@ void test_all()
     test_one<multi_polygon_type, polygon_type>("rt_p19", rt_p19, join_miter, end_flat, 25.5637, 1.0);
     test_one<multi_polygon_type, polygon_type>("rt_p20", rt_p20, join_miter, end_flat, 25.4853, 1.0);
     test_one<multi_polygon_type, polygon_type>("rt_p21", rt_p21, join_miter, end_flat, 17.1716, 1.0);
-#if ! defined(BOOST_GEOMETRY_EXCLUDE)
     test_one<multi_polygon_type, polygon_type>("rt_p22", rt_p22, join_miter, end_flat, 26.5711, 1.0);
-#endif
 
     test_one<multi_polygon_type, polygon_type>("rt_q1", rt_q1, join_miter, end_flat, 27, 1.0);
     test_one<multi_polygon_type, polygon_type>("rt_q2", rt_q2, join_miter, end_flat, 26.4853, 1.0);
     test_one<multi_polygon_type, polygon_type>("rt_q2", rt_q2, join_miter, end_flat, 0.9697, -0.25);
 
     test_one<multi_polygon_type, polygon_type>("rt_r", rt_r, join_miter, end_flat, 21.0761, 1.0);
-#if ! defined(BOOST_GEOMETRY_EXCLUDE)
     test_one<multi_polygon_type, polygon_type>("rt_s1", rt_s1, join_miter, end_flat, 20.4853, 1.0);
-#endif
 
     test_one<multi_polygon_type, polygon_type>("rt_s2", rt_s2, join_miter, end_flat, 24.6495, 1.0);
 
@@ -543,7 +532,7 @@ int test_main(int, char* [])
 #endif
     
 #if defined(BOOST_GEOMETRY_TEST_FAILURES)
-    BoostGeometryWriteExpectedFailures(1, 5);
+    BoostGeometryWriteExpectedFailures(1, 1);
 #endif
 
     return 0;
diff --git a/test/algorithms/set_operations/difference/difference.cpp b/test/algorithms/set_operations/difference/difference.cpp
diff --git a/test/algorithms/set_operations/union/union.cpp b/test/algorithms/set_operations/union/union.cpp
