polygon/ring self-intersection correct

[kleunen]

I have developed a polygon correct approach (remove self-intersection, correct order, correct inners), based on boost geometry:
https://github.com/kleunen/boost_geometry_correct

The approach for calculating and removing self-intersection is vastly simpler and more efficient than existing libraries.
It was developed for tilemaker:
https://github.com/systemed/tilemaker

But it might be useful for including in boost geometry as well

[awulkiew]

Hi and thanks for this proposal!

I think that a function removing self-intersections should be a part of the library.

AFAIK @barendgehrels was working on this and the result was dissolve() algorithm:
https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/extensions/algorithms/dissolve.hpp
Have you compared the two in terms of the quality of the result (whatever this is) and performance?
If the results were different, would both of them be useful in different scenarios or would only one be needed?

In Boost.Geometry we divide algorithms into two parts, coordinate-system agnostic and specific/dependent. The former we call algorithms and the latter strategies. So strategies define some coordinate-system specific building blocks like calculating a side of a point wrt a line or performing an intersection of two segments, etc. Then the coordinate-system agnostic algorithm uses these building blocks.
Can your solution be adapted to work in non-cartesian coordinate system or is it cartesian-only?
Does your algorithm fit into this algorithm-strategy model and if yes what are the coordinate-system specific building blocks?

I'm not sure if naming it correct() is a good idea. This function is for adapting the data to the compile-time parameters of a type, i.e. point order and closure. We probably shouldn't change its behavior at this point because some of the users might unexpectedly notice a drop of performance. But let's leave the naming for later.

[kleunen]

AFAIK @barendgehrels was working on this and the result was dissolve() algorithm:
https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/extensions/algorithms/dissolve.hpp
Have you compared the two in terms of the quality of the result (whatever this is) and performance?
If the results were different, would both of them be useful in different scenarios or would only one be needed?

We tried (from the tilemaker project), but we never got this extension to compile. You can see on the github I generate the examples from @barendgehrels:
https://github.com/kleunen/boost_geometry_correct

TBH i do not understand @barendgehrels approach, so I can not judge on how this compares. He had issues on certain polygons, and this is a generalized approach. Only one approach would make sense.

There is also a discussion going on about this method on JTS:
locationtech/jts#652 (comment)

GEOS also have done an extensive comparison on polygon repair methods:
https://docs.google.com/document/d/19YEQS0goSpZlwaYivS6ZpxJ5gRth2gdG6aY2XVd5fcc/edit
Using this approach, it is possible to apply both NZW or odd-even, so you can have both behaviour regarding handling inners.

I think in your definition (correct me if I am wrong), this is an algorithm, because it uses the following coordinate specific operations:

• Calculate the intersection of two segments
• Checking if two points are equal
• Performing union and intersecting of polygons (different policies apply here).

At the moment it is a full-correct implementation, also correcting orientation, closing ring, removing invalid points. But i don't suggest integrating this, i think the self-intersecting removal is the only thing needed for boost geometry.

I think also a function like this should be part of boost geometry, many projects can benefit from this.

[barendgehrels]

Hi Wouter, thank you for approaching us!

It looks promising. I will look into it in more detail. I didn't look at our dissolve extension for quite a while, sorry if it was broken somehow. We planned to include it at some point, but now, we can maybe better select your approach which looks simpler. And if it is 100% correct - awesome!

What about the "BEER-WARE" license 😃 can it be changed into Boost Software License (your name can stay there, of course)?

[kleunen]

What about the "BEER-WARE" license 😃 can it be changed into Boost Software License (your name can stay there, of course)?

yes, of course. But you owe me a biertje, ok ? :D

It seems the topic became very relevant recently. Projects working on this and some recent scientific papers.
one thing to check if maybe any of the authors of the paper filed any patents.

[barendgehrels]

Sure, I owe you a biertje! 😃

Will you check patents?

Indeed I've troubles myself compiling dissolve too, we changed quite some strategies internally, and the extensions are sometimes not updated accordingly - there is no automatic CI for them. This makes it a bit harder to compare the algorithms.

[kleunen]

Quick search show, there are patents, related to self-intersecting polygons, but some are chinese:

https://patents.google.com/patent/US20200160569A1/en?q=self-intersecting+polygon&oq=self-intersecting+polygon+

https://patents.google.com/patent/CN112053622B/en?q=self-intersecting+polygon&oq=self-intersecting+polygon+

https://patents.google.com/patent/CN107562779B/en?q=self-intersecting+polygon&oq=self-intersecting+polygon+

https://patents.google.com/patent/US20020145616A1/en?q=self-intersecting+polygon&oq=self-intersecting+polygon+

If you can compare the results, i would be interested to hear about your findings.

[barendgehrels]

I did a comparison in terms of results and of performance. In short I think we can use this new algorithm. I like the two variants (correct and correct_odd_even). I agree with @awulkiew that the namings are not convenient, but we can fix that later.

I did not (yet) repair dissolve, but I took changelist eefd70fcb which is around 1.5 year old. That version worked for me (it's the last sha where I touched it).

This is my sample program, it's straightforward and not sophisticated

#include <boost/geometry/extensions/algorithms/dissolve.hpp>
#include <boost/geometry/extensions/multi/algorithms/dissolve.hpp>

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/geometries/multi_polygon.hpp>

#include "boost_geometry_correct/correct.hpp"
#include "boost_geometry_correct/data/CLC2006_180927.wkt.cpp"

#include <iostream>
#include <fstream>
#include <chrono>
#include <random>

template <typename Mp>
void map(std::string const& name, Mp const& mp, Mp const& result)
{
  using point = typename boost::geometry::point_type<Mp>::type;

  std::string const filename = std::string("/tmp/svg/") + name + ".svg";
  std::ofstream svg(filename.c_str());

  boost::geometry::svg_mapper<point> mapper(svg, 400, 400);

  mapper.add(mp);
  mapper.add(result);

  // Original, light yellow
  mapper.map(mp, "fill-opacity:0.3;fill:rgb(255,255,192);stroke:rgb(255,255,0);" "stroke-width:1");
  // Fixed, blue on top
  mapper.map(result, "fill-opacity:0.6;fill:rgb(0,0,255);stroke:rgb(0,0,128);" "stroke-width:2");
}

template <typename Mp>
void measure_performance(std::string const& name, Mp const& mp, char option, int count)
{
  auto const start = std::chrono::steady_clock::now();

  const double remove_spike_threshold = 1E-12;
  for (int i = 0; i < count; i++)
  {
    Mp result;
    switch(option)
    {
      case 'd' : boost::geometry::dissolve(mp, result); break;
      case 'c' : geometry::correct(mp, result, remove_spike_threshold); break;
      case 'o' : geometry::correct_odd_even(mp, result, remove_spike_threshold); break;
    }
  }
  auto const end = std::chrono::steady_clock::now();
  auto const ms = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
  std::cout << name << " " << option << " time: " << ms / 1000.0 << std::endl;
}

void compare_fix_methods(std::string const& name, std::string const& wkt, int count = 0)
{
  using point = boost::geometry::model::d2::point_xy<double>;
  using polygon = boost::geometry::model::polygon<point>;
  using multi_polygon = boost::geometry::model::multi_polygon<polygon>;

  multi_polygon mp, dissolved, corrected, oddeven;

  boost::geometry::read_wkt(wkt, mp);

  boost::geometry::dissolve(mp, dissolved);

  double const remove_spike_threshold = 1E-12;
  geometry::correct(mp, corrected, remove_spike_threshold);
  geometry::correct_odd_even(mp, oddeven, remove_spike_threshold);

  auto valid_str = [](bool v) { return v ? "valid" : "invalid"; };
  std::cout << name << " areas " << boost::geometry::area(mp)

            << " d " << boost::geometry::area(dissolved)
            << " " << valid_str(boost::geometry::is_valid(dissolved))

            << " c " << boost::geometry::area(corrected)
            << " " << valid_str(boost::geometry::is_valid(corrected))

            << " o " << boost::geometry::area(oddeven)
            << " " << valid_str(boost::geometry::is_valid(oddeven))
            << std::endl;

   map(name + "_d", mp, dissolved);
   map(name + "_c", mp, corrected);
   map(name + "_o", mp, oddeven);

   if (count > 0)
   {
     measure_performance(name, mp, 'd', count);
     measure_performance(name, mp, 'c', count);
     measure_performance(name, mp, 'o', count);
   }
   std::cout << std::endl;
}

// Generates a random polygon
std::string get_random(long long seed)
{
  std::mt19937 engine(seed);
  auto real_rand = std::bind(std::uniform_real_distribution<double>(0,1), engine);
  std::string first;
  std::string result = "MULTIPOLYGON(((";
  for (int i = 0; i < 100; i++)
  {
    std::ostringstream out;
    out << real_rand() << " " << real_rand();
    if (i > 0) { result += ","; } else { first = out.str(); }
    result += out.str();
  }
  // Close and finish it
  result += "," + first + ")))";

  return result;
}

int main()
{
  compare_fix_methods("pentagram", "MULTIPOLYGON(((5 0, 2.5 9, 9.5 3.5, 0.5 3.5, 7.5 9, 5 0)))", 10000);
  compare_fix_methods("complex", "MULTIPOLYGON(((55 10, 141 237, 249 23, 21 171, 252 169, 24 89, 266 73, 55 10)))", 10000);
  compare_fix_methods("multiple", "MULTIPOLYGON(((0 0, 10 0, 0 10, 10 10, 0 0, 5 0, 5 10, 0 10, 0 5, 10 5, 10 0, 0 0)))", 10000);
  compare_fix_methods("overlapping", "MULTIPOLYGON (((10 70, 90 70, 90 50, 30 50, 30 30, 50 30, 50 90, 70 90, 70 10, 10 10, 10 70)))", 10000);

  compare_fix_methods("random42", get_random(42), 10);
  compare_fix_methods("random123", get_random(123), 10);
  compare_fix_methods("large", wkt_CLC2006_180927, 10);

  return 0;
}

It can be compiled, for example, like this clang -O3 verify_dissolve.cpp -I ~/git_private/boost_work/ -I ~/git_other -std=c++14 -lstdc++ -lm

The results are (noted as dissolved, correct, correct o/e, first area/validity of results, then times in milliseconds):

pentagram areas 33.5 d 25.6158 valid c 25.6158 valid o 17.7317 valid
pentagram d time: 0.067
pentagram c time: 0.155
pentagram o time: 0.552

complex areas 23248.5 d 21123.3 valid c 20252.5 valid o 17256.4 valid
complex d time: 1.467
complex c time: 0.4
complex o time: 2.712

multiple areas 0 d 37.5 valid c 25 valid o 37.5 valid
multiple d time: 1.597
multiple c time: 1.137
multiple o time: 8.138

overlapping areas 4400 d 4000 valid c 4000 valid o 3600 valid
overlapping d time: 0.095
overlapping c time: 0.173
overlapping o time: 0.959

random42 areas 0.397102 d 0.10757 invalid c 0.317975 valid o 0.370735 valid
random42 d time: 0.025
random42 c time: 0.137
random42 o time: 27.539

random123 areas 0.505316 d 0.301726 invalid c 0.400009 valid o 0.362924 valid
random123 d time: 0.027
random123 c time: 0.167
random123 o time: 33.87

large areas 0.681723 d 0.681723 valid c 0.68125 invalid o 0.681723 invalid
large d time: 24.462
large c time: 8.629
large o time: 19.46

svgs (generated in /tmp/svg) give impressions of the correctness.

The performance depends on the input. This is more often the case. Sometimes dissolve is faster, sometimes correct is faster. The correct_odd_even is usually (but not always) the slowest.

For the two random generated polygons, dissolve doesn't work correctly (visually) and its result is invalid. The version of correct looks wrong too, but is marked as valid. Only the odd_even version looks correctly to me.

The large areas is using wkt_CLC2006_180927 (I changed it into a MULTI, it's the only change I made). It somehow lists the result as invalid (maybe I missed something) for correct and odd_even.

In summary: no version is just better, we can investigate more, and maybe even publish both algorithms. The correct uses repetitive calls of intersects and union (I tried to avoid that in dissolve) and that is a reason why it can be slower. Maybe it can be enhanced (didn't look at that). The results of correct_odd_even look the best and most useful to me. dissolve sometimes fails, and should be repaired. But it is often quicker (at this moment). I know that some users are using it, so we should (probably) keep it unless they are fine with alternatives too.

[barendgehrels]

One of the randomly generated polygons (yellow: input, blue: the result) of:

dissolve

correct

correct_odd_even

[barendgehrels]

I can include more pictures on request, but you can also generate them yourselves.

And it would actually be useful to get more realistic random examples (for example with generated stars on random overlapping places and fix those, I will probably adapt the random_ellipses_stars test for that)

[kleunen]

Yes, odd-even looks the most correct in this case.

This document describes also a lot of test cases which really show the differences in how the union/overlapping regions can be handled:
https://docs.google.com/document/d/19YEQS0goSpZlwaYivS6ZpxJ5gRth2gdG6aY2XVd5fcc/edit

Different approaches to union/intersect will result in different outputs. It is difficult to say what is correct and what is not.

Also, regarding the union, i now perform the "union" myself:
https://github.com/kleunen/boost_geometry_correct/blob/main/correct.hpp#L25-L47

Ideally, i would just like to use a normal union. But this does slow down the approach severely.
So the union/intersect is still something to look at.

[kleunen]

maybe it is not that strange to provide a "decompose into simple polygons" functions. And implement correct based on this, but also allow the user to just get hold of the decomposed polygon him/herself.

[kleunen]

In this case:

the non-zero winding is calculated based on the generated polygon. So you only have an outer, with an inner inside. But I guess the winding should be based on the input polygon. I am not sure how to do this easily, do i have to perform a raytracing for this, or is there some other way ?

https://github.com/kleunen/boost_geometry_correct/blob/main/correct.hpp#L290-L301

I guess maybe some sweep over the generated polygons should be performed

[kleunen]

The is also the possibility to always combine an outer as an outer. In case of the random generated polygon, this would mean, the complete outer perimeter of the polygon is traced and the inside is completely filled. This is fast, because you can just union the generated polygons. No need to detect winding/odd-even.