Implement GEOSMinimumClearance and GEOSMinimumClearanceLine

Includes tests and C-API exposure.

Patch by Daniel Baston <dbaston@maponics.com>
via #65

Closes #776

git-svn-id: http://svn.osgeo.org/geos/trunk@4203 5242fede-7e19-0410-aef8-94bd7d2200fb

NEWS

@@ -6,6 +6,9 @@ Changes in 3.6.0
   - CAPI: GEOSGeom_getPrecision - PHP: Geometry->getPrecision
   - CAPI: GEOSMinimumRotatedRectangle and GEOSMinimumWidth
     (#729, Nyall Dawson)
+  - CAPI: GEOSSTRtree_nearest (#768, Dan Baston)
+  - CAPI: GEOSMinimumClearance and GEOSMinimumClearanceLine
+    (#776, Dan Baston)
 - Improvements:
   - ...
 - C++ API changes:

capi/geos_c.cpp

@@ -448,6 +448,18 @@ GEOSMinimumWidth(const Geometry *g)
     return GEOSMinimumWidth_r( handle, g );
 }
 
+Geometry *
+GEOSMinimumClearanceLine(const Geometry *g)
+{
+    return GEOSMinimumClearanceLine_r( handle, g );
+}
+
+int
+GEOSMinimumClearance(const Geometry *g, double *d)
+{
+    return GEOSMinimumClearance_r( handle, g, d);
+}
+
 Geometry *
 GEOSDifference(const Geometry *g1, const Geometry *g2)
 {

capi/geos_c.h.in

@@ -557,6 +557,13 @@ extern GEOSGeometry GEOS_DLL *GEOSMinimumRotatedRectangle_r(GEOSContextHandle_t
 extern GEOSGeometry GEOS_DLL *GEOSMinimumWidth_r(GEOSContextHandle_t handle,
                                                const GEOSGeometry* g);
 
+extern GEOSGeometry GEOS_DLL *GEOSMinimumClearanceLine_r(GEOSContextHandle_t handle,
+                                                     const GEOSGeometry* g);
+
+extern int GEOS_DLL GEOSMinimumClearance_r(GEOSContextHandle_t handle,
+                                           const GEOSGeometry* g,
+                                           double* distance);
+
 extern GEOSGeometry GEOS_DLL *GEOSDifference_r(GEOSContextHandle_t handle,
                                                const GEOSGeometry* g1,
                                                const GEOSGeometry* g2);
@@ -1496,6 +1503,32 @@ extern GEOSGeometry GEOS_DLL *GEOSMinimumRotatedRectangle(const GEOSGeometry* g)
  */
 extern GEOSGeometry GEOS_DLL *GEOSMinimumWidth(const GEOSGeometry* g);
 
+/* Computes the minimum clearance of a geometry.  The minimum clearance is the smallest amount by which
+ * a vertex could be move to produce an invalid polygon, a non-simple linestring, or a multipoint with
+ * repeated points.  If a geometry has a minimum clearance of 'eps', it can be said that:
+ *
+ * -  No two distinct vertices in the geometry are separated by less than 'eps'
+ * -  No vertex is closer than 'eps' to a line segment of which it is not an endpoint.
+ *
+ * If the minimum clearance cannot be defined for a geometry (such as with a single point, or a multipoint
+ * whose points are identical, a value of Infinity will be calculated.
+ *
+ * @param g the input geometry
+ * @param d a double to which the result can be stored
+ *
+ * @return 0 if no exception occurred
+ *         2 if an exception occurred
+ */
+extern int GEOS_DLL GEOSMinimumClearance(const GEOSGeometry* g, double* d);
+
+/* Returns a LineString whose endpoints define the minimum clearance of a geometry.
+ * If the geometry has no minimum clearance, an empty LineString will be returned.
+ *
+ * @param g the input geometry
+ * @return a LineString, or NULL if an exception occurred.
+ */
+extern GEOSGeometry GEOS_DLL *GEOSMinimumClearanceLine(const GEOSGeometry* g);
+
 extern GEOSGeometry GEOS_DLL *GEOSDifference(const GEOSGeometry* g1, const GEOSGeometry* g2);
 extern GEOSGeometry GEOS_DLL *GEOSSymDifference(const GEOSGeometry* g1, const GEOSGeometry* g2);
 extern GEOSGeometry GEOS_DLL *GEOSBoundary(const GEOSGeometry* g);

capi/geos_ts_c.cpp

@@ -116,6 +116,7 @@ using namespace std;
 #undef VERBOSE_EXCEPTIONS
 
 #include <geos/export.h>
+#include <geos/precision/MinimumClearance.h>
 
 
 // import the most frequently used definitions globally
@@ -2078,6 +2079,72 @@ GEOSMinimumWidth_r(GEOSContextHandle_t extHandle, const Geometry *g)
     return NULL;
 }
 
+Geometry *
+GEOSMinimumClearanceLine_r(GEOSContextHandle_t extHandle, const Geometry *g)
+{
+    if ( 0 == extHandle )
+    {
+        return NULL;
+    }
+
+    GEOSContextHandleInternal_t *handle = 0;
+    handle = reinterpret_cast<GEOSContextHandleInternal_t*>(extHandle);
+    if ( 0 == handle->initialized )
+    {
+        return NULL;
+    }
+
+    try
+    {
+        geos::precision::MinimumClearance mc(g);
+        return mc.getLine().release();
+    }
+    catch (const std::exception &e)
+    {
+        handle->ERROR_MESSAGE("%s", e.what());
+    }
+    catch (...)
+    {
+        handle->ERROR_MESSAGE("Unknown exception thrown");
+    }
+
+    return NULL;
+}
+
+int
+GEOSMinimumClearance_r(GEOSContextHandle_t extHandle, const Geometry *g, double *d)
+{
+    if ( 0 == extHandle )
+    {
+        return 2;
+    }
+
+    GEOSContextHandleInternal_t *handle = 0;
+    handle = reinterpret_cast<GEOSContextHandleInternal_t*>(extHandle);
+    if ( 0 == handle->initialized )
+    {
+        return 2;
+    }
+
+    try
+    {
+        geos::precision::MinimumClearance mc(g);
+        double res = mc.getDistance();
+        *d = res;
+        return 0;
+    }
+    catch (const std::exception &e)
+    {
+        handle->ERROR_MESSAGE("%s", e.what());
+    }
+    catch (...)
+    {
+        handle->ERROR_MESSAGE("Unknown exception thrown");
+    }
+
+    return 2;
+}
+
 
 Geometry *
 GEOSDifference_r(GEOSContextHandle_t extHandle, const Geometry *g1, const Geometry *g2)

include/geos/index/strtree/BoundablePair.h

@@ -78,7 +78,7 @@ class BoundablePair {
 		 *
 		 * @return
 		 */
-		double distance();
+		double distance() const;
 
 		/**
 		 * Gets the minimum possible distance between the Boundables in

include/geos/index/strtree/STRtree.h

@@ -140,8 +140,9 @@ using AbstractSTRtree::query;
 	}
 
 	const void* nearestNeighbour(const geom::Envelope *env, const void* item, ItemDistance* itemDist);
-	const void* nearestNeighbour(BoundablePair* initBndPair);
-	const void* nearestNeighbour(BoundablePair* initBndPair, double maxDistance);
+	std::pair<const void*, const void*> nearestNeighbour(BoundablePair* initBndPair);
+	std::pair<const void*, const void*> nearestNeighbour(ItemDistance* itemDist);
+	std::pair<const void*, const void*> nearestNeighbour(BoundablePair* initBndPair, double maxDistance);
 
 	bool remove(const geom::Envelope *itemEnv, void* item) {
 		return AbstractSTRtree::remove(itemEnv, item);

include/geos/operation/distance/FacetSequence.h

@@ -0,0 +1,66 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2016 Daniel Baston
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: operation/distance/FacetSequence.java (JTS-1.14)
+ *
+ **********************************************************************/
+
+#ifndef GEOS_OPERATION_DISTANCE_FACETSEQUENCE_H
+#define GEOS_OPERATION_DISTANCE_FACETSEQUENCE_H
+
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Coordinate.h>
+
+using namespace geos::geom;
+
+namespace geos {
+    namespace operation {
+        namespace distance {
+            class FacetSequence {
+            private:
+                const CoordinateSequence *pts;
+                const size_t start;
+                const size_t end;
+
+                /* Unlike JTS, we store the envelope in the FacetSequence so that it has a clear owner.  This is
+                 * helpful when making a tree of FacetSequence objects (FacetSequenceTreeBuilder)
+                 * */
+                Envelope env;
+
+                double computeLineLineDistance(const FacetSequence & facetSeq) const;
+
+                double computePointLineDistance(const Coordinate & pt, const FacetSequence & facetSeq) const;
+
+                void computeEnvelope();
+
+            public:
+                const Envelope * getEnvelope() const;
+
+                const Coordinate * getCoordinate(size_t index) const;
+
+                size_t size() const;
+
+                bool isPoint() const;
+
+                double distance(const FacetSequence & facetSeq);
+
+                FacetSequence(const CoordinateSequence *pts, size_t start, size_t end);
+            };
+
+        }
+    }
+}
+
+#endif //GEOS_OPERATION_DISTANCE_FACETSEQUENCE_H

include/geos/operation/distance/FacetSequenceTreeBuilder.h

@@ -0,0 +1,52 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2016 Daniel Baston
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: operation/distance/FacetSequenceTreeBuilder.java (JTS-1.14)
+ *
+ **********************************************************************/
+
+#ifndef GEOS_OPERATION_DISTANCE_FACETSEQUENCETREEBUILDER_H
+#define GEOS_OPERATION_DISTANCE_FACETSEQUENCETREEBUILDER_H
+
+#include <geos/index/strtree/STRtree.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/operation/distance/FacetSequence.h>
+
+using namespace geos::geom;
+using namespace geos::index::strtree;
+using namespace geos::operation::distance;
+
+namespace geos {
+    namespace operation {
+        namespace distance {
+            class FacetSequenceTreeBuilder {
+            private:
+                // 6 seems to be a good facet sequence size
+                static const int FACET_SEQUENCE_SIZE = 6;
+
+                // Seems to be better to use a minimum node capacity
+                static const int STR_TREE_NODE_CAPACITY = 4;
+
+                static void addFacetSequences(const CoordinateSequence* pts, std::vector<FacetSequence*> & sections);
+                static std::vector<FacetSequence*> * computeFacetSequences(const Geometry* g);
+
+            public:
+                static STRtree* build(const Geometry* g);
+            };
+        }
+    }
+}
+
+#endif //GEOS_FACETSEQUENCETREEBUILDER_H

include/geos/operation/distance/Makefile.am

@@ -11,4 +11,6 @@ geos_HEADERS = \
 	ConnectedElementLocationFilter.h \
 	ConnectedElementPointFilter.h \
 	DistanceOp.h \
+	FacetSequence.h \
+	FacetSequenceTreeBuilder.h \
 	GeometryLocation.h

include/geos/precision/Makefile.am

@@ -13,5 +13,6 @@ geos_HEADERS = \
     CommonBitsRemover.h \
     EnhancedPrecisionOp.h \
     GeometryPrecisionReducer.h \
+    MinimumClearance.h \
     PrecisionReducerCoordinateOperation.h \
     SimpleGeometryPrecisionReducer.h

include/geos/precision/MinimumClearance.h

@@ -0,0 +1,60 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2016 Daniel Baston
+ *
+ * This is free software; you can redistribute and/or modify it under
+ * the terms of the GNU Lesser General Public Licence as published
+ * by the Free Software Foundation.
+ * See the COPYING file for more information.
+ *
+ **********************************************************************
+ *
+ * Last port: precision/MinimumClearance.java (JTS-1.14)
+ *
+ **********************************************************************/
+
+#ifndef GEOS_PRECISION_MINIMUMCLEARANCE_H
+#define GEOS_PRECISION_MINIMUMCLEARANCE_H
+
+#include <geos/geom/Geometry.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/CoordinateSequence.h>
+
+namespace geos {
+namespace precision {
+class MinimumClearance {
+    private:
+        const geom::Geometry* inputGeom;
+        double minClearance;
+        std::auto_ptr<geom::CoordinateSequence> minClearancePts;
+
+        void compute();
+    public:
+        MinimumClearance(const geom::Geometry* g);
+
+        /**
+         * Gets the Minimum Clearance distance.
+         *
+         * @return the value of the minimum clearance distance
+         * or <tt>DBL_MAX</tt> if no Minimum Clearance distance exists
+         */
+        double getDistance();
+
+        /**
+         * Gets a LineString containing two points
+         * which are at the Minimum Clearance distance.
+         *
+         * @return the value of the minimum clearance distance
+         * or <tt>LINESTRING EMPTY</tt> if no Minimum Clearance distance exists
+         */
+        std::auto_ptr<geom::LineString> getLine();
+};
+}
+}
+
+#endif
+
+