[geos-commits] [SCM] GEOS branch master updated. edc36747409058cb10aeb8d9505d280dfbfcf697
git at osgeo.org
git at osgeo.org
Wed Feb 20 12:59:46 PST 2019
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, master has been updated
via edc36747409058cb10aeb8d9505d280dfbfcf697 (commit)
from b2347421b7b7faa6d77ecc615ec26d3cdb5d901c (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit edc36747409058cb10aeb8d9505d280dfbfcf697
Author: Martin Davis <mtnclimb at gmail.com>
Date: Wed Feb 20 12:56:13 2019 -0800
Port JTS InteriorPointArea improvements
diff --git a/include/geos/algorithm/InteriorPointArea.h b/include/geos/algorithm/InteriorPointArea.h
index 3a1bfd5..bc90f4c 100644
--- a/include/geos/algorithm/InteriorPointArea.h
+++ b/include/geos/algorithm/InteriorPointArea.h
@@ -3,9 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
- * Copyright (C) 2013 Sandro Santilli <strk at kbt.io>
- * Copyright (C) 2005-2006 Refractions Research Inc.
- * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ * Copyright (C) 2019 Martin Davis <mtnclimb at gmail.com>
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU Lesser General Public Licence as published
@@ -14,7 +12,8 @@
*
**********************************************************************
*
- * Last port: algorithm/InteriorPointArea.java r728 (JTS-1.13+)
+ * Last port: algorithm/InteriorPointArea.java (JTS-1.17+)
+ * https://github.com/locationtech/jts/commit/a140ca30cc51be4f65c950a30b0a8f51a6df75ba
*
**********************************************************************/
@@ -28,60 +27,59 @@
namespace geos {
namespace geom {
class Geometry;
-class LineString;
-class GeometryFactory;
-class GeometryCollection;
+class Polygon;
}
}
-
namespace geos {
namespace algorithm { // geos::algorithm
/** \brief
* Computes a point in the interior of an areal geometry.
+ * The point will lie in the geometry interior
+ * in all except certain pathological cases.
*
* <h2>Algorithm</h2>
+ * For each input polygon:
* <ul>
- * <li>Find a Y value which is close to the centre of
- * the geometry's vertical extent but is different
- * to any of it's Y ordinates.
- * <li>Create a horizontal bisector line using the Y value
- * and the geometry's horizontal extent
- * <li>Find the intersection between the geometry
- * and the horizontal bisector line.
- * The intersection is a collection of lines and points.
- * <li>Pick the midpoint of the largest intersection geometry
+ * <li>Determine a horizontal scan line on which the interior
+ * point will be located.
+ * To increase the chance of the scan line
+ * having non-zero-width intersection with the polygon
+ * the scan line Y ordinate is chosen to be near the centre of the polygon's
+ * Y extent but distinct from all of vertex Y ordinates.
+ * <li>Compute the sections of the scan line
+ * which lie in the interior of the polygon.
+ * <li>Choose the widest interior section
+ * and take its midpoint as the interior point.
* </ul>
+ * The final interior point is chosen as
+ * the one occurring in the widest interior section.
+ * <p>
+ * This algorithm is a tradeoff between performance
+ * and point quality (where points further from the geometry
+ * boundary are considered to be higher quality)
+ * Priority is given to performance.
+ * This means that the computed interior point
+ * may not be suitable for some uses
+ * (such as label positioning).
+ * <p>
+ * The algorithm handles some kinds of invalid/degenerate geometry,
+ * including zero-area and self-intersecting polygons.
+ * <p>
+ * Empty geometry is handled by returning a <code>null</code> point.
*
- * <b>
- * Note: If a fixed precision model is used,
- * in some cases this method may return a point
- * which does not lie in the interior.
- * </b>
+ * <h3>KNOWN BUGS</h3>
+ * <ul>
+ * <li>If a fixed precision model is used, in some cases this method may return
+ * a point which does not lie in the interior.
+ * <li>If the input polygon is <i>extremely</i> narrow the computed point
+ * may not lie in the interior of the polygon.
+ * </ul>
*/
class GEOS_DLL InteriorPointArea {
-private:
-
- bool foundInterior;
-
- const geom::GeometryFactory* factory;
-
- geom::Coordinate interiorPoint;
-
- double maxWidth;
-
- void add(const geom::Geometry* geom);
-
- const geom::Geometry* widestGeometry(const geom::Geometry* geometry);
-
- const geom::Geometry* widestGeometry(const geom::GeometryCollection* gc);
-
- geom::LineString* horizontalBisector(const geom::Geometry* geometry);
-
public:
-
/**
* Creates a new interior point finder
* for an areal geometry.
@@ -96,19 +94,17 @@ public:
* Gets the computed interior point.
*
* @return the coordinate of an interior point
+ * or <code>null</code> if the input geometry is empty
*/
bool getInteriorPoint(geom::Coordinate& ret) const;
private:
+ geom::Coordinate interiorPoint;
+ double maxWidth;
- /** \brief
- * Finds an interior point of a Polygon
- *
- * @param geometry the geometry to analyze
- * @return the midpoint of the largest intersection between the geometry and
- * a line halfway down its envelope
- */
- void addPolygon(const geom::Geometry* geometry);
+ void process(const geom::Geometry* geom);
+
+ void processPolygon(const geom::Polygon* polygon);
};
diff --git a/src/algorithm/InteriorPointArea.cpp b/src/algorithm/InteriorPointArea.cpp
index 5eaa674..ca335e4 100644
--- a/src/algorithm/InteriorPointArea.cpp
+++ b/src/algorithm/InteriorPointArea.cpp
@@ -3,9 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
- * Copyright (C) 2013 Sandro Santilli <strk at kbt.io>
- * Copyright (C) 2005-2006 Refractions Research Inc.
- * Copyright (C) 2001-2002 Vivid Solutions Inc.
+ * Copyright (C) 2019 Martin Davis <mtnclimb at gmail.com>
*
* This is free software; you can redistribute and/or modify it under
* the terms of the GNU Lesser General Public Licence as published
@@ -14,7 +12,8 @@
*
**********************************************************************
*
- * Last port: algorithm/InteriorPointArea.java r728 (JTS-1.13+)
+ * Last port: algorithm/InteriorPointArea.java (JTS-1.17+)
+ * https://github.com/locationtech/jts/commit/a140ca30cc51be4f65c950a30b0a8f51a6df75ba
*
**********************************************************************/
@@ -23,10 +22,12 @@
#include <geos/geom/Geometry.h>
#include <geos/geom/GeometryCollection.h>
#include <geos/geom/Polygon.h>
+#include <geos/geom/LinearRing.h>
#include <geos/geom/LineString.h>
#include <geos/geom/Envelope.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/CoordinateSequenceFactory.h>
+#include <geos/util/Interrupt.h>
#include <vector>
#include <typeinfo>
@@ -38,7 +39,7 @@ using namespace geos::geom;
namespace geos {
namespace algorithm { // geos.algorithm
-// file-statics
+// file statics
namespace {
double
@@ -48,24 +49,30 @@ avg(double a, double b)
}
/**
- * 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.
+ * Finds a safe scan line Y ordinate by projecting
+ * the polygon segments
+ * to the Y axis and finding the
+ * Y-axis interval which contains the centre of the Y extent.
+ * The centre of
+ * this interval is returned as the scan line Y-ordinate.
+ * <p>
+ * Note that in the case of (degenerate, invalid)
+ * zero-area polygons the computed Y value
+ * may be equal to a vertex Y-ordinate.
*
* @author mdavis
*
*/
-class SafeBisectorFinder {
+class ScanLineYOrdinateFinder {
public:
static double
- getBisectorY(const Polygon& poly)
+ getScanLineY(const Polygon& poly)
{
- SafeBisectorFinder finder(poly);
- return finder.getBisectorY();
+ ScanLineYOrdinateFinder finder(poly);
+ return finder.getScanLineY();
}
- SafeBisectorFinder(const Polygon& nPoly)
+
+ ScanLineYOrdinateFinder(const Polygon& nPoly)
: poly(nPoly)
{
// initialize using extremal values
@@ -75,7 +82,7 @@ public:
}
double
- getBisectorY()
+ getScanLineY()
{
process(*poly.getExteriorRing());
for(size_t i = 0; i < poly.getNumInteriorRing(); i++) {
@@ -85,7 +92,6 @@ public:
return bisectY;
}
-
private:
const Polygon& poly;
@@ -117,145 +123,220 @@ private:
}
}
}
+};
+
+class InteriorPointPolygon {
+public:
+ InteriorPointPolygon(const Polygon& poly)
+ : polygon(poly)
+ {
+ interiorPointY = ScanLineYOrdinateFinder::getScanLineY(polygon);
+ }
+
+ bool
+ getInteriorPoint(Coordinate& ret) const
+ {
+ ret = interiorPoint;
+ return true;
+ }
+
+ double getWidth()
+ {
+ return interiorSectionWidth;
+ }
+
+ void
+ process()
+ {
+ /**
+ * This results in returning a null Coordinate
+ */
+ if (polygon.isEmpty()) return;
+ /**
+ * set default interior point in case polygon has zero area
+ */
+ interiorPoint = *polygon.getCoordinate();
+
+ const LinearRing* shell = dynamic_cast<const LinearRing*>(polygon.getExteriorRing());
+ scanRing( *shell );
+ for (size_t i = 0; i < polygon.getNumInteriorRing(); i++) {
+ const LinearRing* hole = dynamic_cast<const LinearRing*>(polygon.getInteriorRingN(i));
+ scanRing( *hole );
+ }
+ findBestMidpoint(crossings);
+ }
+
+private:
+ const Polygon& polygon;
+ double interiorPointY;
+ double interiorSectionWidth = 0.0;
+ vector<double> crossings;
+ Coordinate interiorPoint;
+
+ void scanRing(const LinearRing& ring)
+ {
+ // skip rings which don't cross scan line
+ if ( ! intersectsHorizontalLine( ring.getEnvelopeInternal(), interiorPointY) )
+ return;
+
+ const CoordinateSequence* seq = ring.getCoordinatesRO();
+ for (size_t i = 1; i < seq->size(); i++) {
+ Coordinate ptPrev = seq->getAt(i - 1);
+ Coordinate pt = seq->getAt(i);
+ addEdgeCrossing(ptPrev, pt, interiorPointY, crossings);
+ }
+ }
+
+ void addEdgeCrossing(Coordinate& p0, Coordinate& p1, double scanY, vector<double>& crossings) {
+ // skip non-crossing segments
+ if ( !intersectsHorizontalLine(p0, p1, scanY) )
+ return;
+ if (! isEdgeCrossingCounted(p0, p1, scanY) )
+ return;
+
+ // edge intersects scan line, so add a crossing
+ double xInt = intersection(p0, p1, scanY);
+ crossings.push_back(xInt);
+ //checkIntersectionDD(p0, p1, scanY, xInt);
+ }
+
+ void findBestMidpoint(vector<double>& crossings)
+ {
+ // zero-area polygons will have no crossings
+ if (crossings.size() == 0) return;
+
+ //Assert.isTrue(0 == crossings.size() % 2, "Interior Point robustness failure: odd number of scanline crossings");
+
+ sort(crossings.begin(), crossings.end());
+ /*
+ * Entries in crossings list are expected to occur in pairs representing a
+ * section of the scan line interior to the polygon (which may be zero-length)
+ */
+ for (size_t i = 0; i < crossings.size(); i += 2) {
+ double x1 = crossings[i];
+ // crossings count must be even so this should be safe
+ double x2 = crossings[i + 1];
+
+ double width = x2 - x1;
+ if ( width > interiorSectionWidth ) {
+ interiorSectionWidth = width;
+ double interiorPointX = avg(x1, x2);
+ interiorPoint = Coordinate(interiorPointX, interiorPointY);
+ }
+ }
+ }
- SafeBisectorFinder(SafeBisectorFinder const&); /*= delete*/
- SafeBisectorFinder& operator=(SafeBisectorFinder const&); /*= delete*/
+ static bool
+ isEdgeCrossingCounted(Coordinate& p0, Coordinate& p1, double scanY) {
+ // skip horizontal lines
+ if ( p0.y == p1.y )
+ return false;
+ // handle cases where vertices lie on scan-line
+ // downward segment does not include start point
+ if ( p0.y == scanY && p1.y < scanY )
+ return false;
+ // upward segment does not include endpoint
+ if ( p1.y == scanY && p0.y < scanY )
+ return false;
+ return true;
+ }
+
+ static double
+ intersection(const Coordinate& p0, const Coordinate& p1, double Y)
+ {
+ double x0 = p0.x;
+ double x1 = p1.x;
+
+ if ( x0 == x1 )
+ return x0;
+
+ // Assert: segDX is non-zero, due to previous equality test
+ double segDX = x1 - x0;
+ double segDY = p1.y - p0.y;
+ double m = segDY / segDX;
+ double x = x0 + ((Y - p0.y) / m);
+ return x;
+ }
+
+ static bool
+ intersectsHorizontalLine(const Envelope* env, double y)
+ {
+ if ( y < env->getMinY() )
+ return false;
+ if ( y > env->getMaxY() )
+ return false;
+ return true;
+ }
+
+ static bool
+ intersectsHorizontalLine(const Coordinate& p0, const Coordinate& p1, double y) {
+ // both ends above?
+ if ( p0.y > y && p1.y > y )
+ return false;
+ // both ends below?
+ if ( p0.y < y && p1.y < y )
+ return false;
+ // segment must intersect line
+ return true;
+ }
};
} // anonymous namespace
-
-/*public*/
InteriorPointArea::InteriorPointArea(const Geometry* g)
{
- foundInterior = false;
- maxWidth = 0.0;
- factory = g->getFactory();
- add(g);
+ maxWidth = -1;
+ process(g);
}
-/*public*/
InteriorPointArea::~InteriorPointArea()
{
}
-/*public*/
bool
InteriorPointArea::getInteriorPoint(Coordinate& ret) const
{
- if(! foundInterior) {
+ // GEOS-specific code
+ if (maxWidth < 0)
return false;
- }
ret = interiorPoint;
return true;
}
-/*public*/
+/*private*/
void
-InteriorPointArea::add(const Geometry* geom)
+InteriorPointArea::process(const Geometry* geom)
{
+ if (geom->isEmpty())
+ return;
+
const Polygon* poly = dynamic_cast<const Polygon*>(geom);
if(poly) {
- addPolygon(geom);
+ processPolygon(poly);
return;
}
const GeometryCollection* gc = dynamic_cast<const GeometryCollection*>(geom);
if(gc) {
for(std::size_t i = 0, n = gc->getNumGeometries(); i < n; i++) {
- add(gc->getGeometryN(i));
+ process(gc->getGeometryN(i));
+ GEOS_CHECK_FOR_INTERRUPTS();
}
}
}
/*private*/
void
-InteriorPointArea::addPolygon(const Geometry* geometry)
+InteriorPointArea::processPolygon(const Polygon* polygon)
{
- if(geometry->isEmpty()) {
- return;
+ InteriorPointPolygon intPtPoly(*polygon);
+ intPtPoly.process();
+ double width = intPtPoly.getWidth();
+ if ( width > maxWidth ) {
+ maxWidth = width;
+ intPtPoly.getInteriorPoint(interiorPoint);
}
-
- Coordinate intPt;
- double width;
-
- unique_ptr<LineString> bisector(horizontalBisector(geometry));
- if(bisector->getLength() == 0.0) {
- width = 0;
- intPt = bisector->getCoordinateN(0);
- }
- else {
- unique_ptr<Geometry> intersections(bisector->intersection(geometry));
- const Geometry* widestIntersection = widestGeometry(intersections.get());
- const Envelope* env = widestIntersection->getEnvelopeInternal();
- width = env->getWidth();
- env->centre(intPt);
- }
- if(!foundInterior || width > maxWidth) {
- interiorPoint = intPt;
- maxWidth = width;
- foundInterior = true;
- }
-}
-
-//@return if geometry is a collection, the widest sub-geometry; otherwise,
-//the geometry itself
-const Geometry*
-InteriorPointArea::widestGeometry(const Geometry* geometry)
-{
- const GeometryCollection* gc = dynamic_cast<const GeometryCollection*>(geometry);
- if(gc) {
- return widestGeometry(gc);
- }
- else {
- return geometry;
- }
-}
-
-const Geometry*
-InteriorPointArea::widestGeometry(const GeometryCollection* gc)
-{
- if(gc->isEmpty()) {
- return gc;
- }
- const Geometry* p_widestGeometry = gc->getGeometryN(0);
-
- // scan remaining geom components to see if any are wider
- for(std::size_t i = 1, n = gc->getNumGeometries(); i < n; i++) { // start at 1
- const Envelope* env1(gc->getGeometryN(i)->getEnvelopeInternal());
- const Envelope* env2(p_widestGeometry->getEnvelopeInternal());
- if(env1->getWidth() > env2->getWidth()) {
- p_widestGeometry = gc->getGeometryN(i);
- }
- }
- return p_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 bisectY = SafeBisectorFinder::getBisectorY(*dynamic_cast<const Polygon*>(geometry));
- vector<Coordinate>* cv = new vector<Coordinate>(2);
- (*cv)[0].x = envelope->getMinX();
- (*cv)[0].y = bisectY;
- (*cv)[1].x = envelope->getMaxX();
- (*cv)[1].y = bisectY;
-
- CoordinateSequence* cl = factory->getCoordinateSequenceFactory()->create(cv);
-
- LineString* ret = factory->createLineString(cl);
- return ret;
}
} // namespace geos.algorithm
diff --git a/tests/perf/algorithm/InteriorPointAreaPerfTest.cpp b/tests/perf/algorithm/InteriorPointAreaPerfTest.cpp
index c118439..9ff3913 100644
--- a/tests/perf/algorithm/InteriorPointAreaPerfTest.cpp
+++ b/tests/perf/algorithm/InteriorPointAreaPerfTest.cpp
@@ -59,7 +59,8 @@ public:
* algorithm, and provides a more realistic test.
*/
using geos::precision::SimpleGeometryPrecisionReducer;
- PrecisionModel p_pm(SIZE);
+ double scale = nPts / SIZE;
+ PrecisionModel p_pm(scale);
SimpleGeometryPrecisionReducer reducer(&p_pm);
std::unique_ptr<Geometry> sinePolyCrinkly(reducer.reduce(sinePoly.get()));
sinePoly.reset();
@@ -106,7 +107,6 @@ private:
sw.stop();
cout << poly.getNumPoints() << " points: " << sw.getTotFormatted() << endl;
-
}
std::unique_ptr<Polygon>
@@ -131,6 +131,7 @@ main()
{
InteriorPointAreaPerfTest tester;
+ tester.test(100);
tester.test(1000);
tester.test(10000);
tester.test(100000);
diff --git a/tests/unit/algorithm/InteriorPointAreaTest.cpp b/tests/unit/algorithm/InteriorPointAreaTest.cpp
index fa9a82c..32cd3c8 100644
--- a/tests/unit/algorithm/InteriorPointAreaTest.cpp
+++ b/tests/unit/algorithm/InteriorPointAreaTest.cpp
@@ -3,6 +3,7 @@
* GEOS - Geometry Engine Open Source
* http://geos.osgeo.org
*
+ * Copyright (C) 2019 Martin Davis <mtnclimb at gmail.com>
* Copyright (C) 2011 Sandro Santilli <strk at kbt.io>
*
* This is free software; you can redistribute and/or modify it under
@@ -55,29 +56,24 @@ group test_interiorpointarea_group("geos::algorithm::InteriorPointArea");
//
// http://trac.osgeo.org/geos/ticket/475
-// This is a test for a memory leak more than anything else
+// This test no longer throws, since invalid inputs are now handled
template<>
template<>
void object::test<1>
()
{
- Coordinate centroid;
+ Coordinate result;
- // this polygon is a typical hourglass-shape with a self intersection
+ // invalid polygon - classic hourglass-shape with a self intersection
// without a node
geom.reset(reader.read("POLYGON((6 54, 15 54, 6 47, 15 47, 6 54))"));
- bool threw = false;
+ InteriorPointArea interiorPointArea(geom.get());
+ interiorPointArea.getInteriorPoint(result);
- try {
- InteriorPointArea interior_builder(geom.get());
- interior_builder.getInteriorPoint(centroid);
- }
- catch(...) {
- threw = true;
- }
+ geos::geom::Coordinate expected(6, 54);
- ensure(threw);
+ ensure_equals(result, expected);
}
} // namespace tut
diff --git a/tests/unit/capi/GEOSPointOnSurfaceTest.cpp b/tests/unit/capi/GEOSPointOnSurfaceTest.cpp
index e4899c0..2957558 100644
--- a/tests/unit/capi/GEOSPointOnSurfaceTest.cpp
+++ b/tests/unit/capi/GEOSPointOnSurfaceTest.cpp
@@ -150,7 +150,7 @@ void object::test<4>
wkt_ = GEOSWKTWriter_write(wktw_, geom2_);
- ensure_equals(std::string(wkt_), std::string("POINT (56.528917 25.210417)"));
+ ensure_equals(std::string(wkt_), std::string("POINT (56.528667 25.210167)"));
}
-----------------------------------------------------------------------
Summary of changes:
include/geos/algorithm/InteriorPointArea.h | 94 +++----
src/algorithm/InteriorPointArea.cpp | 313 +++++++++++++--------
tests/perf/algorithm/InteriorPointAreaPerfTest.cpp | 5 +-
tests/unit/algorithm/InteriorPointAreaTest.cpp | 20 +-
tests/unit/capi/GEOSPointOnSurfaceTest.cpp | 2 +-
5 files changed, 254 insertions(+), 180 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list