Union of arbitrary number of overlapping polygons.

[awulkiew]

@barendgehrels is there an algorithm allowing to calculate union of any number of potentially overlapping (multi-)polygons? In buffer maybe? If there is does it work for all coordinate systems?

[barendgehrels]

No, we don't have that out of the box. All our algorithms work on either one or two.
Users can iterate (might be slow) and or recursively combine them, similar (or reusing) what we recently did for the alternative for dissolve (to be continued)

[awulkiew]

I was thinking about buffer because a union of multiple overlapping rings of buffer parts is calculated there. The purpose is not to calculate valid multipolygon from an invalid one (like in the alternative for dissolve) but to calculate a union of some number of polygons like AFAIU it's done in buffer for parts.

[barendgehrels]

True, but it’s done piecewise there. It handles internal overlaps in the same way as exrernal overlaps. Therefore it’s different

[barendgehrels]

I have some code for it, will send it this weekend

[barendgehrels]

This works well:

template <typename Geometry, typename It>
Geometry recursive_union(It begin, It end, std::size_t max_distance) {
  Geometry unioned = {};
  if (begin == end) {
    return unioned;
  }

  const auto dist = std::distance(begin, end);
  if (dist >= max_distance) {
    It mid = begin;
    std::advance(mid, dist / 2);

    const Geometry left = recursive_union<Geometry>(begin, mid, max_distance);
    const Geometry right = recursive_union<Geometry>(mid, end, max_distance);

    bg::union_(left, right, unioned);

    return unioned;
  }

  // Iteratively add this small amount of polygons
  for (It it = begin; it != end; ++it) {
    multipol_t combined;
    bg::union_(unioned, *it, combined);
    unioned = std::move(combined);
  }
  return unioned;
}

and you can call it by for example

vector<multipol_t> vector;
// fill it somehow
// it works better if you sort them too (for example with `geometry::less`)
// perform the union
auto bigunion = recursive_union<multipol_t>(vector.begin(), vector.end(), 16);

[awulkiew]

Interesting, it should indeed be faster than the iterative version. Have you compared the performance of both? If not then I'll do it. I'd probably sort polygons spatially based on their bounding boxes. AFAIR geometry::less works only for points which could be centroids of the boxes or polygons but geometry::less is not "proper" space partitioning. I also have another idea for an algorithm using rtree. So I'll see which one is the fastest and when. Thanks for the algorithm!

[barendgehrels]

Yes, the performance of the recursive union is very very much better then unioning iteratively.

I added this for another project (therefore I first had to ask permission, but it's there, we can reuse it).

Figures:
Recursive (sorted before on center points)

80277 polygons
Elapsed 2725[ms]

Iteratively

80277 polygons
Elapsed 232184[ms]

[barendgehrels]

Two other points:

• Apart from that it is faster, it's also better. Somewhere in the iterative union there are some geometry errors. I will report that later. With recursive union, that doesn't occur (or less frequently)
• Also, it could be used internally within the library. Because if we do a multi-point buffer, over a multipoint of 80000 points, it's terribly slow. With the recursive approach that will be much faster (didn't try that yet, but it will be about the difference of the figures above).

[barendgehrels]

By the way, we could sort on MortonCode or GeoHash too. That would be a cool small addition to the library. It works well.
https://en.wikipedia.org/wiki/Z-order_curve
https://en.wikipedia.org/wiki/Geohash

[awulkiew]

Based on your suggestion I implemented an algorithm avoiding recursive function calls:
https://github.com/boostorg/geometry/blob/develop/include/boost/geometry/algorithms/merge_elements.hpp#L165

[mendes-jose]

The following implementation seems to be quite fast too. Faster than the recursive one in my case (spatially ordered polygons)
It's similar to it, as it uses the dichotomy principle too. But it avoids recursion all together.
It consumes the input vector by storing intermediate unions in it though ⚠️
The complete union will be at to_unify.front() after calling this.

template <typename MutiPolygonType>
void cascadeUnion(std::vector<MutiPolygonType>& to_unify)
{
    size_t step = 1;
    size_t half_step;

    // the outer loop doubles the distance between two
    // polygons to be merged at every iteration
    do
    {
        half_step = step;
        step *= 2;
        size_t i = 0;

        // the inner loop merges polygons at i and i+half_step storing the result at i
        do
        {
            MutiPolygonType unified;
            bg::union_(to_unify.at(i), to_unify.at(i + half_step), unified);

            to_unify.at(i) = std::move(unified);
        
            i += step;
        
        } while (i + half_step < to_unify.size());
    
    } while (step < to_unify.size());
}

[awulkiew]

Thanks for sharing it. It'd be interesting to compare the performance with what I implemented.

[viveksood97]

How could we make this work for a std::vector<boost::geometry::polygon Polygon>?

[mendes-jose]

The example implementation uses template, so it should work with a vector of boost::geometry::model::polygon<> as well as a vector of boost::geometry::model::multi_polygon<>.

A more generic version that would work with other types of containers (e.g. std::list) is given below:

namespace bg = boost::geometry;

namespace boost::geometry
{
    using Point = model::d2::point_xy<double>;
    using Polygon = model::polygon<Point>;
    using MPolygon = model::multi_polygon<Polygon>;
}

template<class It>
void cascadeUnion(It begin, It end)
{
    if (begin == end)
        return;

    size_t step = 1;
    size_t half_step;

    // this outer loop doubles the distance between two
    // polygons to be merged at every iteration
    do
    {
        half_step = step;
        step *= 2;
        size_t i = 0;

        // this inner loop merges polygons at i and i+half_step storing the result at i
        do
        {
            bg::MPolygon unified;
            bg::union_(*(begin + i),
                       *(begin + i + half_step),
                       unified);

            if (unified.empty())
                bg::clear(*(begin + i));
            else
                *(begin + i) = std::move(unified);

            i += step;

        } while (i + half_step < std::distance(begin, end));

    } while (step < std::distance(begin, end));
}

Usage:

std::vector<bg::Polygon> to_merge; // could be a std::list<bg::MPolygon> object or something else

// fill up the to_merge vector with spatially ordered polygons

cascadeUnion(to_merge.begin(), to_merge.end());

// merged polygon is in to_merge.front()