boost geometry intersection give strange results

Question

I want to use the intersection function from boost geometry with a line and a polygon. I expect that the intersection is the part of the line which lies inside the polygon.

Unfortunately boost geometry returns the part of the line which lies outside of the polygon. Is this a bug in boost geometry or is something wrong with my code?

#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/ring.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/segment.hpp>
#include <boost/geometry/geometries/linestring.hpp>
#include <boost/geometry/multi/geometries/multi_point.hpp>
#include <boost/geometry/multi/geometries/multi_linestring.hpp>
#include <boost/geometry/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <boost/geometry/algorithms/intersection.hpp>

namespace bg = boost::geometry;

using value_type = double ;
using cs_type = bg::cs::cartesian;
using point_type = bg::model::point< value_type , 2 , cs_type >;
using polygon_type = bg::model::ring< point_type > ;
using line_string_type = bg::model::linestring< point_type >;
using multi_line_type = bg::model::multi_linestring< line_string_type >;

int main( int argc , char *argv[] )
{
    line_string_type line;
    line.push_back( point_type { 13.37688020921095 , 53.66231710654281 } );
    line.push_back( point_type { 13.3857295713429 , 53.6636835518369 } );
    line.push_back( point_type { 13.39213495232734 , 53.66501934623722 } );
    line.push_back( point_type { 13.39719615524716 , 53.66546436809296 } );
    line.push_back( point_type { 13.40724694386097 , 53.66240690770171 } );

    polygon_type polygon;
    polygon.push_back( point_type { 13.35 , 53.64 } );
    polygon.push_back( point_type { 13.39 , 53.64 } );
    polygon.push_back( point_type { 13.39 , 53.68 } );
    polygon.push_back( point_type { 13.35 , 53.68 } );
    polygon.push_back( point_type { 13.35 , 53.64 } );

    multi_line_type intersection;
    bg::intersection( line , polygon , intersection );

    return 0;
}

Answer 1

Make sure your input geometries satisfy the pre-conditions documented.

You can use bg::correct to fix most issues (such as proper CCW ordering of points in the polygon, closing unclosed polygons etc.):

Live On Coliru

#include <boost/geometry/core/cs.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/ring.hpp>
#include <boost/geometry/geometries/box.hpp>
#include <boost/geometry/geometries/segment.hpp>
#include <boost/geometry/geometries/linestring.hpp>
#include <boost/geometry/multi/geometries/multi_point.hpp>
#include <boost/geometry/multi/geometries/multi_linestring.hpp>
#include <boost/geometry/geometry.hpp>
#include <boost/geometry/geometries/geometries.hpp>
#include <boost/geometry/algorithms/intersection.hpp>

namespace bg = boost::geometry;

using value_type       = double;
using cs_type          = bg::cs::cartesian;
using point_type       = bg::model::point<value_type, 2, cs_type>;
using polygon_type     = bg::model::ring<point_type>;
using line_string_type = bg::model::linestring<point_type>;
using multi_line_type  = bg::model::multi_linestring<line_string_type>;

int main()
{
    line_string_type line;
    line.push_back(point_type{13.37688020921095, 53.66231710654281});
    line.push_back(point_type{13.3857295713429,  53.6636835518369});
    line.push_back(point_type{13.39213495232734, 53.66501934623722});
    line.push_back(point_type{13.39719615524716, 53.66546436809296});
    line.push_back(point_type{13.40724694386097, 53.66240690770171});
    bg::correct(line);

    polygon_type polygon;
    polygon.push_back(point_type{13.35, 53.64});
    polygon.push_back(point_type{13.39, 53.64});
    polygon.push_back(point_type{13.39, 53.68});
    polygon.push_back(point_type{13.35, 53.68});
    polygon.push_back(point_type{13.35, 53.64});
    bg::correct(polygon);

    multi_line_type intersection;
    bg::intersection(line, polygon, intersection);

    std::cout << bg::wkt(intersection);
}

This prints

MULTILINESTRING((13.3769 53.6623,13.3857 53.6637,13.39 53.6646))

With the input un-corrected it would have printed this instead:

MULTILINESTRING((13.39 53.6646,13.3921 53.665,13.3972 53.6655,13.4072 53.6624))