[geos-commits] r3781 - in trunk: include/geos/algorithm src/algorithm tests/unit/capi tests/xmltester/tests/general
svn_geos at osgeo.org
svn_geos at osgeo.org
Fri Mar 8 05:02:34 PST 2013
Author: strk
Date: 2013-03-08 05:02:34 -0800 (Fri, 08 Mar 2013)
New Revision: 3781
Modified:
trunk/include/geos/algorithm/InteriorPointArea.h
trunk/src/algorithm/Centroid.cpp
trunk/src/algorithm/InteriorPointArea.cpp
trunk/tests/unit/capi/GEOSPointOnSurfaceTest.cpp
trunk/tests/xmltester/tests/general/TestInteriorPoint.xml
Log:
Fix GEOSPointOnSurface returning point on boundary (#623)
Ports SafeBisector for InteriorPointArea from JTS
Modified: trunk/include/geos/algorithm/InteriorPointArea.h
===================================================================
--- trunk/include/geos/algorithm/InteriorPointArea.h 2013-03-05 14:55:16 UTC (rev 3780)
+++ trunk/include/geos/algorithm/InteriorPointArea.h 2013-03-08 13:02:34 UTC (rev 3781)
@@ -39,14 +39,19 @@
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.
*
* <h2>Algorithm</h2>
* <ul>
- * <li>Find the intersections between the geometry
- * and the horizontal bisector of the area's envelope
- * <li>Pick the midpoint of the largest intersection (the intersections
- * will be lines and points)
+ * <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
* </ul>
*
* <b>
@@ -77,14 +82,27 @@
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
Modified: trunk/src/algorithm/Centroid.cpp
===================================================================
--- trunk/src/algorithm/Centroid.cpp 2013-03-05 14:55:16 UTC (rev 3780)
+++ trunk/src/algorithm/Centroid.cpp 2013-03-08 13:02:34 UTC (rev 3781)
@@ -25,6 +25,7 @@
#include <geos/geom/GeometryCollection.h>
#include <geos/geom/LineString.h>
+
#include <cmath> // for std::abs
using namespace geos::geom;
Modified: trunk/src/algorithm/InteriorPointArea.cpp
===================================================================
--- trunk/src/algorithm/InteriorPointArea.cpp 2013-03-05 14:55:16 UTC (rev 3780)
+++ trunk/src/algorithm/InteriorPointArea.cpp 2013-03-08 13:02:34 UTC (rev 3781)
@@ -43,9 +43,75 @@
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*/
InteriorPointArea::InteriorPointArea(const Geometry *g)
{
@@ -89,7 +155,7 @@
}
}
-/*public*/
+/*private*/
void
InteriorPointArea::addPolygon(const Geometry *geometry)
{
@@ -137,8 +203,8 @@
}
const Geometry* widestGeometry=gc->getGeometryN(0);
- //Start at 1
- for(std::size_t i=1, n=gc->getNumGeometries(); i<n; i++)
+ // 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(widestGeometry->getEnvelopeInternal());
@@ -149,18 +215,25 @@
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<const Polygon *>(geometry));
vector<Coordinate>*cv=new vector<Coordinate>(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);
Modified: trunk/tests/unit/capi/GEOSPointOnSurfaceTest.cpp
===================================================================
--- trunk/tests/unit/capi/GEOSPointOnSurfaceTest.cpp 2013-03-05 14:55:16 UTC (rev 3780)
+++ trunk/tests/unit/capi/GEOSPointOnSurfaceTest.cpp 2013-03-08 13:02:34 UTC (rev 3781)
@@ -147,7 +147,7 @@
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)" ) );
}
Modified: trunk/tests/xmltester/tests/general/TestInteriorPoint.xml
===================================================================
--- trunk/tests/xmltester/tests/general/TestInteriorPoint.xml 2013-03-05 14:55:16 UTC (rev 3780)
+++ trunk/tests/xmltester/tests/general/TestInteriorPoint.xml 2013-03-08 13:02:34 UTC (rev 3781)
@@ -72,10 +72,17 @@
<desc>A - polygon with horizontal segment at centre (L shape)</desc>
<a> POLYGON ((0 2, 0 4, 6 4, 6 0, 2 0, 2 2, 0 2))
</a>
-<test><op name="getInteriorPoint" arg1="A" > POINT (4 2) </op></test>
+<test><op name="getInteriorPoint" arg1="A" > POINT (3 3) </op></test>
</case>
<case>
+ <desc>A - polygon with horizontal segment at centre (narrower L shape)</desc>
+ <a> POLYGON ((0 2, 0 4, 3 4, 3 0, 2 0, 2 2, 0 2))
+ </a>
+<test><op name="getInteriorPoint" arg1="A" > POINT (2 3) </op></test>
+</case>
+
+<case>
<desc>mA - polygons with holes</desc>
<a> MULTIPOLYGON (((50 260, 240 340, 260 100, 20 60, 90 140, 50 260), (200 280, 140 240, 180 160, 240 140, 200 280)), ((380 280, 300 260, 340 100, 440 80, 380 280), (380 220, 340 200, 400 100, 380 220)))
</a>
More information about the geos-commits
mailing list