Merge pull request #749 from vissarion/feature/robust_area

Area strategy for more accurate computations in cartesian CS

Author: vissarion

include/boost/geometry/extensions/triangulation/strategies/cartesian/in_circle_robust.hpp

@@ -12,7 +12,7 @@
 #ifndef BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_IN_CIRCLE_ROBUST_HPP
 #define BOOST_GEOMETRY_EXTENSIONS_TRIANGULATION_STRATEGIES_CARTESIAN_IN_CIRCLE_ROBUST_HPP
 
-#include<boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp>
+#include<boost/geometry/util/precise_math.hpp>
 
 namespace boost { namespace geometry
 {

include/boost/geometry/extensions/triangulation/strategies/cartesian/side_robust.hpp

@@ -14,7 +14,7 @@
 
 #include <boost/geometry/util/select_most_precise.hpp>
 #include <boost/geometry/util/select_calculation_type.hpp>
-#include <boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp>
+#include <boost/geometry/util/precise_math.hpp>
 
 namespace boost { namespace geometry
 {

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

@@ -112,8 +112,7 @@ public :
     }
 
     template <typename Geometry>
-    static inline typename result_type<Geometry>::type
-        result(state<Geometry>& st)
+    static inline auto result(state<Geometry>& st)
     {
         return st.area();
     }

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

@@ -0,0 +1,117 @@
+// Boost.Geometry
+
+// Copyright (c) 2020, Oracle and/or its affiliates.
+
+// Contributed and/or modified by Vissarion Fysikopoulos, on behalf of Oracle
+
+// Licensed under the Boost Software License version 1.0.
+// http://www.boost.org/users/license.html
+
+#ifndef BOOST_GEOMETRY_STRATEGY_CARTESIAN_AREA_ACCURATE_HPP
+#define BOOST_GEOMETRY_STRATEGY_CARTESIAN_AREA_ACCURATE_HPP
+
+#include <boost/mpl/if.hpp>
+
+//#include <boost/geometry/arithmetic/determinant.hpp>
+#include <boost/geometry/core/access.hpp>
+#include <boost/geometry/core/coordinate_type.hpp>
+#include <boost/geometry/core/coordinate_dimension.hpp>
+#include <boost/geometry/strategy/area.hpp>
+#include <boost/geometry/util/select_most_precise.hpp>
+#include <boost/geometry/util/precise_math.hpp>
+
+namespace boost { namespace geometry
+{
+
+namespace strategy { namespace area
+{
+
+/*!
+\brief Cartesian area calculation
+\ingroup strategies
+\details Calculates cartesian area using the trapezoidal rule and precise
+         summation (useful to increase precision with floating point arithmetic)
+\tparam CalculationType \tparam_calculation
+
+\qbk{
+[heading See also]
+[link geometry.reference.algorithms.area.area_2_with_strategy area (with strategy)]
+}
+
+*/
+template
+<
+    typename CalculationType = void
+>
+class precise_cartesian
+{
+public :
+    template <typename Geometry>
+    struct result_type
+        : strategy::area::detail::result_type
+            <
+                Geometry,
+                CalculationType
+            >
+    {};
+
+    template <typename Geometry>
+    class state
+    {
+        friend class precise_cartesian;
+
+        typedef typename result_type<Geometry>::type return_type;
+
+    public:
+        inline state()
+            : sum1(0)
+            , sum2(0)
+        {
+            // Strategy supports only 2D areas
+            assert_dimension<Geometry, 2>();
+        }
+
+    private:
+        inline return_type area() const
+        {
+            return_type const two = 2;
+            return (sum1 + sum2) / two;
+        }
+
+        return_type sum1;
+        return_type sum2;
+    };
+
+    template <typename PointOfSegment, typename Geometry>
+    static inline void apply(PointOfSegment const& p1,
+                             PointOfSegment const& p2,
+                             state<Geometry>& st)
+    {
+        typedef typename state<Geometry>::return_type return_type;
+
+        auto const det = (return_type(get<0>(p1)) + return_type(get<0>(p2)))
+                * (return_type(get<1>(p1)) - return_type(get<1>(p2)));
+
+        auto const res = boost::geometry::detail::precise_math::two_sum(st.sum1, det);
+
+        st.sum1 = res[0];
+        st.sum2 += res[1];
+    }
+
+    template <typename Geometry>
+    static inline auto result(state<Geometry>& st)
+    {
+        return st.area();
+    }
+
+};
+
+
+}} // namespace strategy::area
+
+
+
+}} // namespace boost::geometry
+
+
+#endif // BOOST_GEOMETRY_STRATEGIES_CARTESIAN_AREA_ACCURATE_HPP

include/boost/geometry/extensions/triangulation/strategies/cartesian/detail/precise_math.hpp renamed to include/boost/geometry/util/precise_math.hpp

@@ -27,10 +27,13 @@
 #include <boost/geometry/geometries/ring.hpp>
 #include <boost/geometry/geometries/polygon.hpp>
 
+#include <boost/geometry/strategy/cartesian/precise_area.hpp>
+
 #include <test_geometries/all_custom_ring.hpp>
 #include <test_geometries/all_custom_polygon.hpp>
 //#define BOOST_GEOMETRY_TEST_DEBUG
 
+#include <boost/multiprecision/cpp_dec_float.hpp>
 #include <boost/variant/variant.hpp>
 
 template <typename Polygon>
@@ -175,6 +178,70 @@ void test_large_integers()
     BOOST_CHECK_CLOSE(int_area, double_area, 0.0001);
 }
 
+struct precise_cartesian : bg::strategies::detail::cartesian_base
+{
+    template <typename Geometry>
+    static auto area(Geometry const&)
+    {
+        return bg::strategy::area::precise_cartesian<>();
+    }
+};
+
+void test_accurate_sum_strategy()
+{
+    typedef bg::model::point<double, 2, bg::cs::cartesian> point_type;
+    typedef bg::model::point
+        <
+            boost::multiprecision::cpp_dec_float_50,
+            2,
+            bg::cs::cartesian
+        > mp_point_type;
+
+    auto const poly0_string = "POLYGON((0 0,0 1,1 0,0 0))";
+
+    bg::model::polygon<point_type> poly0;
+    bg::read_wkt(poly0_string, poly0);
+
+    BOOST_CHECK_CLOSE(bg::area(poly0), 0.5, 0.0001);
+    BOOST_CHECK_CLOSE(bg::area(poly0, precise_cartesian()), 0.5, 0.0001);
+
+    bg::model::polygon<mp_point_type> mp_poly0;
+    bg::read_wkt(poly0_string, mp_poly0);
+
+    BOOST_CHECK_CLOSE(bg::area(mp_poly0), 0.5, 0.0001);
+
+    auto const poly1_string = "POLYGON((0.10000000000000001 0.10000000000000001,\
+            0.20000000000000001 0.20000000000000004,\
+            0.79999999999999993 0.80000000000000004,\
+            1.267650600228229e30 1.2676506002282291e30,\
+            0.10000000000000001 0.10000000000000001))";
+
+    bg::model::polygon<point_type> poly1;
+    bg::read_wkt(poly1_string, poly1);
+
+    BOOST_CHECK_CLOSE(bg::area(poly1), 0, 0.0001);
+    BOOST_CHECK_CLOSE(bg::area(poly1, precise_cartesian()), -0.315, 0.0001);
+
+    bg::model::polygon<mp_point_type> mp_poly1;
+    bg::read_wkt(poly1_string, mp_poly1);
+
+    BOOST_CHECK_CLOSE(bg::area(mp_poly1), 34720783012552.6, 0.0001);
+
+    auto const poly2_string = "POLYGON((1.267650600228229e30 1.2676506002282291e30,\
+            0.8 0.8,0.2 0.2,0.1 0.1,1.267650600228229e30 1.2676506002282291e30))";
+
+    bg::model::polygon<point_type> poly2;
+    bg::read_wkt(poly2_string, poly2);
+
+    BOOST_CHECK_CLOSE(bg::area(poly2), 0, 0.0001);
+    BOOST_CHECK_CLOSE(bg::area(poly2, precise_cartesian()), 0.315, 0.0001);
+
+    bg::model::polygon<mp_point_type> mp_poly2;
+    bg::read_wkt(poly2_string, mp_poly2);
+
+    BOOST_CHECK_CLOSE(bg::area(mp_poly2), 35000000000000, 0.0001);
+}
+
 void test_variant()
 {
     typedef bg::model::point<double, 2, bg::cs::cartesian> double_point_type;
@@ -231,6 +298,8 @@ int test_main(int, char* [])
 
     test_variant();
 
+    test_accurate_sum_strategy();
+
     // test_empty_input<bg::model::d2::point_xy<int> >();
 
     return 0;