Difference with rectilinear multipolygon with integer coordinates produces invalid polygon with disconnected interior

[tinko92]

The code for reproduction is

#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>

namespace bg = boost::geometry;
using point = bg::model::point<int, 2, bg::cs::cartesian>;
using polygon = bg::model::polygon<point>;
using multipoly = bg::model::multi_polygon<polygon>;

int main() {
    auto p1 = bg::from_wkt<multipoly>("MULTIPOLYGON(((1 2,1 3,0 3,0 4,1 4,1 5,3 5,3 0,0 0,0 2,1 2),(1 2,1 1,2 1,2 2,1 2)))");
    auto p2 = bg::from_wkt<polygon>("POLYGON((1 3,1 4,2 4,2 3,1 3))");
    multipoly p3;
    bg::difference(p1, p2, p3);
    std::cout << bg::is_valid(p1) << " , " << bg::is_valid(p2) << "\n";
    std::cout << bg::is_valid(p3) << " , " << bg::wkt(p3) << "\n";
    return 0;
}

The output is

1 , 1
0 , MULTIPOLYGON(((1 4,1 5,3 5,3 0,0 0,0 2,1 2,1 3,0 3,0 4,1 4),(1 4,1 3,2 3,2 4,1 4),(1 2,1 1,2 1,2 2,1 2)))

which covers the right area but is invalid, but it should be two polygons (looks like the first counter example in Figure 12 of the section on polygons in the OGC spec, page 28). If the hole from p1 is removed, a valid output is produced. If it is read as a polygon rather than as a multipolygon, the problem also disappears. It shouldn't be a numerical issue (I think) because the coordinates are integral.

This is what the input looks like (p1 is semi-transparent blue, 2 is semi-transparent red):

And the result in green:

[tinko92]

This issue is fixed on #1350 .

[barendgehrels]

This issue is fixed on #1350 .

Confirmed. Failed before.