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>
test_overlay(std::string const & caseid,Geometry const & g1,Geometry const & g2,double expected_area,bool do_output)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>
update(Ring & ring,double x,double y,std::size_t index)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>
test_case(std::size_t & error_count,std::size_t case_index,std::size_t i,std::size_t j,std::size_t min_vertex_index,std::size_t max_vertex_index,double x,double y,double expectation,MultiPolygon const & poly1,MultiPolygon const & poly2,bool do_output)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>
test_all(std::size_t case_index,std::size_t min_vertex_index,std::size_t max_vertex_index,double expectation,bool do_output)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
test_main(int argc,char ** argv)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
225