Merge pull request #331 from mkaravel/fix/is_valid_numerically_collinear

Fix problem of numerical inconsistency when checking for spikes

include/boost/geometry/algorithms/detail/point_is_spike_or_equal.hpp

@@ -1,9 +1,14 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// 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) 2007-2015 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2008-2015 Bruno Lalande, Paris, France.
+// Copyright (c) 2009-2015 Mateusz Loskot, London, UK.
+// Copyright (c) 2013-2015 Adam Wulkiewicz, Lodz, Poland.
+
+// This file was modified by Oracle on 2015.
+// Modifications copyright (c) 2015 Oracle and/or its affiliates.
+
+// Contributed and/or modified by Menelaos Karavelas, 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
@@ -12,8 +17,6 @@
 #ifndef BOOST_GEOMETRY_ALGORITHMS_DETAIL_POINT_IS_EQUAL_OR_SPIKE_HPP
 #define BOOST_GEOMETRY_ALGORITHMS_DETAIL_POINT_IS_EQUAL_OR_SPIKE_HPP
 
-#include <boost/geometry/arithmetic/arithmetic.hpp>
-#include <boost/geometry/algorithms/detail/convert_point_to_point.hpp>
 #include <boost/geometry/algorithms/detail/recalculate.hpp>
 #include <boost/geometry/policies/robustness/robust_point_type.hpp>
 #include <boost/geometry/strategies/side.hpp>
@@ -28,6 +31,17 @@ namespace boost { namespace geometry
 namespace detail
 {
 
+template <std::size_t Index, typename Point1, typename Point2>
+inline int sign_of_difference(Point1 const& point1, Point2 const& point2)
+{
+    return
+        math::equals(geometry::get<Index>(point1), geometry::get<Index>(point2))
+        ?
+        0
+        :
+        (geometry::get<Index>(point1) > geometry::get<Index>(point2) ? 1 : -1);
+}
+
 // Checks if a point ("last_point") causes a spike w.r.t.
 // the specified two other points (segment_a, segment_b)
 //
@@ -35,7 +49,9 @@ namespace detail
 //  a       lp     b
 //
 // Above, lp generates a spike w.r.t. segment(a,b)
-// So specify last point first, then (a,b) (this is unordered, so unintuitive)
+// So specify last point first, then (a,b)
+// The segment's orientation does matter: if lp is to the right of b
+// no spike is reported
 template <typename Point1, typename Point2, typename Point3>
 static inline bool point_is_spike_or_equal(Point1 const& last_point,
             Point2 const& segment_a,
@@ -46,29 +62,21 @@ static inline bool point_is_spike_or_equal(Point1 const& last_point,
         typename cs_tag<Point1>::type
     >::type side_strategy;
 
-    typedef Point1 vector_type;
-
     int const side = side_strategy::apply(last_point, segment_a, segment_b);
     if (side == 0)
     {
         // Last point is collinear w.r.t previous segment.
         // Check if it is equal
-        vector_type diff1;
-        conversion::convert_point_to_point(last_point, diff1);
-        geometry::subtract_point(diff1, segment_b);
-        int const sgn_x1 = math::sign(geometry::get<0>(diff1));
-        int const sgn_y1 = math::sign(geometry::get<1>(diff1));
+        int const sgn_x1 = sign_of_difference<0>(last_point, segment_b);
+        int const sgn_y1 = sign_of_difference<1>(last_point, segment_b);
         if (sgn_x1 == 0 && sgn_y1 == 0)
         {
             return true;
         }
 
         // Check if it moves forward
-        vector_type diff2;
-        conversion::convert_point_to_point(segment_b, diff2);
-        geometry::subtract_point(diff2, segment_a);
-        int const sgn_x2 = math::sign(geometry::get<0>(diff2));
-        int const sgn_y2 = math::sign(geometry::get<1>(diff2));
+        int const sgn_x2 = sign_of_difference<0>(segment_b, segment_a);
+        int const sgn_y2 = sign_of_difference<1>(segment_b, segment_a);
 
         return sgn_x1 != sgn_x2 || sgn_y1 != sgn_y2;
     }

test/algorithms/is_valid.cpp

@@ -713,6 +713,67 @@ inline void test_open_polygons()
                    0.7808688094430304 -0.6246950475544243,\
                    0.7808688094430304 -0.6246950475544243))",
          AllowDuplicates);
+
+
+    // MySQL report on Sep 30, 2015
+    test::apply
+        ("pg077",
+         "POLYGON((72.8714768817168 -167.0048853643874,9274.40641550926 3433.5957427942167,-58.09039811390054 187.50989457746405,-81.09039811390053 179.50989457746405,-207.99999999999997 135.36742435621204,-208 1,-208 0,-208 -276.9111154485375,49.8714768817168 -176.0048853643874))",
+         true);
+
+    test::apply("pg077-simplified",
+                "POLYGON((-200 0,-207.99999999999997 135.36742435621204,-208 1,-208 0,-208 -276.9111154485375))",
+                true);
+
+    test::apply
+        ("pg078",
+         "POLYGON((0 10,-10 0,0 0,10 0))",
+         true);
+
+    test::apply
+        ("pg078spike1",
+         "POLYGON((0 10,-10 0,0 0,-10 0,10 0))",
+         false);
+
+    test::apply
+        ("pg078spike2",
+         "POLYGON((0 10,-10 0,0 0,-8 0,10 0))",
+         false);
+
+    test::apply
+        ("pg078spike3",
+         "POLYGON((0 10,-10 0,0 0,-11 0,10 0))",
+         false);
+
+    test::apply
+        ("pg078reversed",
+         "POLYGON((0 10,10 0,0 0,-10 0))",
+         false);
+
+    test::apply
+        ("pg079",
+         "POLYGON((10 0,0 10,0 0,0 -10))",
+         true);
+
+    test::apply
+        ("pg079spike1",
+         "POLYGON((10 0,0 10,0 0,0 10,0 -10))",
+         false);
+
+    test::apply
+        ("pg079spike2",
+         "POLYGON((10 0,0 10,0 0,0 8,0 -10))",
+         false);
+
+    test::apply
+        ("pg079spike3",
+         "POLYGON((10 0,0 10,0 0,0 11,0 -10))",
+         false);
+
+    test::apply
+        ("pg079reversed",
+         "POLYGON((10 0,0 -10,0 0,0 10))",
+         false);
 }
 
 template <typename Point>