Boost::geometry::intersection with C++

Question

I'm a PhD student in civil engineering and I recently started doing some coding in C++, basically I'm interested in getting the overlapping or intersection area of two polygons, which represent the projection of two soil particles.

I did a lot of search and I found that boost geometry is the best solution for me. I also did a lot of search for the specific problem I'm facing but I wasn't able to resolve my issue.

Here is the problem, the software I'm using is called PFC3D (particle flow code). I have to use microsoft visual studio 2010 to interact with this software and compile DLL files to run it in PFC.

My code work very well without the overlapping area. Here is the code:

// Includes for overlapping
#include <boost/geometry.hpp>
#include <boost/geometry/core/point_type.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/register/point.hpp>enter code here
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>
typedef boost::geometry::model::polygon<boost::geometry::model::d2::point_xy<double> > polygon;
polygon poly1, poly2;
poly1 {{0.0,  0.0}, {0.0, 1.0}, {1.0, 1.0}, {1.0,  0.0}, {0.05,  0.0}};
poly2 {{0.5, -0.5}, {0.5, 0.5}, {1.5, 0.5}, {1.5, -0.5},  {0.5, -0.5}};
std::deque<polygon> output;
boost::geometry::intersection(poly1, poly2, output);
double area = boost::geometry::area(output);

The error I'm getting is in assigning the poly1 and poly2 coordinates. Hope you can help in this. Thanks!

Answer 1

points in the polygon should be sorted by clockwise. such as

#include <iostream>
#include <algorithm>
#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/algorithms/area.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>

namespace bg = boost::geometry;
namespace bgm = bg::model;
typedef bgm::d2::point_xy<double> point;
typedef bgm::polygon<point> polygon;

void testOverlap_anticlockwise();
void testOverlap_clockwise();


int main() {
    testOverlap_anticlockwise();
    testOverlap_clockwise();
}



void testOverlap_anticlockwise()
{
    polygon poly1;
    poly1.outer().push_back(point(15, 0.176766));
    poly1.outer().push_back(point(0.176767, 15));
    poly1.outer().push_back(point(-5, 16.6844));
    poly1.outer().push_back(point(-5, -5));
    poly1.outer().push_back(point(16.6844, -5));
    poly1.outer().push_back(point(15, 0.176766));

    polygon poly2;
    poly2.outer().push_back(point(15, 15));
    poly2.outer().push_back(point(0, 15));
    poly2.outer().push_back(point(0, 0));
    poly2.outer().push_back(point(15, 0));
    poly2.outer().push_back(point(15, 15));

    std::cout << bg::wkt(poly1) << "\n";
    std::cout << bg::wkt(poly2) << "\n";
    std::deque<polygon> output;
    bg::intersection(poly1, poly2, output);

    double area = 0;
    for (std::deque<polygon>::const_iterator it = output.begin(); it != output.end(); ++it) {
        area += bg::area(*it);
    }
    std::cout << area << std::endl;
}

void testOverlap_clockwise()
{
    polygon poly1;
    poly1.outer().push_back(point(15, 0.176766));
    poly1.outer().push_back(point(0.176767, 15));
    poly1.outer().push_back(point(-5, 16.6844));
    poly1.outer().push_back(point(-5, -5));
    poly1.outer().push_back(point(16.6844, -5));
    poly1.outer().push_back(point(15, 0.176766));
    std::reverse(poly1.outer().begin(), poly1.outer().end());

    polygon poly2;
    poly2.outer().push_back(point(15, 15));
    poly2.outer().push_back(point(0, 15));
    poly2.outer().push_back(point(0, 0));
    poly2.outer().push_back(point(15, 0));
    poly2.outer().push_back(point(15, 15));
    std::reverse(poly2.outer().begin(), poly2.outer().end());

    std::cout << bg::wkt(poly1) << "\n";
    std::cout << bg::wkt(poly2) << "\n";
    std::deque<polygon> output;
    bg::intersection(poly1, poly2, output);

    double area = 0;
    for (std::deque<polygon>::const_iterator it = output.begin(); it != output.end(); ++it) {
        area += bg::area(*it);
    }
    std::cout << area << std::endl;
}

result is POLYGON((15 0.176766,0.176767 15,-5 16.6844,-5 -5,16.6844 -5,15 0.176766)) POLYGON((15 15,0 15,0 0,15 0,15 15)) 0

POLYGON((15 0.176766,16.6844 -5,-5 -5,-5 16.6844,0.176767 15,15 0.176766)) POLYGON((15 15,15 0,0 0,0 15,15 15)) 115.136

Answer 2

Well. identifier { } only ever works if identifier is a typename.

If you wanted uniform initialization, you use { } to commence the constructor parameter list, and wrap each parameter ring in an extra set of { }:

polygon poly1 { { { 0.0, 0.0 }, { 0.0, 1.0 }, { 1.0, 1.0 }, { 1.0, 0.0 }, { 0.05, 0.0 } }  };
polygon poly2 { { { 0.5, -0.5 }, { 0.5, 0.5 }, { 1.5, 0.5 }, { 1.5, -0.5 }, { 0.5, -0.5 } }  };

Next, area doesn't expect a multi-polygon, so write a loop:

double area = 0;
for (auto& p : output)
    area += boost::geometry::area(p);

May I suggest looking at dsv parsing for the input:

polygon poly1, poly2;
bg::read<bg::format_wkt>(poly1, "POLYGON((0 0,0 1,1 1,1 0,0.05 0,0 0))");
bg::read<bg::format_wkt>(poly2, "POLYGON((0.5 -0.5,0.5 0.5,1.5 0.5,1.5 -0.5,0.5 -0.5))");

Live demo: Live On Coliru

#include <iostream>
#include <boost/geometry.hpp>
#include <boost/geometry/io/io.hpp>
#include <boost/geometry/algorithms/area.hpp>
#include <boost/geometry/geometries/point.hpp>
#include <boost/geometry/geometries/point_xy.hpp>
#include <boost/geometry/geometries/polygon.hpp>

namespace bg = boost::geometry;
namespace bgm = bg::model;
typedef bgm::polygon<bgm::d2::point_xy<double> > polygon;

int main() {
    polygon poly1, poly2;
    bg::read<bg::format_wkt>(poly1, "POLYGON((0 0,0 1,1 1,1 0,0.05 0,0 0))");
    bg::read<bg::format_wkt>(poly2, "POLYGON((0.5 -0.5,0.5 0.5,1.5 0.5,1.5 -0.5,0.5 -0.5))");

    std::cout << bg::wkt(poly1) << "\n";
    std::cout << bg::wkt(poly2) << "\n";
    std::deque<polygon> output;
    bg::intersection(poly1, poly2, output);


    double area = 0;
    for (auto& p : output)
        area += bg::area(p);
}