1 // Boost.Geometry (aka GGL, Generic Geometry Library)
2 // Robustness Test
3
4 // Copyright (c) 2012-2015 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 #if defined(_MSC_VER)
11 # pragma warning( disable : 4244 )
12 # pragma warning( disable : 4267 )
13 #endif
14
15 #include <fstream>
16 #include <sstream>
17
18 #include <boost/foreach.hpp>
19 #include <boost/program_options.hpp>
20 #include <boost/random/linear_congruential.hpp>
21 #include <boost/random/uniform_int.hpp>
22 #include <boost/random/uniform_real.hpp>
23 #include <boost/random/variate_generator.hpp>
24 #include <boost/timer.hpp>
25
26 #include <boost/geometry.hpp>
27 #include <boost/geometry/geometries/geometries.hpp>
28 #include <boost/geometry/geometries/point_xy.hpp>
29
30 #include <boost/geometry/io/svg/svg_mapper.hpp>
31 #include <boost/geometry/extensions/algorithms/midpoints.hpp>
32
33 #include <common/common_settings.hpp>
34 #include <common/make_square_polygon.hpp>
35
36
37 namespace bg = boost::geometry;
38
39 template <typename Geometry1, typename Geometry2, typename Geometry3>
create_svg(std::string const & filename,Geometry1 const & mp,Geometry2 const & ls,Geometry3 const & difference,Geometry3 const & intersection)40 void create_svg(std::string const& filename
41 , Geometry1 const& mp
42 , Geometry2 const& ls
43 , Geometry3 const& difference
44 , Geometry3 const& intersection
45 )
46 {
47 typedef typename boost::geometry::point_type<Geometry1>::type point_type;
48
49
50 std::ofstream svg(filename.c_str());
51 boost::geometry::svg_mapper<point_type> mapper(svg, 800, 800);
52
53 boost::geometry::model::box<point_type> box;
54 bg::envelope(mp, box);
55 bg::buffer(box, box, 1.0);
56 mapper.add(box);
57 bg::envelope(ls, box);
58 bg::buffer(box, box, 1.0);
59 mapper.add(box);
60
61 mapper.map(mp, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:3");
62 mapper.map(ls, "stroke-opacity:0.9;stroke:rgb(0,0,0);stroke-width:1");
63
64 mapper.map(intersection,"opacity:0.6;stroke:rgb(0,128,0);stroke-width:5");
65 mapper.map(difference, "opacity:0.6;stroke:rgb(255,0,255);stroke-width:5"); //;stroke-dasharray:1,7;stroke-linecap:round
66 }
67
68
69 template <typename Linestring, typename Generator, typename Settings>
make_random_linestring(Linestring & line,Generator & generator,Settings const & settings)70 inline void make_random_linestring(Linestring& line, Generator& generator, Settings const& settings)
71 {
72 using namespace boost::geometry;
73
74 typedef typename point_type<Linestring>::type point_type;
75 typedef typename coordinate_type<Linestring>::type coordinate_type;
76
77 coordinate_type x, y;
78 x = generator();
79 y = generator();
80
81 int count = 3 + generator() % 6;
82 int d = 0; // direction
83
84 for (int i = 0; i <= count && x <= settings.field_size; i++, x++, d = 1 - d)
85 {
86 append(line, make<point_type>(x, y + d));
87 append(line, make<point_type>(x, y + 1 - d));
88 }
89
90 if (d == 0 && generator() % 4 < 3 && y >= 2)
91 {
92 d = 1 - d;
93 x--;
94 y -= 2;
95 count = 3 + generator() % 6;
96 for (int i = 0; i <= count && x >= 0; i++, x--, d = 1 - d)
97 {
98 append(line, make<point_type>(x, y + d));
99 append(line, make<point_type>(x, y + 1 - d));
100 }
101 }
102
103 //if (settings.triangular)
104 //{
105 // // Remove a point, generator says which
106 // int c = generator() % 4;
107 // if (c >= 1 && c <= 3)
108 // {
109 // ring.erase(ring.begin() + c);
110 // }
111 //}
112 }
113
114 template<typename Geometry>
115 class inside_check
116 {
117 Geometry const& m_geo;
118 bool& m_result;
119 public :
120
inside_check(Geometry const & geo,bool & result)121 inside_check(Geometry const& geo, bool& result)
122 : m_geo(geo)
123 , m_result(result)
124 {}
125
operator =(inside_check<Geometry> const & input)126 inline inside_check<Geometry> operator=(inside_check<Geometry> const& input)
127 {
128 return inside_check<Geometry>(input.m_geo, input.m_result);
129 }
130
131 template <typename Point>
operator ()(Point const & p)132 inline void operator()(Point const& p)
133 {
134 if (! bg::covered_by(p, m_geo))
135 {
136 if (m_result)
137 {
138 std::cout << "Does not fulfill inside check" << std::endl;
139 }
140 m_result = false;
141 }
142 }
143 };
144
145
146 template<typename Geometry>
147 class outside_check
148 {
149 Geometry const& m_geo;
150 bool& m_result;
151 public :
outside_check(Geometry const & geo,bool & result)152 outside_check(Geometry const& geo, bool& result)
153 : m_geo(geo)
154 , m_result(result)
155 {}
156
operator =(outside_check<Geometry> const & input)157 inline outside_check<Geometry> operator=(outside_check<Geometry> const& input)
158 {
159 return outside_check<Geometry>(input.m_geo, input.m_result);
160 }
161
162 template <typename Point>
operator ()(Point const & p)163 inline void operator()(Point const& p)
164 {
165 if (bg::within(p, m_geo))
166 {
167 if (m_result)
168 {
169 std::cout << "Does not fulfill outside check" << std::endl;
170 }
171 m_result = false;
172 }
173 }
174 };
175
176 template<typename Segment>
177 class border2_check
178 {
179 Segment const& m_segment;
180 bool& m_result;
181
182 public :
border2_check(Segment const & seg,bool & result)183 border2_check(Segment const& seg, bool& result)
184 : m_segment(seg)
185 , m_result(result)
186 {}
187
operator =(border2_check<Segment> const & input)188 inline border2_check<Segment> operator=(border2_check<Segment> const& input)
189 {
190 return border2_check<Segment>(input.m_segment, input.m_result);
191 }
192
193 template <typename Segment2>
operator ()(Segment2 const & segment)194 inline void operator()(Segment2 const& segment)
195 {
196 // Create copies (TODO: find out why referring_segment does not compile)
197 typedef typename bg::point_type<Segment2>::type pt;
198 typedef bg::model::segment<pt> segment_type;
199
200 typedef bg::strategy::intersection::relate_cartesian_segments
201 <
202 bg::policies::relate::segments_intersection_points
203 <
204 segment_type,
205 segment_type,
206 bg::segment_intersection_points<pt>
207 >
208 > policy;
209
210 segment_type seg1, seg2;
211 bg::convert(m_segment, seg1);
212 bg::convert(segment, seg2);
213 bg::segment_intersection_points<pt> is = policy::apply(seg1, seg2);
214
215 if (is.count == 2)
216 {
217 if (m_result)
218 {
219 std::cout << "Does not fulfill border2 check" << std::endl;
220 }
221 m_result = true;
222 }
223 }
224 };
225
226 template<typename Geometry>
227 class border_check
228 {
229 Geometry const& m_geo;
230 bool& m_result;
231 public :
border_check(Geometry const & geo,bool & result)232 border_check(Geometry const& geo, bool& result)
233 : m_geo(geo)
234 , m_result(result)
235 {}
236
operator =(border_check<Geometry> const & input)237 inline border_check<Geometry> operator=(border_check<Geometry> const& input)
238 {
239 return border_check<Geometry>(input.m_geo, input.m_result);
240 }
241
242 template <typename Segment>
operator ()(Segment const & s)243 inline void operator()(Segment const& s)
244 {
245 bool on_border = false;
246 border2_check<Segment> checker(s, on_border);
247 bg::for_each_segment(m_geo, checker);
248
249 if (on_border)
250 {
251 if (m_result)
252 {
253 std::cout << "Does not fulfill border check" << std::endl;
254 }
255 m_result = false;
256 }
257 }
258 };
259
260
261 template <typename MultiPolygon, typename Linestring, typename Settings>
verify(std::string const & caseid,MultiPolygon const & mp,Linestring const & ls,Settings const & settings)262 bool verify(std::string const& caseid, MultiPolygon const& mp, Linestring const& ls, Settings const& settings)
263 {
264 bg::model::multi_linestring<Linestring> difference, intersection;
265 bg::difference(ls, mp, difference);
266 bg::intersection(ls, mp, intersection);
267
268 //typedef typename bg::length_result_type<point_type>::type length_type;
269
270 bool result = true;
271
272 // 1) Check by length
273 typedef double length_type;
274 length_type len_input = bg::length(ls);
275 length_type len_difference = bg::length(difference);
276 length_type len_intersection = bg::length(intersection);
277 if (! bg::math::equals(len_input, len_difference + len_intersection))
278 {
279 std::cout << "Input: " << len_input
280 << " difference: " << len_difference
281 << " intersection: " << len_intersection
282 << std::endl;
283
284 std::cout << "Does not fulfill length check" << std::endl;
285
286 result = false;
287 }
288
289 // 2) Check by within and covered by
290 inside_check<MultiPolygon> ic(mp, result);
291 bg::for_each_point(intersection, ic);
292
293 outside_check<MultiPolygon> oc(mp, result);
294 bg::for_each_point(difference, oc);
295
296 border_check<MultiPolygon> bc(mp, result);
297 bg::for_each_segment(difference, bc);
298
299 // 3) check also the mid-points from the difference to remove false positives
300 BOOST_FOREACH(Linestring const& d, difference)
301 {
302 Linestring difference_midpoints;
303 bg::midpoints(d, false, std::back_inserter(difference_midpoints));
304 outside_check<MultiPolygon> ocm(mp, result);
305 bg::for_each_point(difference_midpoints, ocm);
306 }
307
308
309 bool svg = settings.svg;
310 bool wkt = settings.wkt;
311 if (! result)
312 {
313 std::cout << "ERROR " << caseid << std::endl;
314 std::cout << bg::wkt(mp) << std::endl;
315 std::cout << bg::wkt(ls) << std::endl;
316 svg = true;
317 wkt = true;
318 }
319
320 if (svg)
321 {
322 std::ostringstream filename;
323 filename << caseid << "_"
324 << typeid(typename bg::coordinate_type<Linestring>::type).name()
325 << ".svg";
326 create_svg(filename.str(), mp, ls, difference, intersection);
327 }
328
329 if (wkt)
330 {
331 std::ostringstream filename;
332 filename << caseid << "_"
333 << typeid(typename bg::coordinate_type<Linestring>::type).name()
334 << ".wkt";
335 std::ofstream stream(filename.str().c_str());
336 stream << bg::wkt(mp) << std::endl;
337 stream << bg::wkt(ls) << std::endl;
338 }
339
340 return result;
341 }
342
343 template <typename MultiPolygon, typename Generator>
test_linear_areal(MultiPolygon & result,int & index,Generator & generator,int level,common_settings const & settings)344 bool test_linear_areal(MultiPolygon& result, int& index,
345 Generator& generator,
346 int level, common_settings const& settings)
347 {
348 MultiPolygon p, q;
349
350 // Generate two boxes
351 if (level == 0)
352 {
353 p.resize(1);
354 q.resize(1);
355 make_square_polygon(p.front(), generator, settings);
356 make_square_polygon(q.front(), generator, settings);
357 bg::correct(p);
358 bg::correct(q);
359 }
360 else
361 {
362 bg::correct(p);
363 bg::correct(q);
364 if (! test_linear_areal(p, index, generator, level - 1, settings)
365 || ! test_linear_areal(q, index, generator, level - 1, settings))
366 {
367 return false;
368 }
369 }
370
371 typedef typename boost::range_value<MultiPolygon>::type polygon;
372
373 MultiPolygon mp;
374 bg::detail::union_::union_insert
375 <
376 polygon
377 >(p, q, std::back_inserter(mp));
378
379 bg::unique(mp);
380 bg::simplify(mp, result, 0.01);
381 bg::correct(mp);
382
383 // Generate a linestring
384 typedef typename bg::point_type<MultiPolygon>::type point_type;
385 typedef bg::model::linestring<point_type> linestring_type;
386 linestring_type ls;
387 make_random_linestring(ls, generator, settings);
388
389 std::ostringstream out;
390 out << "recursive_la_" << index++ << "_" << level;
391 return verify(out.str(), mp, ls, settings);
392 }
393
394
395 template <typename T, bool Clockwise, bool Closed>
test_all(int seed,int count,int level,common_settings const & settings)396 void test_all(int seed, int count, int level, common_settings const& settings)
397 {
398 boost::timer t;
399
400 typedef boost::minstd_rand base_generator_type;
401
402 base_generator_type generator(seed);
403
404 boost::uniform_int<> random_coordinate(0, settings.field_size - 1);
405 boost::variate_generator<base_generator_type&, boost::uniform_int<> >
406 coordinate_generator(generator, random_coordinate);
407
408 typedef bg::model::polygon
409 <
410 bg::model::d2::point_xy<T>, Clockwise, Closed
411 > polygon;
412 typedef bg::model::multi_polygon<polygon> mp;
413
414
415 int index = 0;
416 for(int i = 0; i < count; i++)
417 {
418 mp p;
419 test_linear_areal<mp>(p, index, coordinate_generator, level, settings);
420 }
421 std::cout
422 << "geometries: " << index
423 << " type: " << typeid(T).name()
424 << " time: " << t.elapsed() << std::endl;
425 }
426
main(int argc,char ** argv)427 int main(int argc, char** argv)
428 {
429 try
430 {
431 namespace po = boost::program_options;
432 po::options_description description("=== recursive_polygons_linear_areal ===\nAllowed options");
433
434 int count = 1;
435 int seed = static_cast<unsigned int>(std::time(0));
436 int level = 3;
437 bool ccw = false;
438 bool open = false;
439 common_settings settings;
440 std::string form = "box";
441
442 description.add_options()
443 ("help", "Help message")
444 ("seed", po::value<int>(&seed), "Initialization seed for random generator")
445 ("count", po::value<int>(&count)->default_value(1), "Number of tests")
446 ("diff", po::value<bool>(&settings.also_difference)->default_value(false), "Include testing on difference")
447 ("level", po::value<int>(&level)->default_value(3), "Level to reach (higher->slower)")
448 ("form", po::value<std::string>(&form)->default_value("box"), "Form of the polygons (box, triangle)")
449 ("ccw", po::value<bool>(&ccw)->default_value(false), "Counter clockwise polygons")
450 ("open", po::value<bool>(&open)->default_value(false), "Open polygons")
451 ("size", po::value<int>(&settings.field_size)->default_value(10), "Size of the field")
452 ("wkt", po::value<bool>(&settings.wkt)->default_value(false), "Create a WKT of the inputs, for all tests")
453 ("svg", po::value<bool>(&settings.svg)->default_value(false), "Create a SVG for all tests")
454 ;
455
456 po::variables_map varmap;
457 po::store(po::parse_command_line(argc, argv, description), varmap);
458 po::notify(varmap);
459
460 if (varmap.count("help")
461 || (form != "box" && form != "triangle"))
462 {
463 std::cout << description << std::endl;
464 return 1;
465 }
466
467 settings.triangular = form != "box";
468
469 if (ccw && open)
470 {
471 test_all<double, false, false>(seed, count, level, settings);
472 }
473 else if (ccw)
474 {
475 test_all<double, false, true>(seed, count, level, settings);
476 }
477 else if (open)
478 {
479 test_all<double, true, false>(seed, count, level, settings);
480 }
481 else
482 {
483 test_all<double, true, true>(seed, count, level, settings);
484 }
485
486 #if defined(HAVE_TTMATH)
487 // test_all<ttmath_big, true, true>(seed, count, max, svg, level);
488 #endif
489 }
490 catch(std::exception const& e)
491 {
492 std::cout << "Exception " << e.what() << std::endl;
493 }
494 catch(...)
495 {
496 std::cout << "Other exception" << std::endl;
497 }
498
499 return 0;
500 }
501