mardi 3 avril 2018

GEOS OverlayOp intersection operation

I am using GEOS 3.6.2 to compute an intersection between two polygons. I was able to construct my polygons, but when I try to compute the intersection it won't work, and compiling my program in Debug mode, I get the error message:

The inferior stopped because it received a signal from the operating system.

Signal name : SIGSEG

Signal meaning : Segmentation fault

Any idea where I'm wrong?

Here is my code:

#include <geos/geom/Polygon.h>
#include <geos/geom/LinearRing.h>
#include <geos/geom/CoordinateSequenceFactory.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Geometry.h>
#include <geos/operation/overlay/OverlayOp.h>
#include <iostream>
#include <array>

////////////////////////////////////////////////////////////////////////////////
geos::geom::Polygon* MakePoly(std::vector<std::vector<int>> const& polyCoords)
{
    geos::geom::GeometryFactory* factory = geos::geom::GeometryFactory::create().get();
    geos::geom::CoordinateSequence* temp = factory->getCoordinateSequenceFactory()->create((std::size_t) 0, 0);

    std::vector<std::vector<int>>::const_iterator it_x = polyCoords.begin();
    int size = it_x->size();

    for (int i=0; i<size; i++)
    {
        temp->add(geos::geom::Coordinate(polyCoords[0][i], polyCoords[1][i]));
    }

    geos::geom::LinearRing *shell=factory->createLinearRing(temp);


    //NULL in this case could instead be a collection of one or more holes
    //in the interior of the polygon
    return factory->createPolygon(shell,NULL);
}


////////////////////////////////////////////////////////////////////////////////
int main()
{
    // Create geometry.
    std::vector<std::vector<int>> polyCoords1 = {
            {1, 1, 2, 2, 1, 1, 4, 5, 4, 1},
            {1, 2, 2, 4, 4, 5, 5, 3, 1, 1}
    };

    geos::geom::Polygon* poly1 = MakePoly(polyCoords1);

    std::vector<std::vector<int>> polyCoords2 = {
            {4, 4, 6, 6, 4},
            {1, 5, 5, 1, 1}
    };

    geos::geom::Polygon* poly2 = MakePoly(polyCoords2);

    // Actually perform the operation.
    geos::operation::overlay::OverlayOp intersection(poly1, poly2);

    // Extracting the geometry of the intersection (position of the error).
    geos::geom::Geometry* intersectionGeo = intersection.getResultGeometry( geos::operation::overlay::OverlayOp::OpCode::opINTERSECTION );

    std::cout<<intersectionGeo->getArea()<<std::endl;
}

Aucun commentaire:

Enregistrer un commentaire