boost::polygon boolean subtract results in additional lines

Question

I have a rectangle with a circle that is located inside the rectangle. Now I do a boolean subtract operation where the circle is subtracted from the larger rectangle using boost::polygon and the following (pseudo) code:

boost::polygon::polygon_with_holes_data<int> rect;
boost::polygon::polygon_with_holes_data<int> circle;
// fill both with data here
boost::polygon::polygon_set_data<int> result=rect-circle;

After that "result" contains a bunch of coordinate-pairs where each of them are marked with an additional value -1 or 1. Then I fetch the resulting data the following way:

std::vector<boost::polygon::polygon_data<int>> outputData;
result.get<std::vector<boost::polygon::polygon_data<int>>>(outputData);

Now "outputData" contains a list of coordinates which form the following result:

So the hole is incorporated within one complete, closed polygon coordinate array where I would have expected something like boost::polygon::polygon_with_holes_data where the outline and the holes would have been separated.

So my question: am I doing something wrong here? How do I get the outline and the hole separated so that I can get rid of this extra line which connects them? Is there a function to fetch boost::polygon::polygon_with_holes_data out of the returned boost::polygon::polygon_set_data or what else is the magic here?

Thanks :-)

Can you include the source data? It may be invalid. E.g. not closed, wrong orientation etc. If your input doesn't satisfy the pre-conditions, the behaviour of the algorithms will be undefined (GIGO) — sehe
@sehe: let's say I subtract a triangle {{-78020,-72350},{110650,-7660},{7090,69510},{-78020,-72350}} from a rectangle {{-147530,127670},{163140,127670},{163140,-143280},{-147530,-143280},{-147530,127670}} - hopefully there is no typo in my coordinates. Both are closed and both result n one single array of 11 output coordinate elements which again contain this additional line from outline to hole. — Elmi

sehe sehe · Accepted Answer · 2020-07-01T00:19:03

I took the trouble to generate some data myself. I don't know what data you use, nor how you render it. My suspicion is something is wrong in those processes:

Let's create a polygon_data factories:

template <typename T>
auto make_rect(T x1, T y1, T x2, T y2) {
    auto r = bp::rectangle_data(x1, y1, x2, y2);
    bp::polygon_data<T> p;
    bp::assign(p, r);
    return p;
}

So far so good. For the circle I used your new discovery from your other question:

template <typename T>
auto make_circle(T x, T y, T radius) {
    std::vector<bp::polygon_data<T> > ps;

    bp::assign(ps, bp::rectangle_data(x-1, y-1, x+1, y+1));
    bp::resize(ps, radius, true, 512);

    assert(ps.size() == 1);
    return ps.front();
}

Now let's start out with the straight-forward:

int main() {
    auto rect       = make_rect(0, 0, 1000, 1000);
    auto circle     = make_circle(500, 500, 400);
    auto difference = rect - circle;

    std::cout << std::fixed;
    std::cout << "rect area: " << area(rect)   << "\n";
    std::cout << "circle area: " << area(circle) << "\n";
    std::cout << "difference area: " << area(difference) << "\n";

Prints

rect area: 1000000.000000
circle area: 505481.000000
difference area: 494519.000000

Convert And Check

Now to see whether any of these geometries would render correctly, I know no other technique than via Boost Geometry¹

So, let's convert to a polygon_with_holes_data and do some checks:

bp::polygon_with_holes_data<int> result(rect.begin(), rect.end(),
        &circle, &circle+1);

std::cout << "result area: " << area(result) << "\n";

assert(area(rect) - area(circle) == area(difference));
assert(bp::equivalence(result, rect - circle));

The asserts pass, and this prints the expected:

result area: 494519.000000

Which correctly matches the difference area.

Rendering SVG

Now we adapt in Boost Geometry and write a SVG:

namespace bg = boost::geometry;
std::cout << bg::wkt(rect);
std::cout << bg::wkt(result);

{
    using C = bg::coordinate_type<decltype(result)>::type;
    using PT = bg::model::d2::point_xy<C>;
    std::ofstream svg("output.svg");
    boost::geometry::svg_mapper<PT> mapper(svg, 400, 400);
    mapper.add(result);
    mapper.map(result, "fill-opacity:0.3;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
}

This writes the WKT and output.svg:

Looks fine here.

What Could Be Wrong? How To Debug?

Like I mentioned in my comments, many library/tools have restrictions/invariants on the input they take. Things like the orientation of outer vs inner rings, self intersections and the presence of closing points on polygon rings.

Now, if you don't know what data is acceptable, you can hit the docs, or you can use the library API to diagnose. For example in Boost Geometry I like to make a quick-and-dirty diagnostic function:

template<typename G>
void auto_correct(G& geo) {
    std::string reason;
    while (!bg::is_valid(geo, reason)) {
        std::cout << "Trying to auto-correct broken input data: " << reason << "\n";
        bg::correct(geo);
    }
}

You can add wkt printing to debug what changes are actually applied.

See e.g. boost read_wkt produced invalid polygon for examples of diagnostics and fixes.

(Note: don't use this in production, the corrections may not do what you expect and the loop may even be infinite for some inputs)

Live Demo

As always: Live On Coliru

#include <boost/polygon/polygon.hpp>
#include <boost/polygon/gtl.hpp>
#include <boost/polygon/rectangle_data.hpp>
#include <boost/polygon/polygon_set_data.hpp>
#include <boost/polygon/polygon_data.hpp>
#include <iostream>

namespace bp = boost::polygon;

#ifndef NO_GEO
#include <fstream>
#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/geometries/adapted/boost_polygon.hpp>
namespace bg = boost::geometry;

/* 
 * template<typename G>
 * void auto_correct(G& geo) {
 *     std::string reason;
 *     while (!bg::is_valid(geo, reason)) {
 *         std::cout << "Trying to auto-correct broken input data: " << reason << "\n";
 *         bg::correct(geo);
 *     }
 * }
 */
#endif

template <typename T>
auto make_rect(T x1, T y1, T x2, T y2) {
    auto r = bp::rectangle_data(x1, y1, x2, y2);
    bp::polygon_data<T> p;
    bp::assign(p, r);
    return p;
}

template <typename T>
auto make_circle(T x, T y, T radius) {
    std::vector<bp::polygon_data<T> > ps;

    bp::assign(ps, bp::rectangle_data(x-1, y-1, x+1, y+1));
    bp::resize(ps, radius, true, 512);

    assert(ps.size() == 1);
    return ps.front();
}

int main() {
    auto rect       = make_rect(0, 0, 1000, 1000);
    auto circle     = make_circle(500, 500, 400);
    auto difference = rect - circle;

    std::cout << std::fixed;
    std::cout << "rect   area: " << area(rect)   << "\n";
    std::cout << "circle area: " << area(circle) << "\n";
    std::cout << "difference area: " << area(difference) << "\n";

    bp::polygon_with_holes_data<int> result(rect.begin(), rect.end(),
            &circle, &circle+1);

    std::cout << "result area: " << area(result) << "\n";

    assert(area(rect) - area(circle) == area(difference));
    assert(bp::equivalence(result, rect - circle));

#ifndef NO_GEO
    std::cout << bg::wkt(rect) << "\n";
    std::cout << bg::wkt(result) << "\n";

    {
        using C = bg::coordinate_type<decltype(result)>::type;
        using PT = bg::model::d2::point_xy<C>;
        std::ofstream svg("output.svg");
        boost::geometry::svg_mapper<PT> mapper(svg, 400, 400);
        mapper.add(result);
        mapper.map(result, "fill-opacity:0.3;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
    }
#endif
}