[coordinate_types] support for Boost.Rational and Boost.Multiprecision

including a unit test

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

@@ -79,7 +79,8 @@ struct direction_code_impl<cartesian_tag>
         }
 
         calc_t const sv = arithmetic::side_value(line, point);
-        return sv == 0 ? 0 : sv > 0 ? 1 : -1;
+        static calc_t const zero = 0;
+        return sv == zero ? 0 : sv > zero ? 1 : -1;
     }
 };
 

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

@@ -23,7 +23,7 @@ namespace detail { namespace overlay
 
 template <typename Point1, typename Point2, typename E>
 inline bool approximately_equals(Point1 const& a, Point2 const& b,
-                                 E const& multiplier)
+                                 E const& epsilon_multiplier)
 {
     using coor_t = typename select_coordinate_type<Point1, Point2>::type;
     using calc_t = typename geometry::select_most_precise<coor_t, E>::type;
@@ -34,7 +34,7 @@ inline bool approximately_equals(Point1 const& a, Point2 const& b,
     calc_t const& b1 = geometry::get<1>(b);
 
     math::detail::equals_factor_policy<calc_t> policy(a0, b0, a1, b1);
-    policy.factor *= multiplier;
+    policy.multiply_epsilon(epsilon_multiplier);
 
     return math::detail::equals_by_policy(a0, b0, policy)
         && math::detail::equals_by_policy(a1, b1, policy);

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

@@ -39,7 +39,8 @@ struct sweep_equal_policy
     static inline bool equals(P const& p1, P const& p2)
     {
         // Points within a kilo epsilon are considered as equal
-        return approximately_equals(p1, p2, 1000.0);
+        using coor_t = typename coordinate_type<P>::type;
+        return approximately_equals(p1, p2, coor_t(1000));
     }
 
     template <typename T>
@@ -47,7 +48,8 @@ struct sweep_equal_policy
     {
         // This threshold is an arbitrary value
         // as long as it is than the used kilo-epsilon
-        return value > 1.0e-3;
+        T const limit = T(1) / T(1000);
+        return value > limit;
     }
 };
 

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

@@ -128,7 +128,8 @@ struct base_turn_handler
         }
 
         auto 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;
+        static decltype(dm.measure) const zero = 0;
+        return dm.measure == zero ? 0 : dm.measure > zero ? 1 : -1;
     }
 
     template

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

@@ -316,7 +316,15 @@ public :
         // then take a point (or more) further back.
         // The limit of offset avoids theoretical infinite loops.
         // In practice it currently walks max 1 point back in all cases.
-        double const tolerance = 1.0e9;
+        // Use the coordinate type, but if it is too small (e.g. std::int16), use a double
+        using ct_type = typename geometry::select_most_precise
+            <
+                typename geometry::coordinate_type<Point>::type,
+                double
+            >::type;
+
+        ct_type const tolerance = 1000000000;
+
         int offset = 0;
         while (approximately_equals(point_from, turn.point, tolerance)
                && offset > -10)

include/boost/geometry/algorithms/detail/sections/sectionalize.hpp

@@ -783,22 +783,24 @@ inline void enlarge_sections(Sections& sections, Strategy const&)
 
     // It makes section a tiny bit too large, which might cause (a small number)
     // of more comparisons
-    for (typename boost::range_iterator<Sections>::type it = boost::begin(sections);
+    for (auto& section : sections)
-        it != boost::end(sections);
-        ++it)
     {
 #if defined(BOOST_GEOMETRY_USE_RESCALING)
         detail::sectionalize::expand_by_epsilon
             <
                 typename Strategy::cs_tag
-            >::apply(it->bounding_box);
+            >::apply(section.bounding_box);
 
 #else
         // Expand the box to avoid missing any intersection. The amount is
         // should be larger than epsilon. About the value itself: the smaller
         // it is, the higher the risk to miss intersections. The larger it is,
         // the more comparisons are made. So it should be on the high side.
-        detail::buffer::buffer_box(it->bounding_box, 0.001, it->bounding_box);
+        using gt = decltype(section.bounding_box);
+        using ct = typename geometry::coordinate_type<gt>::type;
+        ct const tolerance = ct(1) / ct(1000);
+        detail::buffer::buffer_box(section.bounding_box, tolerance,
+                                   section.bounding_box);
 #endif
     }
 }

include/boost/geometry/policies/robustness/segment_ratio.hpp

@@ -119,7 +119,8 @@ struct possibly_collinear<Type, false>
     template <typename Ratio, typename Threshold>
     static inline bool apply(Ratio const& ratio, Threshold)
     {
-        return ratio.denominator() == 0;
+        static Type const zero = 0;
+        return ratio.denominator() == zero;
     }
 };
 
@@ -216,14 +217,14 @@ public :
     {
         // Minimal normalization
         // 1/-4 => -1/4, -1/-4 => 1/4
-        if (m_denominator < 0)
+        if (m_denominator < zero_instance())
         {
             m_numerator = -m_numerator;
             m_denominator = -m_denominator;
         }
 
         m_approximation =
-            m_denominator == 0 ? 0
+            m_denominator == zero_instance() ? 0
             : (
                 boost::numeric_cast<floating_point_type>(m_numerator) * scale()
                 / boost::numeric_cast<floating_point_type>(m_denominator)
@@ -235,12 +236,12 @@ public :
     inline bool on_segment() const
     {
         // e.g. 0/4 or 4/4 or 2/4
-        return m_numerator >= 0 && m_numerator <= m_denominator;
+        return m_numerator >= zero_instance() && m_numerator <= m_denominator;
     }
     inline bool in_segment() const
     {
         // e.g. 1/4
-        return m_numerator > 0 && m_numerator < m_denominator;
+        return m_numerator > zero_instance() && m_numerator < m_denominator;
     }
     inline bool on_end() const
     {
@@ -250,7 +251,7 @@ public :
     inline bool left() const
     {
         // e.g. -1/4
-        return m_numerator < 0;
+        return m_numerator < zero_instance();
     }
     inline bool right() const
     {
@@ -336,6 +337,11 @@ private :
         static floating_point_type const fp_scale{1000000.0};
         return fp_scale;
     }
+
+    static inline Type zero_instance()
+    {
+        return 0;
+    }
 };
 
 

include/boost/geometry/strategy/cartesian/side_non_robust.hpp

@@ -87,7 +87,7 @@ public:
                 );
         PromotedType const zero = PromotedType();
 
-        return sv == 0 ? 0 : sv > zero ? 1 : -1;
+        return sv == zero ? 0 : sv > zero ? 1 : -1;
     }
 
 };

include/boost/geometry/util/math.hpp

@@ -139,6 +139,12 @@ struct equals_factor_policy
         return factor;
     }
 
+    template <typename E>
+    void multiply_epsilon(E const& multiplier)
+    {
+        factor *= multiplier;
+    }
+
     T factor;
 };
 
@@ -153,6 +159,8 @@ struct equals_factor_policy<T, false>
     {
         return T(1);
     }
+
+    void multiply_epsilon(T const& ) {}
 };
 
 template <typename Type,

test/algorithms/detail/approximately_equals.cpp

@@ -96,5 +96,10 @@ int test_main(int, char* [])
     test_all<bg::model::point<double, 2, bg::cs::cartesian>>(m, 23);
     test_all<bg::model::point<float, 2, bg::cs::cartesian>>(m, 23);
 
+    // MP is not a floating point type, therefore approximately_equal
+    // is behaves as exact and returns only true if they are identical.
+    // This takes 334 steps (then "d" above is 0.0)
+    test_all<bg::model::point<mp_test_type, 2, bg::cs::cartesian>>(m, 334);
+
     return 0;
 }

test/algorithms/overlay/get_clusters.cpp

@@ -93,7 +93,7 @@ int test_main(int, char* [])
     using fp = bg::model::point<float, 2, bg::cs::cartesian>;
     using dp = bg::model::point<double, 2, bg::cs::cartesian>;
     using ep = bg::model::point<long double, 2, bg::cs::cartesian>;
-    test_get_clusters<fp>(1.0e-8f);
+    test_get_clusters<fp>(1.0e-4);
     test_get_clusters<dp>(1.0e-13);
     test_get_clusters<ep>(1.0e-16);
     return 0;

test/algorithms/overlay/multi_overlay_cases.hpp

@@ -26,6 +26,27 @@ static std::string case_multi_simplex[2] =
     "MULTIPOLYGON(((3 0,0 3,4 5,3 0)))"
 };
 
+// To support exact behavior in integer coordinates (rectangular and diagonal)
+static std::string case_multi_rectangular[2] =
+{
+    "MULTIPOLYGON(((100 100,100 200,200 200,200 100,100 100)),"
+                 "((300 100,300 200,400 200,400 100,300 100),"
+                   "(325 125,375 125,375 175,325 175,325 125)))",
+    "MULTIPOLYGON(((150 50,150 150,350 150,350 50,150 50)))"
+};
+
+static std::string case_multi_diagonal[2] =
+{
+    "MULTIPOLYGON(((40 0,0 40,40 80,80 40,40 0),(50 20,60 30,50 40,40 30,50 20)))",
+    "MULTIPOLYGON(((80 0,40 40,80 80,120 40,80 0),(80 30,90 40,80 50,70 40,80 30)))"
+};
+
+static std::string case_multi_hard[2] =
+{
+    "MULTIPOLYGON(((0 0,0 4,2 4,2 3,4 3,4 0,0 0)))",
+    "MULTIPOLYGON(((2 7,4 7,4 2.99959993362426758,2 3,2 7)))"
+};
+
 // To mix multi/single
 static std::string case_single_simplex = "POLYGON((3 0,0 3,4 5,3 0))";
 

test/algorithms/set_operations/union/Jamfile

@@ -24,5 +24,6 @@ test-suite boost-geometry-algorithms-union
     [ run union_multi.cpp         : : : <define>BOOST_GEOMETRY_TEST_ONLY_ONE_TYPE
                                         : algorithms_union_multi ]
     [ run union_pl_pl.cpp         : : : : algorithms_union_pl_pl ]
-	[ run union_tupled.cpp        : : : : algorithms_union_tupled ]
+    [ run union_tupled.cpp        : : : : algorithms_union_tupled ]
+    [ run union_other_types.cpp   : : : : algorithms_union_other_types ]
     ;

test/algorithms/set_operations/union/union_multi.cpp

@@ -201,6 +201,10 @@ void test_areal()
     TEST_UNION(case_140_multi, 2, 1, -1, 64.953);
     TEST_UNION(case_141_multi, 1, 0, -1, 100.0);
 
+    TEST_UNION(case_multi_rectangular, 1, 1, -1, 33125);
+    TEST_UNION(case_multi_diagonal, 1, 2, -1, 5350);
+    TEST_UNION(case_multi_hard, 1, 0, -1, 22);
+
     test_one<Polygon, MultiPolygon, MultiPolygon>("case_recursive_boxes_1",
         case_recursive_boxes_1[0], case_recursive_boxes_1[1],
         1, 1, 16, 97.0);
@@ -409,7 +413,7 @@ void test_areal()
         ticket_12118[0], ticket_12118[1],
         1, -1, 27, 2221.38713);
 
-#if defined(BOOST_GEOMETRY_TEST_FAILURES) || ! defined(BOOST_GEOMETRY_USE_RESCALING)
+#if ! defined(BOOST_GEOMETRY_USE_RESCALING) || defined(BOOST_GEOMETRY_TEST_FAILURES)
     // No output if rescaling is done
     test_one<Polygon, MultiPolygon, MultiPolygon>("ticket_12125",
         ticket_12125[0], ticket_12125[1],

test/algorithms/set_operations/union/union_other_types.cpp

@@ -0,0 +1,154 @@
+// Boost.Geometry (aka GGL, Generic Geometry Library)
+// Unit Test
+
+// Copyright (c) 2021 Barend Gehrels, Amsterdam, the Netherlands.
+
+// 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_ROBUSTNESS
+
+#include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+
+#include <geometry_test_common.hpp>
+#include <algorithms/overlay/multi_overlay_cases.hpp>
+
+#include <boost/geometry/algorithms/correct.hpp>
+#include <boost/geometry/algorithms/union.hpp>
+#include <boost/geometry/io/wkt/wkt.hpp>
+#include <boost/geometry/geometries/geometries.hpp>
+
+#include <algorithms/overlay/multi_overlay_cases.hpp>
+
+#include <boost/geometry/geometries/point_xy.hpp>
+#include <boost/geometry/geometries/multi_polygon.hpp>
+
+#include <boost/geometry/util/rational.hpp>
+
+#include <set>
+
+enum class exclude { all, rectangular, diagonal, hard, fp };
+
+template <typename Geometry, typename Expected>
+void test_one(std::string const& case_id,
+              std::string const& wkt1, std::string const& wkt2,
+              bool debug,
+              Expected const& expected_area,
+              Expected const& expected_max = -1)
+{
+    using coor_t = typename bg::coordinate_type<Geometry>::type;
+    Geometry g1, g2, clip;
+
+    bg::read_wkt(wkt1, g1);
+    bg::read_wkt(wkt2, g2);
+
+    bg::correct(g1);
+    bg::correct(g2);
+
+    bg::union_(g1, g2, clip);
+
+    auto const area = bg::area(clip);
+    if (debug)
+    {
+      std::cout << "AREA: " << std::setprecision(64) << area
+                << " expected " << expected_area
+                << " types coordinate " << string_from_type<coor_t>::name()
+                << " area " << typeid(decltype(area)).name()
+                << " expected " << typeid(Expected).name()
+                << " size " << sizeof(coor_t)
+                << std::endl;
+    }
+
+    // Check areas, they always have to be specified in integer for this test
+    // and therefore the checking (including a tolerance) is different
+    bool const ok = expected_max == -1
+                    ? bg::math::equals(area, expected_area)
+                    : bg::math::larger_or_equals(area, expected_area)
+                      && bg::math::smaller_or_equals(area, expected_max);
+    BOOST_CHECK_MESSAGE(ok,
+            "union: " << case_id
+            << " area: expected: " << expected_area
+            << " detected: " << area
+            << " type: " << (string_from_type<coor_t>::name())
+            << " (" << (typeid(coor_t).name()) << ")");
+}
+
+template <typename Point>
+void test_areal(std::set<exclude> const& exclude = {}, bool debug = false)
+{
+    using polygon = bg::model::polygon<Point>;
+    using multi_polygon = bg::model::multi_polygon<polygon>;
+
+    // Intended tests: only 3:
+    // - simple case having only horizontal/vertical lines ("rectangular")
+    // - simple case on integer grid but also having diagonals ("diagonal")
+    // - case going wrong for <float> ("hard")
+
+    if (exclude.count(exclude::rectangular)
+        + exclude.count(exclude::all) == 0)
+    {
+        test_one<multi_polygon>("case_multi_rectangular",
+            case_multi_rectangular[0], case_multi_rectangular[1], debug, 33125);
+    }
+    if (exclude.count(exclude::diagonal)
+        + exclude.count(exclude::all) == 0)
+    {
+        test_one<multi_polygon>("case_multi_diagonal",
+            case_multi_diagonal[0], case_multi_diagonal[1], debug, 5350);
+    }
+    if (exclude.count(exclude::hard)
+        + exclude.count(exclude::fp)
+        + exclude.count(exclude::all) == 0)
+    {
+        test_one<multi_polygon>("case_multi_hard",
+            case_multi_hard[0], case_multi_hard[1], debug, 21, 23);
+    }
+}
+
+int test_main(int, char* [])
+{
+    namespace bm = boost::multiprecision;
+
+    using bg::model::d2::point_xy;
+
+    // Standard floating point types
+    test_areal<point_xy<float>>({exclude::hard});
+    test_areal<point_xy<double>>({});
+    test_areal<point_xy<long double>>({});
+
+    // Standard integer types
+    test_areal<point_xy<std::int16_t>>({exclude::fp});
+    test_areal<point_xy<std::int32_t>>({exclude::fp});
+    test_areal<point_xy<std::int64_t>>({exclude::fp});
+
+    // Boost multi precision (integer)
+    test_areal<point_xy<bm::int128_t>>({exclude::fp});
+    test_areal<point_xy<bm::checked_int128_t>>({exclude::fp});
+
+    // Boost multi precision (floating point)
+    test_areal<point_xy<bm::number<bm::cpp_bin_float<5>>>>();
+    test_areal<point_xy<bm::number<bm::cpp_bin_float<10>>>>();
+    test_areal<point_xy<bm::number<bm::cpp_bin_float<50>>>>();
+    test_areal<point_xy<bm::number<bm::cpp_bin_float<100>>>>();
+
+    test_areal<point_xy<bm::number<bm::cpp_dec_float<50>>>>({});
+
+    // Boost multi precision (rational)
+    test_areal<point_xy<bm::cpp_rational>>({exclude::fp});
+    test_areal<point_xy<bm::checked_cpp_rational>>({exclude::fp});
+
+    // Boost multi precision float128 wrapper, is currently NOT supported
+    // and it is limited to certain compilers anyway
+    // test_areal<point_xy<bm::float128>>();
+
+    // Boost rational (tests compilation)
+    // (the rectangular case is correct; other input might give wrong results)
+    // The int16 version throws a <zero denominator> exception
+    test_areal<point_xy<boost::rational<std::int16_t>>>({exclude::all});
+    test_areal<point_xy<boost::rational<std::int32_t>>>({exclude::fp});
+    test_areal<point_xy<boost::rational<std::int64_t>>>({exclude::fp});
+
+    return 0;
+}

test/robustness/overlay/areal_areal/general_intersection_precision.cpp

@@ -38,6 +38,7 @@ static std::string case_b[2] =
     "MULTIPOLYGON(((-1 -1,-1 8,8 8,8 -1,-1 -1),(2 7,2 3,4 3,4 7,2 7)))"
     };
 
+// Union should deliver 14.0
 static std::string case_c[2] =
     {
     "MULTIPOLYGON(((0 0,0 4,2 4,2 3,4 3,4 0,0 0)))",
@@ -49,7 +50,7 @@ struct test_settings
     bool verbose{false};
     bool do_output{false};
 
-    // Settings currently not modifiable, and still giving quite some errors
+    // Settings currently not modifiable
     double start_bound{1.0e-2};
     double step_factor{50.0};  // on each side -> 100 steps per factor
     int max_factor{10000};
@@ -61,7 +62,6 @@ bool test_overlay(std::string const& caseid,
         double expected_area,
         test_settings const& settings)
 {
-
     typedef typename boost::range_value<Geometry>::type geometry_out;
     typedef bg::detail::overlay::overlay
         <
@@ -96,8 +96,8 @@ bool test_overlay(std::string const& caseid,
     overlay::apply(g1, g2, robust_policy, std::back_inserter(result),
                    strategy, visitor);
 
-    const double detected_area = bg::area(result);
+    auto const detected_area = bg::area(result);
-    if (std::fabs(detected_area - expected_area) > 0.01)
+    if (std::fabs(detected_area - expected_area) > 0.1)
     {
         if (settings.do_output)
         {

test/string_from_type.hpp

@@ -24,7 +24,8 @@
 #include <boost/multiprecision/cpp_bin_float.hpp>
 
 template <typename T>
-struct string_from_type {};
+struct string_from_type
+{ static std::string name() { return "?"; }  };
 
 template <> struct string_from_type<void>
 { static std::string name() { return "v"; }  };
@@ -47,7 +48,8 @@ template <> struct string_from_type<int>
 template <> struct string_from_type<long>
 { static std::string name() { return "l"; }  };
 
-template <> struct string_from_type<boost::multiprecision::cpp_bin_float_100>
+template <typename Backend>
+struct string_from_type<boost::multiprecision::number<Backend>>
 { static std::string name() { return "m"; }  };
 
 // this is what g++ and clang++ use