1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 // Unit Test
3 //
4 // Copyright (c) 2009-2015 Barend Gehrels, Amsterdam, the Netherlands.
5 // Use, modification and distribution is subject to the Boost Software License,
6 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
7 // http://www.boost.org/LICENSE_1_0.txt)
8
9 #ifndef BOOST_GEOMETRY_TEST_OVERLAY_P_Q_HPP
10 #define BOOST_GEOMETRY_TEST_OVERLAY_P_Q_HPP
11
12 #include <iostream>
13 #include <fstream>
14 #include <sstream>
15 #include <iomanip>
16
17 #include <boost/typeof/typeof.hpp>
18
19 //#define BOOST_GEOMETRY_ROBUSTNESS_USE_DIFFERENCE
20
21 #include <geometry_test_common.hpp>
22
23 // For mixing int/float
24 #if defined(_MSC_VER)
25 #pragma warning( disable : 4244 )
26 #pragma warning( disable : 4267 )
27 #endif
28
29
30 #include <boost/geometry.hpp>
31 #include <boost/geometry/geometries/geometries.hpp>
32 #include <boost/geometry/geometries/point_xy.hpp>
33 #include <boost/geometry/io/svg/svg_mapper.hpp>
34
35 #include <boost/geometry/algorithms/detail/overlay/debug_turn_info.hpp>
36 #include <boost/geometry/algorithms/intersects.hpp>
37 #include <boost/geometry/algorithms/is_valid.hpp>
38 #include <boost/geometry/algorithms/touches.hpp>
39
40 struct p_q_settings
41 {
42 bool svg;
43 bool also_difference;
44 bool validity;
45 bool wkt;
46 bool verify_area;
47 double tolerance;
48
p_q_settingsp_q_settings49 p_q_settings()
50 : svg(false)
51 , also_difference(false)
52 , validity(false)
53 , wkt(false)
54 , verify_area(false)
55 , tolerance(1.0e-3) // since rescaling to integer the tolerance should be less. Was originally 1.0e-6
56 {}
57 };
58
59 template <typename Geometry>
p_q_area(Geometry const & g)60 inline typename bg::default_area_result<Geometry>::type p_q_area(Geometry const& g)
61 {
62 try
63 {
64 return bg::area(g);
65 }
66 catch(bg::empty_input_exception const&)
67 {
68 return 0;
69 }
70 }
71
72 struct verify_area
73 {
74 template <typename Iterator>
check_ringverify_area75 static inline bool check_ring(Iterator begin, Iterator end)
76 {
77 for (Iterator it = begin; it != end; ++it)
78 {
79 double const area = bg::area(*it);
80 if (fabs(area) < 0.01)
81 {
82 return false;
83 }
84 }
85 return true;
86 }
87
88 template <typename Interiors>
check_ringsverify_area89 static inline bool check_rings(Interiors const& rings)
90 {
91 return check_ring(boost::begin(rings), boost::end(rings));
92 }
93
94 template <typename Iterator>
check_polysverify_area95 static inline bool check_polys(Iterator begin, Iterator end)
96 {
97 for (Iterator it = begin; it != end; ++it)
98 {
99 // If necessary, exterior_ring can be checked too
100 if (! check_rings(bg::interior_rings(*it)))
101 {
102 return false;
103 }
104 }
105 return true;
106 }
107
108 template <typename Geometry>
applyverify_area109 static inline bool apply(Geometry const& g)
110 {
111 return check_polys(boost::begin(g), boost::end(g));
112 }
113 };
114
115 template <typename OutputType, typename CalculationType, typename G1, typename G2>
test_overlay_p_q(std::string const & caseid,G1 const & p,G2 const & q,p_q_settings const & settings)116 static bool test_overlay_p_q(std::string const& caseid,
117 G1 const& p, G2 const& q,
118 p_q_settings const& settings)
119 {
120 bool result = true;
121
122 typedef typename bg::coordinate_type<G1>::type coordinate_type;
123 typedef typename bg::point_type<G1>::type point_type;
124
125 bg::model::multi_polygon<OutputType> out_i, out_u, out_d1, out_d2;
126
127 CalculationType area_p = p_q_area(p);
128 CalculationType area_q = p_q_area(q);
129 CalculationType area_d1 = 0, area_d2 = 0;
130
131 bg::intersection(p, q, out_i);
132 CalculationType area_i = p_q_area(out_i);
133
134 bg::union_(p, q, out_u);
135 CalculationType area_u = p_q_area(out_u);
136
137 double sum = (area_p + area_q) - area_u - area_i;
138
139 bool wrong = std::abs(sum) > settings.tolerance;
140
141 if (settings.also_difference)
142 {
143 bg::difference(p, q, out_d1);
144 bg::difference(q, p, out_d2);
145 area_d1 = p_q_area(out_d1);
146 area_d2 = p_q_area(out_d2);
147 double sum_d1 = (area_u - area_q) - area_d1;
148 double sum_d2 = (area_u - area_p) - area_d2;
149 bool wrong_d1 = std::abs(sum_d1) > settings.tolerance;
150 bool wrong_d2 = std::abs(sum_d2) > settings.tolerance;
151
152 if (wrong_d1 || wrong_d2)
153 {
154 wrong = true;
155 }
156 }
157
158 if (settings.validity)
159 {
160 std::string message;
161 if (! bg::is_valid(out_u, message))
162 {
163 std::cout << "Union is not valid: " << message << std::endl;
164 wrong = true;
165 }
166 if (! bg::is_valid(out_i, message))
167 {
168 std::cout << "Intersection is not valid: " << message << std::endl;
169 wrong = true;
170 }
171 if (settings.also_difference)
172 {
173 if (! bg::is_valid(out_d1, message))
174 {
175 std::cout << "Difference (p-q) is not valid: " << message << std::endl;
176 wrong = true;
177 }
178 if (! bg::is_valid(out_d2, message))
179 {
180 std::cout << "Difference (q-p) is not valid: " << message << std::endl;
181 wrong = true;
182 }
183 }
184
185 if (settings.verify_area && ! verify_area::apply(out_u))
186 {
187 std::cout << "Union/interior area incorrect" << std::endl;
188 wrong = true;
189 }
190 if (settings.verify_area && ! verify_area::apply(out_i))
191 {
192 std::cout << "Intersection/interior area incorrect" << std::endl;
193 wrong = true;
194 }
195 }
196
197 if (true)
198 {
199 if ((area_i > 0 && bg::touches(p, q))
200 || (area_i <= 0 && bg::intersects(p, q) && ! bg::touches(p, q)))
201 {
202 std::cout << "Wrong 'touch'! "
203 << " Intersection area: " << area_i
204 << " Touch gives: " << std::boolalpha << bg::touches(p, q)
205 << std::endl;
206 wrong = true;
207 }
208 }
209
210 bool svg = settings.svg;
211
212 if (wrong || settings.wkt)
213 {
214 if (wrong)
215 {
216 result = false;
217 svg = true;
218 }
219 bg::unique(out_i);
220 bg::unique(out_u);
221
222 std::cout
223 << "type: " << string_from_type<CalculationType>::name()
224 << " id: " << caseid
225 << " area i: " << area_i
226 << " area u: " << area_u
227 << " area p: " << area_p
228 << " area q: " << area_q
229 << " sum: " << sum;
230
231 if (settings.also_difference)
232 {
233 std::cout
234 << " area d1: " << area_d1
235 << " area d2: " << area_d2;
236 }
237 std::cout
238 << std::endl
239 << std::setprecision(9)
240 << " p: " << bg::wkt(p) << std::endl
241 << " q: " << bg::wkt(q) << std::endl
242 << " i: " << bg::wkt(out_i) << std::endl
243 << " u: " << bg::wkt(out_u) << std::endl
244 ;
245
246 }
247
248 if(svg)
249 {
250 std::ostringstream filename;
251 filename << "overlay_" << caseid << "_"
252 << string_from_type<coordinate_type>::name()
253 << string_from_type<CalculationType>::name()
254 << ".svg";
255
256 std::ofstream svg(filename.str().c_str());
257
258 bg::svg_mapper<point_type> mapper(svg, 500, 500);
259
260 mapper.add(p);
261 mapper.add(q);
262
263 // Input shapes in green/blue
264 mapper.map(p, "fill-opacity:0.5;fill:rgb(153,204,0);"
265 "stroke:rgb(153,204,0);stroke-width:3");
266 mapper.map(q, "fill-opacity:0.3;fill:rgb(51,51,153);"
267 "stroke:rgb(51,51,153);stroke-width:3");
268
269 if (settings.also_difference)
270 {
271 for (BOOST_AUTO(it, out_d1.begin()); it != out_d1.end(); ++it)
272 {
273 mapper.map(*it,
274 "opacity:0.8;fill:none;stroke:rgb(255,128,0);stroke-width:4;stroke-dasharray:1,7;stroke-linecap:round");
275 }
276 for (BOOST_AUTO(it, out_d2.begin()); it != out_d2.end(); ++it)
277 {
278 mapper.map(*it,
279 "opacity:0.8;fill:none;stroke:rgb(255,0,255);stroke-width:4;stroke-dasharray:1,7;stroke-linecap:round");
280 }
281 }
282 else
283 {
284 mapper.map(out_i, "fill-opacity:0.1;stroke-opacity:0.4;fill:rgb(255,0,128);"
285 "stroke:rgb(255,0,0);stroke-width:4");
286 mapper.map(out_u, "fill-opacity:0.1;stroke-opacity:0.4;fill:rgb(255,0,0);"
287 "stroke:rgb(255,0,255);stroke-width:4");
288 }
289 }
290 return result;
291 }
292
293 #endif // BOOST_GEOMETRY_TEST_OVERLAY_P_Q_HPP
294