libs/geometry/test/robustness/overlay/linear_areal/recursive_polygons_linear_areal.cpp

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