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