diff --git a/include/geos/operation/overlayng/OverlayNG.h b/include/geos/operation/overlayng/OverlayNG.h index 9a5c23a5dd..7fed36bba9 100644 --- a/include/geos/operation/overlayng/OverlayNG.h +++ b/include/geos/operation/overlayng/OverlayNG.h @@ -57,9 +57,16 @@ namespace overlayng { // geos.operation.overlayng * * DIFFERENCE - all points which lie in the first geometry but not the second * * SYMDIFFERENCE - all points which lie in one geometry but not both * - * Input geometries may have different dimension. - * Input collections must be homogeneous - * (all elements must have the same dimension). + * The requirements for overlay input are: + * * Input collections must be homogeneous + * (all elements must have the same dimension). + * * Inputs may be simple link GeometryCollections. + * A GeometryCollection is simple if it can be flattened into a valid Multi-geometry; + * i.e. it is homogeneous and does not contain any overlapping Polygons. + * * In general, inputs must be valid geometries. + * However, polygonal inputs may contain the following two kinds of "mild" invalid topology: + * (i) rings which self-touch at discrete points (sometimes called inverted shells and exverted holes). + * (ii) rings which touch along line segments (i.e. topology collapse). * * The precision model used for the computation can be supplied * independent of the precision model of the input geometry. @@ -82,7 +89,7 @@ namespace overlayng { // geos.operation.overlayng * since the intersection clipping optimization can * interact with the snapping to alter the result. * - * TOptionally the overlay computation can process using strict mode + * Optionally the overlay computation can process using strict mode * (via setStrictMode(boolean). In strict mode result semantics are: * * - Lines and Points resulting from topology collapses are not included @@ -101,7 +108,6 @@ namespace overlayng { // geos.operation.overlayng * * The original JTS overlay semantics correspond to non-strict mode. * - * * If a robustness error occurs, a TopologyException is thrown. * These are usually caused by numerical rounding causing the noding * output to not be fully noded. @@ -403,4 +409,3 @@ class GEOS_DLL OverlayNG { } // namespace geos.operation.overlayng } // namespace geos.operation } // namespace geos - diff --git a/include/geos/operation/overlayng/OverlayPoints.h b/include/geos/operation/overlayng/OverlayPoints.h index 87597c813a..2e9c153038 100644 --- a/include/geos/operation/overlayng/OverlayPoints.h +++ b/include/geos/operation/overlayng/OverlayPoints.h @@ -82,10 +82,6 @@ class GEOS_DLL OverlayPoints { std::map> buildPointMap(const Geometry* geom); - Coordinate roundCoord(const Point* pt, const PrecisionModel* pm) const; - - - public: /** @@ -120,4 +116,3 @@ class GEOS_DLL OverlayPoints { } // namespace geos.operation.overlayng } // namespace geos.operation } // namespace geos - diff --git a/src/algorithm/locate/IndexedPointInAreaLocator.cpp b/src/algorithm/locate/IndexedPointInAreaLocator.cpp index da233e20ee..0182903efe 100644 --- a/src/algorithm/locate/IndexedPointInAreaLocator.cpp +++ b/src/algorithm/locate/IndexedPointInAreaLocator.cpp @@ -52,11 +52,19 @@ IndexedPointInAreaLocator::IntervalIndexedGeometry::init(const geom::Geometry& g // pre-compute size of segment vector std::size_t nsegs = 0; for(const geom::LineString* line : lines) { + //-- only include rings of Polygons or LinearRings + if (! line->isClosed()) + continue; + nsegs += line->getCoordinatesRO()->size() - 1; } index = decltype(index)(10, nsegs); for(const geom::LineString* line : lines) { + //-- only include rings of Polygons or LinearRings + if (! line->isClosed()) + continue; + addLine(line->getCoordinatesRO()); } } @@ -90,12 +98,6 @@ IndexedPointInAreaLocator::buildIndex(const geom::Geometry& g) IndexedPointInAreaLocator::IndexedPointInAreaLocator(const geom::Geometry& g) : areaGeom(g) { - const std::type_info& areaGeomId = typeid(areaGeom); - if(areaGeomId != typeid(geom::Polygon) - && areaGeomId != typeid(geom::MultiPolygon) - && areaGeomId != typeid(geom::LinearRing)) { - throw util::IllegalArgumentException("Argument must be Polygonal or LinearRing"); - } } geom::Location diff --git a/src/operation/overlayng/OverlayMixedPoints.cpp b/src/operation/overlayng/OverlayMixedPoints.cpp index de2067d5b7..f79087f8bc 100644 --- a/src/operation/overlayng/OverlayMixedPoints.cpp +++ b/src/operation/overlayng/OverlayMixedPoints.cpp @@ -29,15 +29,48 @@ #include #include - - namespace geos { // geos namespace operation { // geos.operation namespace overlayng { // geos.operation.overlayng - using namespace geos::geom; +/** + * @brief Extracts and rounds coordinates from a geometry + * + */ +class CoordinateExtractingFilter: public geom::CoordinateFilter { +public: + CoordinateExtractingFilter(CoordinateArraySequence& p_pts, const PrecisionModel& p_pm) + : pts(p_pts), pm(p_pm) + {} + + /** + * Destructor. + * Virtual dctor promises appropriate behaviour when someone will + * delete a derived-class object via a base-class pointer. + * http://www.parashift.com/c++-faq-lite/virtual-functions.html#faq-20.7 + */ + ~CoordinateExtractingFilter() override {} + + /** + * Performs a filtering operation with or on coord in "read-only" mode. + * @param coord The "read-only" Coordinate to which + * the filter is applied. + */ + void + filter_ro(const geom::Coordinate* coord) override + { + Coordinate p(*coord); + pm.makePrecise(p); + pts.add(p); + } + +private: + CoordinateArraySequence& pts; + const PrecisionModel& pm; +}; + /*public*/ OverlayMixedPoints::OverlayMixedPoints(int p_opCode, const Geometry* geom0, const Geometry* geom1, const PrecisionModel* p_pm) : opCode(p_opCode) @@ -234,16 +267,9 @@ std::unique_ptr OverlayMixedPoints::extractCoordinates(const Geometry* points, const PrecisionModel* p_pm) const { std::unique_ptr coords(new CoordinateArraySequence()); - std::size_t n = points->getNumGeometries(); - for (std::size_t i = 0; i < n; i++) { - const Point* point = static_cast(points->getGeometryN(i)); - if (point->isEmpty()) { - continue; - } - Coordinate coord; - OverlayUtil::round(point, p_pm, coord); - coords->add(coord, true); - } + + CoordinateExtractingFilter filter(*coords, *p_pm); + points->apply_ro(&filter); return coords; } @@ -275,12 +301,6 @@ OverlayMixedPoints::extractLines(const Geometry* geom) const return list; } - - - - - - } // namespace geos.operation.overlayng } // namespace geos.operation } // namespace geos diff --git a/src/operation/overlayng/OverlayPoints.cpp b/src/operation/overlayng/OverlayPoints.cpp index 8c1c1bf91a..c328687b56 100644 --- a/src/operation/overlayng/OverlayPoints.cpp +++ b/src/operation/overlayng/OverlayPoints.cpp @@ -27,6 +27,47 @@ namespace geos { // geos namespace operation { // geos.operation namespace overlayng { // geos.operation.overlayng +struct PointExtractingFilter: public GeometryComponentFilter { + + PointExtractingFilter(std::map>& p_ptMap, const PrecisionModel* p_pm) + : ptMap(p_ptMap), pm(p_pm) + {} + + void + filter_ro(const Geometry* geom) + { + if (geom->getGeometryTypeId() != GEOS_POINT) return; + + const Point* pt = static_cast(geom); + // don't add empty points + if (pt->isEmpty()) return; + + Coordinate p = roundCoord(pt, pm); + /** + * Only add first occurrence of a point. + * This provides the merging semantics of overlay + */ + if (ptMap.find(p) == ptMap.end()) { + std::unique_ptr newPt(pt->getFactory()->createPoint(p)); + ptMap[p] = std::move(newPt); + } + } + + static Coordinate + roundCoord(const Point* pt, const PrecisionModel* p_pm) + { + const Coordinate* p = pt->getCoordinate(); + if (OverlayUtil::isFloating(p_pm)) + return *p; + Coordinate p2 = *p; + p_pm->makePrecise(p2); + return p2; + } + +private: + std::map>& ptMap; + const PrecisionModel* pm; +}; /*public static*/ std::unique_ptr @@ -126,54 +167,16 @@ OverlayPoints::computeUnion(std::map>& map0, } } - /*private*/ std::map> OverlayPoints::buildPointMap(const Geometry* geom) { std::map> map; - for (std::size_t i = 0; i < geom->getNumGeometries(); i++) { - const Geometry* elt = geom->getGeometryN(i); - if (elt->getGeometryTypeId() != GEOS_POINT) { - throw util::IllegalArgumentException("Non-point geometry input to point overlay"); - } - // don't add empty points - if (elt->isEmpty()) continue; - - const Point* pt = static_cast(elt); - Coordinate p = roundCoord(pt, pm); - /** - * Only add first occurrence of a point. - * This provides the merging semantics of overlay - */ - if (map.find(p) == map.end()) { - std::unique_ptr newPt(pt->getFactory()->createPoint(p)); - map[p] = std::move(newPt); - } - } + PointExtractingFilter filter(map, pm); + geom->apply_ro(&filter); return map; } - -/*private*/ -Coordinate -OverlayPoints::roundCoord(const Point* pt, const PrecisionModel* p_pm) const -{ - const Coordinate* p = pt->getCoordinate(); - if (OverlayUtil::isFloating(p_pm)) - return *p; - Coordinate p2 = *p; - p_pm->makePrecise(p2); - return p2; -} - - - - - - - - } // namespace geos.operation.overlayng } // namespace geos.operation } // namespace geos diff --git a/tests/unit/operation/overlayng/OverlayNGGeometryCollectionTest.cpp b/tests/unit/operation/overlayng/OverlayNGGeometryCollectionTest.cpp new file mode 100644 index 0000000000..55d32d5154 --- /dev/null +++ b/tests/unit/operation/overlayng/OverlayNGGeometryCollectionTest.cpp @@ -0,0 +1,155 @@ +// +// Test Suite for geos::operation::overlayng::OverlayNG class for GeometryCollections. + +#include +#include + +// geos +#include + +// std +#include + +using namespace geos::geom; +using namespace geos::operation::overlayng; +using geos::io::WKTReader; +using geos::io::WKTWriter; + +namespace tut { +// +// Test Group +// + +// Common data used by all tests +struct test_overlaynggc_data { + + WKTReader r; + WKTWriter w; + + void + testOverlay(const std::string& a, const std::string& b, const std::string& expected, int opCode) + { + PrecisionModel pm; + std::unique_ptr geom_a = r.read(a); + std::unique_ptr geom_b = r.read(b); + std::unique_ptr geom_expected = r.read(expected); + std::unique_ptr geom_result = OverlayNG::overlay(geom_a.get(), geom_b.get(), opCode, &pm); + // std::string wkt_result = w.write(geom_result.get()); + // std::cout << std::endl << wkt_result << std::endl; + ensure_equals_geometry(geom_expected.get(), geom_result.get()); + } + + void + testIntersection(const std::string& a, const std::string& b, const std::string& expected) + { + testOverlay(a, b, expected, OverlayNG::INTERSECTION); + } + + void + testUnion(const std::string& a, const std::string& b, const std::string& expected) + { + testOverlay(a, b, expected, OverlayNG::UNION); + } +}; + +typedef test_group group; +typedef group::object object; + +group test_overlaynggc_group("geos::operation::overlayng::OverlayNGGeometryCollection"); + +// +// Test Cases +// + +// testSimpleA_mP +template<> +template<> +void object::test<1> () +{ + std::string a = "POLYGON ((0 0, 0 1, 1 1, 0 0))"; + std::string b = "GEOMETRYCOLLECTION ( MULTIPOINT ((0 0), (99 99)) )"; + testIntersection(a, b, + "POINT (0 0)"); + testUnion(a, b, + "GEOMETRYCOLLECTION (POINT (99 99), POLYGON ((0 0, 0 1, 1 1, 0 0)))"); +} + +// testSimpleP_mP +template<> +template<> +void object::test<2> () +{ + std::string a = "POINT(0 0)"; + std::string b = "GEOMETRYCOLLECTION ( MULTIPOINT ((0 0), (99 99)) )"; + testIntersection(a, b, + "POINT (0 0)"); + testUnion(a, b, + "MULTIPOINT ((0 0), (99 99))"); +} + +// testSimpleP_mL +template<> +template<> +void object::test<3> () +{ + std::string a = "POINT(5 5)"; + std::string b = "GEOMETRYCOLLECTION ( MULTILINESTRING ((1 9, 9 1), (1 1, 9 9)) )"; + testIntersection(a, b, + "POINT (5 5)"); + testUnion(a, b, + "MULTILINESTRING ((1 1, 5 5), (1 9, 5 5), (5 5, 9 1), (5 5, 9 9))"); +} + +// testSimpleP_mA +template<> +template<> +void object::test<4> () +{ + std::string a = "POINT(5 5)"; + std::string b = "GEOMETRYCOLLECTION ( MULTIPOLYGON (((1 1, 1 5, 5 5, 5 1, 1 1)), ((9 9, 9 5, 5 5, 5 9, 9 9))) )"; + testIntersection(a, b, + "POINT (5 5)"); + testUnion(a, b, + "MULTIPOLYGON (((1 1, 1 5, 5 5, 5 1, 1 1)), ((9 9, 9 5, 5 5, 5 9, 9 9)))"); +} + +// testSimpleP_AA +template<> +template<> +void object::test<5> () +{ + std::string a = "POINT(5 5)"; + std::string b = "GEOMETRYCOLLECTION ( POLYGON ((1 1, 1 5, 5 5, 5 1, 1 1)), POLYGON ((9 9, 9 5, 5 5, 5 9, 9 9)) )"; + testIntersection(a, b, + "POINT (5 5)"); + testUnion(a, b, + "MULTIPOLYGON (((1 1, 1 5, 5 5, 5 1, 1 1)), ((9 9, 9 5, 5 5, 5 9, 9 9)))"); + } + +// testSimpleL_AA +template<> +template<> +void object::test<6> () +{ + std::string a = "LINESTRING (0 0, 10 10)"; + std::string b = "GEOMETRYCOLLECTION ( POLYGON ((1 1, 1 5, 5 5, 5 1, 1 1)), POLYGON ((9 9, 9 5, 5 5, 5 9, 9 9)) )"; + testIntersection(a, b, + "MULTILINESTRING ((1 1, 5 5), (5 5, 9 9))"); + testUnion(a, b, + "GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1), LINESTRING (9 9, 10 10), POLYGON ((1 1, 1 5, 5 5, 5 1, 1 1)), POLYGON ((5 5, 5 9, 9 9, 9 5, 5 5)))"); + } + +// testSimpleA_AA +template<> +template<> +void object::test<7> () +{ + std::string a = "POLYGON ((2 8, 8 8, 8 2, 2 2, 2 8))"; + std::string b = "GEOMETRYCOLLECTION ( POLYGON ((1 1, 1 5, 5 5, 5 1, 1 1)), POLYGON ((9 9, 9 5, 5 5, 5 9, 9 9)) )"; + testIntersection(a, b, + "MULTIPOLYGON (((2 2, 2 5, 5 5, 5 2, 2 2)), ((5 5, 5 8, 8 8, 8 5, 5 5)))"); + testUnion(a, b, + "POLYGON ((1 1, 1 5, 2 5, 2 8, 5 8, 5 9, 9 9, 9 5, 8 5, 8 2, 5 2, 5 1, 1 1))"); +} + +} // namespace tut diff --git a/tests/xmltester/tests/general/TestNGOverlayGC.xml b/tests/xmltester/tests/general/TestNGOverlayGC.xml new file mode 100644 index 0000000000..d50eb44c71 --- /dev/null +++ b/tests/xmltester/tests/general/TestNGOverlayGC.xml @@ -0,0 +1,127 @@ + +Tests for OverlayNG operations with simple GeometryCollection inputs. +Covers topological situations with no precision collapse. +Uses a floating precision model. + + + + AgA - simple overlapping + + POLYGON ((2 8, 8 8, 8 2, 2 2, 2 8)) + + + GEOMETRYCOLLECTION ( POLYGON ((1 1, 1 5, 5 5, 5 1, 1 1)), POLYGON ((9 9, 9 5, 5 5, 5 9, 9 9)) ) + + + + MULTIPOLYGON (((2 2, 2 5, 5 5, 5 2, 2 2)), ((5 5, 5 8, 8 8, 8 5, 5 5))) + + + + + POLYGON ((1 1, 1 5, 2 5, 2 8, 5 8, 5 9, 9 9, 9 5, 8 5, 8 2, 5 2, 5 1, 1 1)) + + + + + MULTIPOLYGON (((5 8, 5 5, 2 5, 2 8, 5 8)), ((8 2, 5 2, 5 5, 8 5, 8 2))) + + + + + MULTIPOLYGON (((5 8, 5 5, 2 5, 2 8, 5 8)), ((8 8, 5 8, 5 9, 9 9, 9 5, 8 5, 8 8)), ((8 2, 5 2, 5 5, 8 5, 8 2)), ((2 2, 5 2, 5 1, 1 1, 1 5, 2 5, 2 2))) + + + + + + AgmP - simple partially overlapping + + POLYGON ((0 0, 0 1, 1 1, 0 0)) + + + GEOMETRYCOLLECTION ( MULTIPOINT ((0 0), (99 99)) ) + + + + POINT (0 0) + + + + + GEOMETRYCOLLECTION (POINT (99 99), POLYGON ((0 0, 0 1, 1 1, 0 0))) + + + + + POLYGON ((0 1, 1 1, 0 0, 0 1)) + + + + + GEOMETRYCOLLECTION (POLYGON ((0 1, 1 1, 0 0, 0 1)), POINT (99 99)) + + + + + + LgcA - simple overlapping + + LINESTRING (0 0, 10 10) + + + GEOMETRYCOLLECTION ( POLYGON ((1 1, 1 5, 5 5, 5 1, 1 1)), POLYGON ((9 9, 9 5, 5 5, 5 9, 9 9)) ) + + + + MULTILINESTRING ((1 1, 5 5), (5 5, 9 9)) + + + + + GEOMETRYCOLLECTION (LINESTRING (0 0, 1 1), LINESTRING (9 9, 10 10), POLYGON ((1 1, 1 5, 5 5, 5 1, 1 1)), POLYGON ((5 5, 5 9, 9 9, 9 5, 5 5))) + + + + + MULTILINESTRING ((0 0, 1 1), (9 9, 10 10)) + + + + + GEOMETRYCOLLECTION (POLYGON ((1 5, 5 5, 5 1, 1 1, 1 5)), POLYGON ((9 5, 5 5, 5 9, 9 9, 9 5)), LINESTRING (0 0, 1 1), LINESTRING (9 9, 10 10)) + + + + + + PgcA - simple covered + + POINT(5 5) + + + GEOMETRYCOLLECTION ( MULTIPOLYGON (((1 1, 1 5, 5 5, 5 1, 1 1)), ((9 9, 9 5, 5 5, 5 9, 9 9))) ) + + + + POINT (5 5) + + + + + MULTIPOLYGON (((1 1, 1 5, 5 5, 5 1, 1 1)), ((9 9, 9 5, 5 5, 5 9, 9 9))) + + + + + POINT EMPTY + + + + + MULTIPOLYGON (((1 5, 5 5, 5 1, 1 1, 1 5)), ((9 5, 5 5, 5 9, 9 9, 9 5))) + + + + +