Geometry Boost List: Re: [geometry] Multipolygon from vector of rings

From: Tinko Bartels (tinkobartels_at_[hidden])
Date: 2021-04-04 10:13:58

Hi,

I'm not entirely sure I understand the question correctly, so I'll try to
cover some cases of what could be meant, starting with the most trivial just
to rule out misunderstandings.

If your geometries are interpreted as rings in the Boost.Geometry sense
(filled features) and they are nonintersecting in the DE-9IM-sense, which, I
think, is disjoint, i.e. the areas they cover, including boundaries, do not
share a single point, then all your rings are exterior rings.

If that nonintersection property is only meant to say that they are all
closed curves around some areas whose boundaries do not intersect but whose
interiors maybe do, then I understand the problem like this: You need to
know which rings are contained in which other rings. E.g., if we have 3
rings A, B, C where C is contained in B and B is contained in A, then we
have two polygons one with A as exterior ring and B as interior ring, and a
second polygon with C as exterior ring and no interior rings. If all rings
are clockwise and your coordinates are cartesian, then, using the means
available in Boost.Geometry, my first idea would be something like this
(warning: this probably leaves a lot of room for optimization):

#include <numeric>
#include <vector>
#include <algorithm>
#include <iostream>
#include <fstream>

#include <boost/geometry.hpp>

namespace bg = boost::geometry;

using point = bg::model::d2::point_xy<double>;
using box = bg::model::box<point>;
using ring = bg::model::ring<point>;
using polygon = bg::model::polygon<point>;
using multi_polygon = bg::model::multi_polygon<polygon>;

multi_polygon make_multi_polygon(std::vector<ring> const& rings)
{
    std::vector<box> envelopes;
    envelopes.reserve(rings.size());
    for(const auto& r : rings)
        envelopes.push_back(bg::return_envelope<box>(r));
    std::vector<std::size_t> indices(rings.size());
    std::iota(indices.begin(), indices.end(), 0);
    std::sort(indices.begin(), indices.end(),
              [&envelopes](std::size_t a, std::size_t b)
                {
                    return bg::area(envelopes[a]) > bg::area(envelopes[b]);
                });
    std::vector<bool> interior(rings.size(), false);
    std::vector<std::size_t> parent(rings.size(), -1);
    for(int i = 0; i < rings.size(); ++i)
    {
        auto ri = indices[i];
        auto env = envelopes[ri];
        for(int j = i - 1; j >= 0; --j)
        {
            auto rj = indices[j];
            if( !bg::within(env, envelopes[rj]) ) continue;
            if( bg::within(rings[ri], rings[rj] ) )
            {
                if(interior[j]) break;
                interior[i] = true;
                parent[i] = j;
            }
        }
    }

    multi_polygon result;
    for(int i = 0; i < rings.size(); ++i)
    {
        if(interior[i]) continue;
        auto ri = indices[i];
        result.push_back(polygon());
        auto& poly_i = result.back();
        poly_i.outer() = rings[ri];
        for(int j = i + 1; j < rings.size(); ++j)
            if(parent[j] == i)
            {
                poly_i.inners().push_back(rings[indices[j]]);
                bg::reverse(poly_i.inners().back());
            }
    }
    return result;
}

int main()
{
    std::vector<ring> input{
        {{-12,-12}, {-12,42},{12,42},{12,-12},{-12,-12}},
        {{-10,-10}, {-10,10},{10,10}, {10,-10},{-10,-10}},
        {{-5,-5}, {-5,5},{5,5}, {5,-5},{-5,-5}},
        {{-1,-1}, {-1,1},{1,1}, {1,-1},{-1,-1}},
        {{-.1,-.1}, {-.1,.1},{.1,.1}, {.1,-.1},{-.1,-.1}},
        {{-10,20}, {-10,40},{10,40}, {10,20},{-10,20}},
        {{-8,22}, {-8,38}, {-4,38}, {-4,22}, {-8,22}},
        {{2,22}, {2,38}, {8,38}, {8,22}, {2,22}}
    };
    auto merged = make_multi_polygon(input);
    std::ofstream svg("my_map.svg");
    boost::geometry::svg_mapper<point> mapper(svg, 800, 800);
    mapper.add(merged);
    mapper.map(merged,
"fill-opacity:0.3;fill:rgb(51,51,153);stroke:rgb(51,51,153);stroke-width:1");
    return 0;
}

I don't know if that's sufficiently efficient for your usecase. I guess the
box within-tests make it at least quadratic in the number of rings in the
average case (which could probably be solved by using an rtree). If the
within-tests for the rings dominate the runtime, which I'd expect if the
rings are large compared to the number of rings, then I have no idea right
now how to improve that using algorithms in Boost.Geometry. Maybe a
sweep-line based approach could be faster. If it's non-cartesian, then I'm
not sure how the envelops in BG work there, but I think the general approach
should still be mostly applicable.

--
Sent from: http://boost-geometry.203548.n3.nabble.com/