From 3ab8ecce402ef9de35fea5e362a3a73910c81e45 Mon Sep 17 00:00:00 2001 From: strk Date: Fri, 8 Mar 2013 13:02:34 +0000 Subject: [PATCH] Fix GEOSPointOnSurface returning point on boundary (#623) Ports SafeBisector for InteriorPointArea from JTS git-svn-id: http://svn.osgeo.org/geos/trunk@3781 5242fede-7e19-0410-aef8-94bd7d2200fb --- include/geos/algorithm/InteriorPointArea.h | 30 ++++++-- src/algorithm/Centroid.cpp | 1 + src/algorithm/InteriorPointArea.cpp | 87 ++++++++++++++++++++-- tests/unit/capi/GEOSPointOnSurfaceTest.cpp | 2 +- .../xmltester/tests/general/TestInteriorPoint.xml | 9 ++- 5 files changed, 114 insertions(+), 15 deletions(-) diff --git a/include/geos/algorithm/InteriorPointArea.h b/include/geos/algorithm/InteriorPointArea.h index 8447a33d..7ca3064c 100644 --- a/include/geos/algorithm/InteriorPointArea.h +++ b/include/geos/algorithm/InteriorPointArea.h @@ -39,14 +39,19 @@ namespace geos { namespace algorithm { // geos::algorithm /** \brief - * Computes a point in the interior of an area geometry. + * Computes a point in the interior of an areal geometry. * *

Algorithm

* * * @@ -77,14 +82,27 @@ private: public: + /** + * Creates a new interior point finder + * for an areal geometry. + * + * @param g an areal geometry + */ InteriorPointArea(const geom::Geometry *g); ~InteriorPointArea(); + /** + * Gets the computed interior point. + * + * @return the coordinate of an interior point + */ bool getInteriorPoint(geom::Coordinate& ret) const; +private: + /** \brief - * Finds a reasonable point at which to label a Geometry. + * Finds an interior point of a Polygon * * @param geometry the geometry to analyze * @return the midpoint of the largest intersection between the geometry and diff --git a/src/algorithm/Centroid.cpp b/src/algorithm/Centroid.cpp index bae77925..4871b403 100644 --- a/src/algorithm/Centroid.cpp +++ b/src/algorithm/Centroid.cpp @@ -25,6 +25,7 @@ #include #include + #include // for std::abs using namespace geos::geom; diff --git a/src/algorithm/InteriorPointArea.cpp b/src/algorithm/InteriorPointArea.cpp index fd87263f..8d504760 100644 --- a/src/algorithm/InteriorPointArea.cpp +++ b/src/algorithm/InteriorPointArea.cpp @@ -43,7 +43,73 @@ namespace { double avg(double a, double b){return (a+b)/2.0;} -} + /** + * Finds a safe bisector Y ordinate + * by projecting to the Y axis + * and finding the Y-ordinate interval + * which contains the centre of the Y extent. + * The centre of this interval is returned as the bisector Y-ordinate. + * + * @author mdavis + * + */ + class SafeBisectorFinder + { + public: + static double getBisectorY(const Polygon& poly) + { + SafeBisectorFinder finder(poly); + return finder.getBisectorY(); + } + SafeBisectorFinder(const Polygon& nPoly) + : poly(nPoly) + { + // initialize using extremal values + hiY = poly.getEnvelopeInternal()->getMaxY(); + loY = poly.getEnvelopeInternal()->getMinY(); + centreY = avg(loY, hiY); + } + + double getBisectorY() + { + process(*poly.getExteriorRing()); + for (size_t i = 0; i < poly.getNumInteriorRing(); i++) { + process(*poly.getInteriorRingN(i)); + } + double bisectY = avg(hiY, loY); + return bisectY; + } + + + private: + const Polygon& poly; + + double centreY; + double hiY; + double loY; + + void process(const LineString& line) { + const CoordinateSequence* seq = line.getCoordinatesRO(); + for (int i = 0, s = seq->size(); i < s; i++) { + double y = seq->getY(i); + updateInterval(y); + } + } + + void updateInterval(double y) { + if (y <= centreY) { + if (y > loY) + loY = y; + } + else if (y > centreY) { + if (y < hiY) { + hiY = y; + } + } + } + }; + +} // anonymous namespace /*public*/ @@ -89,7 +155,7 @@ InteriorPointArea::add(const Geometry *geom) } } -/*public*/ +/*private*/ void InteriorPointArea::addPolygon(const Geometry *geometry) { @@ -137,8 +203,8 @@ InteriorPointArea::widestGeometry(const GeometryCollection* gc) { } const Geometry* widestGeometry=gc->getGeometryN(0); - //Start at 1 - for(std::size_t i=1, n=gc->getNumGeometries(); igetNumGeometries(); igetGeometryN(i)->getEnvelopeInternal()); const Envelope *env2(widestGeometry->getEnvelopeInternal()); @@ -149,18 +215,25 @@ InteriorPointArea::widestGeometry(const GeometryCollection* gc) { return widestGeometry; } +/* private */ LineString* InteriorPointArea::horizontalBisector(const Geometry *geometry) { const Envelope *envelope=geometry->getEnvelopeInternal(); + + /** + * Original algorithm. Fails when geometry contains a horizontal + * segment at the Y midpoint. + */ // Assert: for areas, minx <> maxx - double avgY=avg(envelope->getMinY(),envelope->getMaxY()); + //double avgY=avg(envelope->getMinY(),envelope->getMaxY()); + double bisectY = SafeBisectorFinder::getBisectorY(*dynamic_cast(geometry)); vector*cv=new vector(2); (*cv)[0].x = envelope->getMinX(); - (*cv)[0].y = avgY; + (*cv)[0].y = bisectY; (*cv)[1].x = envelope->getMaxX(); - (*cv)[1].y = avgY; + (*cv)[1].y = bisectY; CoordinateSequence *cl = factory->getCoordinateSequenceFactory()->create(cv); diff --git a/tests/unit/capi/GEOSPointOnSurfaceTest.cpp b/tests/unit/capi/GEOSPointOnSurfaceTest.cpp index 0849ff4c..92100d8a 100644 --- a/tests/unit/capi/GEOSPointOnSurfaceTest.cpp +++ b/tests/unit/capi/GEOSPointOnSurfaceTest.cpp @@ -147,7 +147,7 @@ namespace tut wkt_ = GEOSWKTWriter_write(wktw_, geom2_); - ensure_equals(std::string(wkt_), std::string( "POINT (56.528833 25.210333)" ) ); + ensure_equals(std::string(wkt_), std::string( "POINT (56.528917 25.210417)" ) ); } diff --git a/tests/xmltester/tests/general/TestInteriorPoint.xml b/tests/xmltester/tests/general/TestInteriorPoint.xml index fc5944bb..0d765cac 100644 --- a/tests/xmltester/tests/general/TestInteriorPoint.xml +++ b/tests/xmltester/tests/general/TestInteriorPoint.xml @@ -72,7 +72,14 @@ A - polygon with horizontal segment at centre (L shape) POLYGON ((0 2, 0 4, 6 4, 6 0, 2 0, 2 2, 0 2)) - POINT (4 2) + POINT (3 3) + + + + A - polygon with horizontal segment at centre (narrower L shape) + POLYGON ((0 2, 0 4, 3 4, 3 0, 2 0, 2 2, 0 2)) + + POINT (2 3) -- 2.11.4.GIT