Boost::Geometry intersection not working as expected

Question

I have a little test application to intersect a few rectangles with boost::geometry.

typedef boost::geometry::model::point
<
double, 2, boost::geometry::cs::cartesian
> point;

typedef boost::geometry::model::polygon<point > polygon;

polygon Intersect(polygon p1, polygon p2) {
    std::vector < polygon > result;
    boost::geometry::intersection(p1, p2, result);
    return result.front();
}

polygon IntersectionTest() {
    polygon one, two, three, four;

    boost::geometry::read_wkt("POLYGON((35 25, 35 35, 15 35, 15 25, 35 25))", one);
    boost::geometry::read_wkt("POLYGON((45 30, 45 50, 25 50, 25 30, 45 30))", two);
    boost::geometry::read_wkt("POLYGON((50  0, 50 40, 10 40, 10  0, 50  0))", three);
    boost::geometry::read_wkt("POLYGON((40 20, 40 60,  0 60,  0 20, 40 20))", four);

    return Intersect(Intersect(Intersect(one, two), three), four);
}

I must do something wrong, because I expect the result to be something like (35 30, 35 40, 25 40, 25 30, 35 30) yet I get 10 points long polygons containing points like 50 0 while intersections of parallel rectangles should always be rectangles with 4+1 points and 50 0 is at the edge so it shouldn't be in the intersection at all. If I put it into an SVG, the rectangles seems to be as I expect them.

What can be wrong? If it is a bug in boost::geometry how can I make sure it is? (I'm currently using 1.48.) If it is a bug is there a way to circumvent the problem?

score 3 · Accepted Answer · answered Nov 28 '11 at 10:35

The polygon class has the following template parameters:

template
<
    typename Point,
    bool ClockWise = true,
    bool Closed = true,
    template<typename, typename> class PointList = std::vector,
    template<typename, typename> class RingList = std::vector,
    template<typename> class PointAlloc = std::allocator,
    template<typename> class RingAlloc = std::allocator
>
class polygon {...}

As you see the second template parameter defaults to true. That means that the points that define the polygon are considered to be / should be in clockwise order.

This is your problem actually.

If you look at your WKT form of the geometries you wrote the points in counterclockwise order.

So you should either:

typedef boost::geometry::model::polygon<point,false> polygon;

Or write the points in the WKT string in clockwise order. I.e:

boost::geometry::read_wkt("POLYGON((35 25, 15 25, 15 35, 35 35, 35 25))", one);
boost::geometry::read_wkt("POLYGON((45 30, 25 30, 25 50, 45 50, 45 30))", two);

With this, the result will be as you expected it.

As a personal comment it would be probably be a nice thing if read_wkt would enforce the correct direction on reading...

isn't 35 25, 15 25, 15 35, 35 35, 35 25 counter clockwise? – Bedir Yilmaz Aug 18 '23 at 08:12 — Bedir Yilmaz, Aug 18 '23 at 08:12

score 2 · Answer 2 · answered Jan 19 '12 at 17:19

I agree with the answer of ds27680.

A third option is to call boost::geometry::correct(a geometry) after read_wkt, and users are encouraged to do this, if not sure about the orientation.

read_wkt does indeed not enforce the correct orientation. The reason is that the correct function requires an area computation, which can be saved if users know the polygon is in the right order (usually it is known).

Boost::Geometry intersection not working as expected

2 Answers2