Search code examples
c++boostboost-geometry

How to return intersecting linestrings in boost::geometry::intersection(ring1, ring2, vector_of_linestring)?


I have 2 rings A and B, and I want to use boost::geometry::intersection() to return linestrings (the orange arrow ones):

2 rings intersection screenshot

But my code only returns the intersecting points P1 and P2. Which part should I modify?

#include <boost/geometry.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/register/ring.hpp>
#include <boost/geometry/geometries/register/point.hpp>
#include <iostream>

namespace bg = boost::geometry;
typedef bg::model::point<double, 2, bg::cs::cartesian> point_t;
typedef bg::model::linestring<point_t> linestring_t;
typedef bg::model::ring<point_t> ring_t;
typedef bg::model::polygon<point_t> polygon_t;

typedef bg::model::multi_point<point_t> mpoint_t;
typedef bg::model::multi_linestring<linestring_t> mlinestring_t;
typedef bg::model::multi_polygon<polygon_t> mpolygon_t;

int main() {
    point_t ptA0(0, 0);
    point_t ptA1(10, 0);
    point_t ptA2(10, 10);
    point_t ptA3(0, 10);
    ring_t ringA;
    bg::append(ringA, ptA0);
    bg::append(ringA, ptA1);
    bg::append(ringA, ptA2);
    bg::append(ringA, ptA3);
    bg::append(ringA, ptA0);
    
    point_t ptB0(5, -5);
    point_t ptB1(15, -5);
    point_t ptB2(15, 5);
    point_t ptB3(5, 5);
    ring_t ringB;
    bg::append(ringB, ptB0);
    bg::append(ringB, ptB1);
    bg::append(ringB, ptB2);
    bg::append(ringB, ptB3);
    bg::append(ringB, ptB0);  

    std::vector<linestring_t> resline;      
    bg::intersection(ringB, ringA, resline);
    for (int i = 0; i < resline.size(); i++) {
        std::cout << bg::dsv(resline[i]) << std::endl;
    }   
    return 0;
}

Below is the output: (10, 5) is P1 and (5, 0) is P2, which is not what I expected.

((10, 5), (10, 5)) 
((5, 0), (5, 0))

Solution

  • I've tried to get a handle on this problem using just the DE-9IM that Boost Geometry implements: https://godbolt.org/z/KWzvzExr7 which outputs https://pastebin.ubuntu.com/p/9ck6gcPK5P/

    ---- void do_test(Input, Input) [with Input = boost::geometry::model::ring<boost::geometry::model::point<int, 2, boost::geometry::cs::cartesian> >]
    a: POLYGON((0 0,0 10,10 10,10 0,0 0))
    b: POLYGON((5 -5,5 5,15 5,15 -5,5 -5))
    MULTIPOINT((10 5),(5 0))
    MULTIPOLYGON(((10 5,10 0,5 0,5 5,10 5)))
    Segment of ringB LINESTRING(5 -5,5 5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(5 5,15 5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(15 5,15 -5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(15 -5,5 -5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
     ---- void do_test(Input, Input) [with Input = boost::geometry::model::linestring<boost::geometry::model::point<int, 2, boost::geometry::cs::cartesian> >]
    a: LINESTRING(0 0,0 10,10 10,10 0,0 0)
    b: LINESTRING(5 -5,5 5,15 5,15 -5,5 -5)
    MULTIPOINT((10 5),(5 0))
    Segment of ringB LINESTRING(5 -5,5 5) relates to ringA 
    ['F', 'F', '1']
    ['F', 'F', 'F']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(5 5,15 5) relates to ringA 
    ['F', 'F', '1']
    ['F', 'F', 'F']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(15 5,15 -5) relates to ringA 
    ['F', 'F', '1']
    ['F', 'F', 'F']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(15 -5,5 -5) relates to ringA 
    ['F', 'F', '1']
    ['F', 'F', 'F']
    ['F', 'F', '2']
     ---- void do_test(Input, Input) [with Input = boost::geometry::model::polygon<boost::geometry::model::point<int, 2, boost::geometry::cs::cartesian> >]
    a: POLYGON((0 0,0 10,10 10,10 0,0 0))
    b: POLYGON((5 -5,5 5,15 5,15 -5,5 -5))
    MULTIPOINT((10 5),(5 0))
    MULTIPOLYGON(((10 5,10 0,5 0,5 5,10 5)))
    Segment of ringB LINESTRING(5 -5,5 5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(5 5,15 5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(15 5,15 -5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
    Segment of ringB LINESTRING(15 -5,5 -5) relates to ringA 
    ['F', 'F', '2']
    ['F', 'F', '1']
    ['F', 'F', '2']
    

    but as far as I can see all the results are at least fully disjoint. So you might have to zoom in on the points instead of the segments:

    Live On Wandbox

    #include <boost/geometry.hpp>
    #include <boost/geometry/geometries/ring.hpp>
    #include <boost/geometry/iterators/segment_iterator.hpp>
    #include <iostream>
    
    namespace bg = boost::geometry;
    using point_t  = bg::model::point<int, 2, bg::cs::cartesian>;
    using ring_t   = bg::model::ring<point_t>;
    using mpoint_t = bg::model::multi_point<point_t>;
    
    void do_test(ring_t ringA, ring_t ringB)
    {
        std::string reason;
        if (!bg::is_valid(ringA, reason)) { std::cout << "warning ringA: " << reason << "\n"; bg::correct(ringA); }
        if (!bg::is_valid(ringB, reason)) { std::cout << "warning ringB: " << reason << "\n"; bg::correct(ringB); }
    
        std::cout << "ringA: " << bg::wkt(ringA) << std::endl;
        std::cout << "ringB: " << bg::wkt(ringB) << std::endl;
    
        mpoint_t result;
        bg::intersection(ringA, ringB, result);
        std::cout << bg::wkt(result) << std::endl;
    
        for (auto seg : boost::make_iterator_range(bg::segments_begin(ringB),
                                                   bg::segments_end(ringB))) {
            for (auto& p : result)
                if (bg::intersects(p, seg))
                    std::cout << bg::wkt(p) << " intersects " << bg::wkt(seg) << "\n";
        }
    }
    
    int main() {
        do_test({{0, 0}, {0, 10}, {10, 10}, {10, 0}/*, {0, 0}*/},
                {{5, -5}, {5, 5}, {15, 5}, {15, -5}/*, {5, -5}*/});
    }
    

    Prints

    warning ringA: Geometry is defined as closed but is open
    warning ringB: Geometry is defined as closed but is open
    ringA: POLYGON((0 0,0 10,10 10,10 0,0 0))
    ringB: POLYGON((5 -5,5 5,15 5,15 -5,5 -5))
    MULTIPOINT((10 5),(5 0))
    POINT(5 0) intersects LINESTRING(5 -5,5 5)
    POINT(10 5) intersects LINESTRING(5 5,15 5)