Geometry Boost List: [ggl] Buggy polygon difference case

Subject: [ggl] Buggy polygon difference case
From: Pal Tengerecki (tengereckipal)
Date: 2011-09-20 04:55:42

Hi,

I isolated a case of two polygons where the output of boost::geometry::difference seems to be incorect. We compute the difference of
"MULTIPOLYGON(((18000 18000,18000 -18000,-18000 -18000,-18000 18000,18000 18000),(18000 18000,15928 1521,16000 0,18000 18000)))"
and
"MULTIPOLYGON(((15711 3028,18000 18000,15928 1521,15711 3028)))".
The answer produced by boost is
"MULTIPOLYGON(((15928 1521,16000 0,18000 18000,15711 3028,15928 1521)))"

which is definitely wrong. The first polygon has a single outer and a single inner loop, touching the outer loop in a vertex. The second polygon is a triangle touching the inner loop if the first one in an edge, as well as the outer loop of the first one in a vertex.

I do not think that precision of the number type has a role, as the input has only integers, and no nontrivial intersections can be computed from the input data. The other thing I am not sure in, is the correctness of the representation of the first polygon. It could be given also by a single self-intersecting outer loop. However, this was extracted from the result of another boost geometry boolean operation, so it should entirely satisfy the corresponding boost concepts.

Is that a bug, or how should I tweak my code to result a correct answer?

Thanks,

??Pal

Full example code:

#include <iostream>
#include <fstream>
#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
#include <boost/geometry/multi/geometries/multi_polygon.hpp>
#include <boost/geometry/domains/gis/io/wkt/wkt.hpp>
#include <boost/geometry/extensions/io/svg/svg_mapper.hpp>

template <typename Geometry1>
void create_svg(std::string const& filename, Geometry1 const& a) {
??typedef typename boost::geometry::point_type<Geometry1>::type point_type;
??std::ofstream svg(filename.c_str());
??boost::geometry::svg_mapper<point_type> mapper(svg, 400, 400);
??mapper.add(a);?
??mapper.map(a, "fill-opacity:0.5;fill:rgb(153,204,0);stroke:rgb(153,204,0);stroke-width:2");
}

int main() {
??typedef boost::geometry::model::d2::point_xy<double> point;
??typedef boost::geometry::model::polygon<point> polygon;
??typedef boost::geometry::model::multi_polygon<polygon> multi_polygon;

??multi_polygon mp1;
??std::string s1 = "MULTIPOLYGON(((18000 18000,18000 -18000,-18000 -18000,-18000 18000,18000 18000),(18000 18000,15928 1521,16000 0,18000 18000)))";
??boost::geometry::read_wkt(s1, mp1);
??create_svg("mp1.svg", mp1);
??boost::geometry::correct(mp1);
??std::cout << boost::geometry::wkt(mp1) << std::endl;

??multi_polygon mp2;
??std::string s2 = "MULTIPOLYGON(((15711 3028,18000 18000,15928 1521,15711 3028)))";
??boost::geometry::read_wkt(s2, mp2);
??create_svg("mp2.svg", mp2);
??boost::geometry::correct(mp2);
??std::cout << boost::geometry::wkt(mp2) << std::endl;

??multi_polygon output;
??boost::geometry::difference(mp1, mp2, output);
??create_svg("output.svg", output);
??std::cout << boost::geometry::wkt(output) << std::endl;
}
-------------- next part --------------
An HTML attachment was scrubbed...
URL: http://lists.osgeo.org/pipermail/ggl/attachments/20110919/959f20cf/attachment.html