ceph/src/boost/libs/geometry/test/robustness/overlay/areal_areal/general_intersection_precision.cpp

   1 // Boost.Geometry
   2 // Robustness Test
   3
   4 // Copyright (c) 2019 Barend Gehrels, Amsterdam, the Netherlands.
   5
   6 // Use, modification and distribution is subject to the Boost Software License,
   7 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
   8 // http://www.boost.org/LICENSE_1_0.txt)
   9
  10 #include <iostream>
  11 #include <iomanip>
  12 #include <fstream>
  13 #include <sstream>
  14 #include <string>
  15
  16 #include <boost/type_traits/is_same.hpp>
  17
  18 #include <geometry_test_common.hpp>
  19
  20 #include <boost/geometry.hpp>
  21 #include <boost/geometry/geometries/geometries.hpp>
  22
  23 // Basic case. Union should deliver 22.0
  24 static std::string case_a[2] =
  25     {
  26     "MULTIPOLYGON(((0 0,0 4,2 4,2 3,4 3,4 0,0 0)))",
  27     "MULTIPOLYGON(((2 7,4 7,4 3,2 3,2 7)))"
  28     };
  29
  30 // Case with an interior ring. Union should deliver 73.0
  31 static std::string case_b[2] =
  32     {
  33     "MULTIPOLYGON(((0 0,0 4,2 4,2 3,4 3,4 0,0 0)))",
  34     "MULTIPOLYGON(((-1 -1,-1 8,8 8,8 -1,-1 -1),(2 7,2 3,4 3,4 7,2 7)))"
  35     };
  36
  37 static std::string case_c[2] =
  38     {
  39     "MULTIPOLYGON(((0 0,0 4,2 4,2 3,4 3,4 0,0 0)))",
  40     "MULTIPOLYGON(((1 1,0 1,0 3,1 3,1 1)))"
  41     };
  42
  43 template <bg::overlay_type OverlayType, typename Geometry>
  44 bool test_overlay(std::string const& caseid,
  45         Geometry const& g1, Geometry const& g2,
  46         double expected_area,
  47         bool do_output)
  48 {
  49
  50     typedef typename boost::range_value<Geometry>::type geometry_out;
  51     typedef bg::detail::overlay::overlay
  52         <
  53             Geometry, Geometry,
  54             bg::detail::overlay::do_reverse<bg::point_order<Geometry>::value>::value,
  55             OverlayType == bg::overlay_difference
  56             ? ! bg::detail::overlay::do_reverse<bg::point_order<Geometry>::value>::value
  57             : bg::detail::overlay::do_reverse<bg::point_order<Geometry>::value>::value,
  58             bg::detail::overlay::do_reverse<bg::point_order<Geometry>::value>::value,
  59             geometry_out,
  60             OverlayType
  61         > overlay;
  62
  63     typedef typename bg::strategy::intersection::services::default_strategy
  64         <
  65             typename bg::cs_tag<Geometry>::type
  66         >::type strategy_type;
  67
  68     strategy_type strategy;
  69
  70     typedef typename bg::rescale_overlay_policy_type
  71     <
  72         Geometry,
  73         Geometry
  74     >::type rescale_policy_type;
  75
  76     rescale_policy_type robust_policy
  77         = bg::get_rescale_policy<rescale_policy_type>(g1, g2);
  78
  79     Geometry result;
  80     bg::detail::overlay::overlay_null_visitor visitor;
  81     overlay::apply(g1, g2, robust_policy, std::back_inserter(result),
  82                    strategy, visitor);
  83
  84     const double detected_area = bg::area(result);
  85     if (std::fabs(detected_area - expected_area) > 0.01)
  86     {
  87         if (do_output)
  88         {
  89             std::cout << "ERROR: " << caseid << std::setprecision(18)
  90                       << " detected=" << detected_area
  91                       << " expected=" << expected_area << std::endl
  92                       << "    " <<  bg::wkt(g1) << std::endl
  93                       << "    " << bg::wkt(g2) << std::endl;
  94         }
  95         return false;
  96     }
  97     return true;
  98 }
  99
 100 template <typename Ring>
 101 void update(Ring& ring, double x, double y, std::size_t index)
 102 {
 103     if (index >= ring.size())
 104     {
 105         return;
 106     }
 107     bg::set<0>(ring[index], bg::get<0>(ring[index]) + x);
 108     bg::set<1>(ring[index], bg::get<1>(ring[index]) + y);
 109     if (index == 0)
 110     {
 111         ring.back() = ring.front();
 112     }
 113 }
 114
 115 template <bg::overlay_type OverlayType, typename MultiPolygon>
 116 std::size_t test_case(std::size_t& error_count,
 117         std::size_t case_index, std::size_t i, std::size_t j,
 118         std::size_t min_vertex_index, std::size_t max_vertex_index,
 119         double x, double y, double expectation,
 120         MultiPolygon const& poly1, MultiPolygon const& poly2,
 121         bool do_output)
 122 {
 123     std::size_t n = 0;
 124     for (std::size_t k = min_vertex_index; k <= max_vertex_index; k++, ++n)
 125     {
 126         MultiPolygon poly2_adapted = poly2;
 127
 128         switch (case_index)
 129         {
 130             case 2 :
 131                 update(bg::interior_rings(poly2_adapted.front()).front(), x, y, k);
 132                 break;
 133             default :
 134                 update(bg::exterior_ring(poly2_adapted.front()), x, y, k);
 135                 break;
 136         }
 137
 138         std::ostringstream out;
 139         out << "case_" << i << "_" << j << "_" << k;
 140         if (! test_overlay<OverlayType>(out.str(), poly1, poly2_adapted, expectation, do_output))
 141         {
 142             if (error_count == 0)
 143             {
 144                 // First failure is always reported
 145                 test_overlay<OverlayType>(out.str(), poly1, poly2_adapted, expectation, true);
 146             }
 147             error_count++;
 148         }
 149     }
 150     return n;
 151 }
 152
 153 template <typename T, bool Clockwise, bg::overlay_type OverlayType>
 154 std::size_t test_all(std::size_t case_index, std::size_t min_vertex_index,
 155                      std::size_t max_vertex_index,
 156                      double expectation, bool do_output)
 157 {
 158     typedef bg::model::point<T, 2, bg::cs::cartesian> point_type;
 159     typedef bg::model::polygon<point_type, Clockwise> polygon;
 160     typedef bg::model::multi_polygon<polygon> multi_polygon;
 161
 162     const std::string& first = case_a[0];
 163
 164     const std::string& second
 165             = case_index == 1 ? case_a[1]
 166             : case_index == 2 ? case_b[1]
 167             : case_index == 3 ? case_c[1]
 168             : "";
 169
 170     multi_polygon poly1;
 171     bg::read_wkt(first, poly1);
 172
 173     multi_polygon poly2;
 174     bg::read_wkt(second, poly2);
 175
 176     std::size_t error_count = 0;
 177     std::size_t n = 0;
 178     for (double factor = 1.0; factor > 1.0e-10; factor /= 10.0)
 179     {
 180         std::size_t i = 0;
 181         double const bound = 1.0e-5 * factor;
 182         double const step = 1.0e-6 * factor;
 183         for (double x = -bound; x <= bound; x += step, ++i)
 184         {
 185             std::size_t j = 0;
 186             for (double y = -bound; y <= bound; y += step, ++j, ++n)
 187             {
 188                 n += test_case<OverlayType>(error_count,
 189                                             case_index, i, j,
 190                                        min_vertex_index, max_vertex_index,
 191                                        x, y, expectation,
 192                                        poly1, poly2, do_output);
 193             }
 194         }
 195     }
 196
 197     std::cout << case_index
 198             << " #cases: " << n << " #errors: " << error_count << std::endl;
 199     BOOST_CHECK_EQUAL(error_count, 0u);
 200
 201     return error_count;
 202 }
 203
 204 int test_main(int argc, char** argv)
 205 {
 206     typedef double coordinate_type;
 207
 208     const double factor = argc > 1 ? atof(argv[1]) : 2.0;
 209     const bool do_output = argc > 2 && atol(argv[2]) == 1;
 210
 211     std::cout << std::endl << " -> TESTING "
 212               << string_from_type<coordinate_type>::name()
 213               << " " << factor
 214               << std::endl;
 215
 216     test_all<coordinate_type, true, bg::overlay_union>(1, 0, 3, 22.0, do_output);
 217     test_all<coordinate_type, true, bg::overlay_union>(2, 0, 3, 73.0, do_output);
 218     test_all<coordinate_type, true, bg::overlay_intersection>(3, 1, 2, 2.0, do_output);
 219     test_all<coordinate_type, true, bg::overlay_union>(3, 1, 2, 14.0, do_output);
 220
 221     return 0;
 222  }
 223
 224