Automatically generated, fast, robust floating-point predicates (WIP)

extensions/example/Jamfile

@@ -13,3 +13,4 @@ project boost-geometry-examples-extensions
     ;
 
 build-project gis ;
+build-project generic_robust_predicates ;

extensions/example/generic_robust_predicates/Jamfile

@@ -0,0 +1,12 @@
+# Boost.Geometry (aka GGL, Generic Geometry Library)
+
+# Use, modification and distribution is subject to the Boost Software License,
+# Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+# http://www.boost.org/LICENSE_1_0.txt)
+
+
+project boost-geometry-example-extensions-generic_robust_predicates
+    :
+    ;
+
+exe static_side_2d : static_side_2d.cpp ;

extensions/example/generic_robust_predicates/static_side_2d.cpp

@@ -0,0 +1,79 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+
+#define BOOST_GEOMETRY_NO_BOOST_TEST
+
+#include <iostream>
+
+#include <boost/geometry/core/cs.hpp>
+#include <boost/geometry/geometries/point.hpp>
+
+#include <boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp>
+
+
+namespace bg = boost::geometry;
+using point = bg::model::point<double, 2, bg::cs::cartesian>;
+
+template <typename CalculationType>
+struct side_robust_with_static_filter
+{
+private:
+    using ct = CalculationType;
+    using expression = bg::detail::generic_robust_predicates::orient2d;
+    using filter = bg::detail::generic_robust_predicates::stage_a_static
+            <
+                expression,
+                ct
+            >;
+    filter m_filter;
+public:
+    side_robust_with_static_filter(ct x_max, ct y_max, ct x_min, ct y_min)
+        : m_filter(x_max, y_max, x_max, y_max, x_max, y_max,
+                   x_min, y_min, x_min, y_min, x_min, y_min) {};
+
+    template
+    <
+        typename P1,
+        typename P2,
+        typename P  
+    >
+    inline int apply(P1 const& p1, P2 const& p2, P const& p) const
+    {
+        int sign = m_filter.apply(bg::get<0>(p1),
+                                  bg::get<1>(p1),
+                                  bg::get<0>(p2),
+                                  bg::get<1>(p2),
+                                  bg::get<0>(p),
+                                  bg::get<1>(p));
+        if(sign != bg::detail::generic_robust_predicates::sign_uncertain)
+        {
+            return sign;
+        }
+        else
+        {
+            //fallback if filter fails.
+            return bg::strategy::side::side_robust<double>::apply(p1, p2, p);
+        }
+    }
+};
+
+int main(int argc, char** argv)
+{
+    point p1(0.0, 0.0);
+    point p2(1.0, 1.0);
+    point p (0.0, 1.0);
+    side_robust_with_static_filter<double> static_strategy(2.0, 2.0, 1.0, 1.0);
+    std::cout << "Side value: " << static_strategy.apply(p1, p2, p) << "\n"; 
+    return 0;
+}

extensions/test/Jamfile

@@ -28,3 +28,4 @@ build-project gis ;
 build-project iterators ;
 build-project nsphere ;
 build-project triangulation ;
+build-project generic_robust_predicates ;

extensions/test/generic_robust_predicates/Jamfile

@@ -0,0 +1,12 @@
+# Boost.Geometry (aka GGL, Generic Geometry Library)
+#
+# Use, modification and distribution is subject to the Boost Software License,
+# Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+# http://www.boost.org/LICENSE_1_0.txt)
+
+test-suite boost-geometry-extensions-generic_robust_predicates
+    :
+    [ run approximate.cpp ]
+    [ run side3d.cpp : : : <debug-symbols>off ]
+    ;
+

extensions/test/generic_robust_predicates/approximate.cpp

@@ -0,0 +1,66 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Unit Test
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <array>
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/approximate.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expression_tree.hpp>
+
+template <typename CalculationType>
+void test_all()
+{
+    using bg::detail::generic_robust_predicates::_1;
+    using bg::detail::generic_robust_predicates::_2;
+    using bg::detail::generic_robust_predicates::_3;
+    using bg::detail::generic_robust_predicates::_4;
+    using bg::detail::generic_robust_predicates::sum;
+    using bg::detail::generic_robust_predicates::difference;
+    using bg::detail::generic_robust_predicates::product;
+    using bg::detail::generic_robust_predicates::max;
+    using bg::detail::generic_robust_predicates::abs;
+    using bg::detail::generic_robust_predicates::post_order;
+    using bg::detail::generic_robust_predicates::approximate_value;
+    using bg::detail::generic_robust_predicates::get_approx;
+    using bg::detail::generic_robust_predicates::approximate_interim;
+    using ct = CalculationType;
+    ct r1 = approximate_value<sum<_1, _2>, ct>(
+            std::array<ct, 2>{1.0, 2.0});
+    BOOST_CHECK_EQUAL(3.0, r1);
+    ct r2 = approximate_value<max<abs<_1>, abs<_2>>, ct>(
+            std::array<ct, 2>{-10.0, 2.0});
+    BOOST_CHECK_EQUAL(10.0, r2);
+
+    using expression = product
+        <
+            difference<_1, _2>,
+            difference<_3, _4>
+        >;
+    using evals = post_order<expression>;
+    std::array<ct, boost::mp11::mp_size<evals>::value> r;
+    std::array<ct, 4> input {5.0, 3.0, 2.0, 8.0};
+    approximate_interim<evals, evals, ct>(r, input);
+    ct r3 = get_approx<evals, typename expression::left, ct>(r, input);
+    BOOST_CHECK_EQUAL(2.0, r3);
+    ct r4 = get_approx<evals, typename expression::right, ct>(r, input);
+    BOOST_CHECK_EQUAL(-6.0, r4);
+    ct r5 = get_approx<evals, expression, ct>(r, input);
+    BOOST_CHECK_EQUAL(-12.0, r5);
+}
+
+
+int test_main(int, char* [])
+{
+    test_all<double>();
+    return 0;
+}

extensions/test/generic_robust_predicates/side3d.cpp

@@ -0,0 +1,76 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Unit Test
+
+// Copyright (c) 2020 Tinko Bartels, Berlin, Germany.
+
+// Contributed and/or modified by Tinko Bartels,
+//   as part of Google Summer of Code 2020 program.
+
+// Use, modification and distribution is subject to the Boost Software License,
+// Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
+// http://www.boost.org/LICENSE_1_0.txt)
+
+#include <array>
+
+#include <geometry_test_common.hpp>
+
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/expressions.hpp>
+#include <boost/geometry/extensions/generic_robust_predicates/strategies/cartesian/detail/stage_a.hpp>
+
+template <typename CalculationType>
+void test_all()
+{
+    using bg::detail::generic_robust_predicates::orient3d;
+    using bg::detail::generic_robust_predicates::stage_a_semi_static;
+    using bg::detail::generic_robust_predicates::stage_a_almost_static;
+    using bg::detail::generic_robust_predicates::stage_a_static;
+    using ct = CalculationType;
+    using semi_static = stage_a_semi_static<orient3d, ct>;
+    BOOST_CHECK_EQUAL(1,
+                      semi_static::apply(1, 0, 0,
+                                         0, 1, 0,
+                                         1, 1, 0,
+                                         0, 0, 1));
+    BOOST_CHECK_EQUAL(-1,
+                      semi_static::apply(1, 0, 0,       
+                                         0, 1, 0,                
+                                         1, 1, 0,                
+                                         0, 0, -1));
+    BOOST_CHECK_EQUAL(0,
+                      semi_static::apply(1, 0, 0,       
+                                         0, 1, 0,                
+                                         1, 1, 0,                
+                                         0, 0, 0));
+    stage_a_static<orient3d, ct> stat(10, 10, 10,
+                                      10, 10, 10,
+                                      10, 10, 10,
+                                      10, 10, 10,
+                                       0,  0,  0,
+                                       0,  0,  0,
+                                       0,  0,  0,
+                                       0,  0,  0);
+    BOOST_CHECK_EQUAL(1, stat.apply(1, 0, 0,
+                                    0, 1, 0,
+                                    1, 1, 0,
+                                    0, 0, 1));
+
+    stage_a_static<orient3d, ct> pessimistic(1e40, 1e40, 1e40,
+                                             1e40, 1e40, 1e40,
+                                             1e40, 1e40, 1e40,
+                                             1e40, 1e40, 1e40, 
+                                             0,    0,    0,  
+                                             0,    0,    0, 
+                                             0,    0,    0,
+                                             0,    0,    0);
+    BOOST_CHECK_EQUAL(bg::detail::generic_robust_predicates::sign_uncertain,
+                      pessimistic.apply(1, 0, 0,
+                                        0, 1, 0,
+                                        1, 1, 0,
+                                        0, 0, 1));
+}
+
+int test_main(int, char* [])
+{
+    test_all<double>();
+    return 0;
+}

include/boost/geometry/algorithms/detail/buffer/buffer_inserter.hpp

@@ -1,6 +1,6 @@
 // Boost.Geometry (aka GGL, Generic Geometry Library)
 
-// Copyright (c) 2012-2014 Barend Gehrels, Amsterdam, the Netherlands.
+// Copyright (c) 2012-2020 Barend Gehrels, Amsterdam, the Netherlands.
 
 // This file was modified by Oracle on 2017.
 // Modifications copyright (c) 2017 Oracle and/or its affiliates.
@@ -31,16 +31,13 @@
 
 #include <boost/geometry/strategies/buffer.hpp>
 #include <boost/geometry/strategies/side.hpp>
-#include <boost/geometry/algorithms/detail/make/make.hpp>
+#include <boost/geometry/algorithms/detail/direction_code.hpp>
 #include <boost/geometry/algorithms/detail/buffer/buffered_piece_collection.hpp>
 #include <boost/geometry/algorithms/detail/buffer/line_line_intersection.hpp>
 
-#include <boost/geometry/algorithms/assign.hpp>
 #include <boost/geometry/algorithms/num_interior_rings.hpp>
 #include <boost/geometry/algorithms/simplify.hpp>
 
-#include <boost/geometry/arithmetic/infinite_line_functions.hpp>
-
 #include <boost/geometry/views/detail/normalized_view.hpp>
 
 
@@ -113,23 +110,11 @@ struct buffer_range
             JoinStrategy const& join_strategy,
             EndStrategy const& end_strategy,
             RobustPolicy const& ,
-            SideStrategy const& side_strategy) // side strategy
+            SideStrategy const& side_strategy)
     {
-        output_point_type intersection_point;
-        geometry::assign_zero(intersection_point);
-
-        geometry::strategy::buffer::join_selector join
-                = get_join_type(penultimate_input, previous_input, input, side_strategy);
-        if (join == geometry::strategy::buffer::join_convex)
-        {
-            // Calculate the intersection-point formed by the two sides.
-            // It might be that the two sides are not convex, but continue
-            // or spikey, we then change the join-type
-            join = line_line_intersection::apply(
-                        perp1, perp2, prev_perp1, prev_perp2,
-                        intersection_point);
-
-        }
+        geometry::strategy::buffer::join_selector const join
+                = get_join_type(penultimate_input, previous_input, input,
+                                side_strategy);
 
         switch(join)
         {
@@ -163,6 +148,11 @@ struct buffer_range
                 {
                     // The corner is convex, we create a join
                     // TODO (future) - avoid a separate vector, add the piece directly
+                    output_point_type const
+                        intersection_point
+                            = line_line_intersection::apply(perp1, perp2,
+                                    prev_perp1, prev_perp2);
+
                     std::vector<output_point_type> range_out;
                     if (join_strategy.apply(intersection_point,
                                 previous_input, prev_perp2, perp1,
@@ -177,14 +167,14 @@ struct buffer_range
         }
     }
 
-    static inline bool similar_direction(output_point_type const& p0,
+    // Returns true if collinear point p2 continues after p0 and p1.
+    // If it turns back (spike), it returns false.
+    static inline bool same_direction(output_point_type const& p0,
             output_point_type const& p1,
             output_point_type const& p2)
     {
-        typedef model::infinite_line<coordinate_type> line_type;
-        line_type const p = detail::make::make_infinite_line<coordinate_type>(p0, p1);
-        line_type const q = detail::make::make_infinite_line<coordinate_type>(p1, p2);
-        return arithmetic::similar_direction(p, q);
+        typedef typename cs_tag<output_point_type>::type cs_tag;
+        return direction_code<cs_tag>(p0, p1, p2) == 1;
     }
 
     template <typename SideStrategy>
@@ -197,8 +187,7 @@ struct buffer_range
         int const side = side_strategy.apply(p0, p1, p2);
         return side == -1 ? geometry::strategy::buffer::join_convex
             :  side == 1  ? geometry::strategy::buffer::join_concave
-            :  similar_direction(p0, p1, p2)
-                          ? geometry::strategy::buffer::join_continue
+            :  same_direction(p0, p1, p2) ? geometry::strategy::buffer::join_continue
             : geometry::strategy::buffer::join_spike;
     }
 

include/boost/geometry/algorithms/detail/buffer/buffer_policies.hpp

@@ -139,8 +139,8 @@ struct buffer_turn_operation
     {}
 };
 
-// Version for buffer including type of location, is_opposite, and helper variables
-template <typename Point, typename RobustPoint, typename SegmentRatio>
+// Version of turn_info for buffer with its turn index and other helper variables
+template <typename Point, typename SegmentRatio>
 struct buffer_turn_info
     : public detail::overlay::turn_info
         <
@@ -153,9 +153,6 @@ struct buffer_turn_info
 
     std::size_t turn_index;
 
-    // Still necessary for turn-in-original calculations
-    RobustPoint robust_point;
-
     // Information if turn can be used. It is not traversable if it is within
     // another piece, or within the original (depending on deflation),
     // or (for deflate) if there are not enough points to traverse it.