[geos-commits] r4203 - in trunk: . capi include/geos/index/strtree include/geos/operation/distance include/geos/precision src/index/strtree src/operation/distance src/precision tests/unit tests/unit/capi
Sandro Santilli
strk at keybit.net
Mon Apr 25 11:06:50 PDT 2016
Author: strk
Date: 2016-04-25 11:06:50 -0700 (Mon, 25 Apr 2016)
New Revision: 4203
Added:
trunk/include/geos/operation/distance/FacetSequence.h
trunk/include/geos/operation/distance/FacetSequenceTreeBuilder.h
trunk/include/geos/precision/MinimumClearance.h
trunk/src/operation/distance/FacetSequence.cpp
trunk/src/operation/distance/FacetSequenceTreeBuilder.cpp
trunk/src/precision/MinimumClearance.cpp
trunk/tests/unit/capi/GEOSMinimumClearanceTest.cpp
Modified:
trunk/NEWS
trunk/capi/geos_c.cpp
trunk/capi/geos_c.h.in
trunk/capi/geos_ts_c.cpp
trunk/include/geos/index/strtree/BoundablePair.h
trunk/include/geos/index/strtree/STRtree.h
trunk/include/geos/operation/distance/Makefile.am
trunk/include/geos/precision/Makefile.am
trunk/src/index/strtree/BoundablePair.cpp
trunk/src/index/strtree/STRtree.cpp
trunk/src/operation/distance/Makefile.am
trunk/src/precision/Makefile.am
trunk/tests/unit/Makefile.am
Log:
Implement GEOSMinimumClearance and GEOSMinimumClearanceLine
Includes tests and C-API exposure.
Patch by Daniel Baston <dbaston at maponics.com>
via https://github.com/libgeos/libgeos/pull/65
Closes #776
Modified: trunk/NEWS
===================================================================
--- trunk/NEWS 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/NEWS 2016-04-25 18:06:50 UTC (rev 4203)
@@ -6,6 +6,9 @@
- 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:
Modified: trunk/capi/geos_c.cpp
===================================================================
--- trunk/capi/geos_c.cpp 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/capi/geos_c.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -449,6 +449,18 @@
}
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)
{
return GEOSDifference_r( handle, g1, g2 );
Modified: trunk/capi/geos_c.h.in
===================================================================
--- trunk/capi/geos_c.h.in 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/capi/geos_c.h.in 2016-04-25 18:06:50 UTC (rev 4203)
@@ -557,6 +557,13 @@
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 *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);
Modified: trunk/capi/geos_ts_c.cpp
===================================================================
--- trunk/capi/geos_ts_c.cpp 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/capi/geos_ts_c.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -116,6 +116,7 @@
#undef VERBOSE_EXCEPTIONS
#include <geos/export.h>
+#include <geos/precision/MinimumClearance.h>
// import the most frequently used definitions globally
@@ -2078,7 +2079,73 @@
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)
{
Modified: trunk/include/geos/index/strtree/BoundablePair.h
===================================================================
--- trunk/include/geos/index/strtree/BoundablePair.h 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/include/geos/index/strtree/BoundablePair.h 2016-04-25 18:06:50 UTC (rev 4203)
@@ -78,7 +78,7 @@
*
* @return
*/
- double distance();
+ double distance() const;
/**
* Gets the minimum possible distance between the Boundables in
Modified: trunk/include/geos/index/strtree/STRtree.h
===================================================================
--- trunk/include/geos/index/strtree/STRtree.h 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/include/geos/index/strtree/STRtree.h 2016-04-25 18:06:50 UTC (rev 4203)
@@ -140,8 +140,9 @@
}
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);
Added: trunk/include/geos/operation/distance/FacetSequence.h
===================================================================
--- trunk/include/geos/operation/distance/FacetSequence.h (rev 0)
+++ trunk/include/geos/operation/distance/FacetSequence.h 2016-04-25 18:06:50 UTC (rev 4203)
@@ -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
Added: trunk/include/geos/operation/distance/FacetSequenceTreeBuilder.h
===================================================================
--- trunk/include/geos/operation/distance/FacetSequenceTreeBuilder.h (rev 0)
+++ trunk/include/geos/operation/distance/FacetSequenceTreeBuilder.h 2016-04-25 18:06:50 UTC (rev 4203)
@@ -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
Modified: trunk/include/geos/operation/distance/Makefile.am
===================================================================
--- trunk/include/geos/operation/distance/Makefile.am 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/include/geos/operation/distance/Makefile.am 2016-04-25 18:06:50 UTC (rev 4203)
@@ -11,4 +11,6 @@
ConnectedElementLocationFilter.h \
ConnectedElementPointFilter.h \
DistanceOp.h \
+ FacetSequence.h \
+ FacetSequenceTreeBuilder.h \
GeometryLocation.h
Modified: trunk/include/geos/precision/Makefile.am
===================================================================
--- trunk/include/geos/precision/Makefile.am 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/include/geos/precision/Makefile.am 2016-04-25 18:06:50 UTC (rev 4203)
@@ -13,5 +13,6 @@
CommonBitsRemover.h \
EnhancedPrecisionOp.h \
GeometryPrecisionReducer.h \
+ MinimumClearance.h \
PrecisionReducerCoordinateOperation.h \
SimpleGeometryPrecisionReducer.h
Added: trunk/include/geos/precision/MinimumClearance.h
===================================================================
--- trunk/include/geos/precision/MinimumClearance.h (rev 0)
+++ trunk/include/geos/precision/MinimumClearance.h 2016-04-25 18:06:50 UTC (rev 4203)
@@ -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
+
+
Modified: trunk/src/index/strtree/BoundablePair.cpp
===================================================================
--- trunk/src/index/strtree/BoundablePair.cpp 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/src/index/strtree/BoundablePair.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -35,10 +35,11 @@
return boundable2;
}
-double BoundablePair::distance() {
+double BoundablePair::distance() const {
// if items, compute exact distance
- if (isLeaves())
+ if (isLeaves()) {
return itemDistance->distance((ItemBoundable*) boundable1, (ItemBoundable*) boundable2);
+ }
// otherwise compute distance between bounds of boundables
const geom::Envelope* e1 = (const geom::Envelope*) boundable1->getBounds();
@@ -99,11 +100,9 @@
std::vector<Boundable*> *children = ((AbstractNode*) bndComposite)->getChildBoundables();
for(std::vector<Boundable*>::iterator it = children->begin(); it != children->end(); ++it) {
Boundable* child = *it;
- BoundablePair* bp = new BoundablePair(child, bndOther, itemDistance);
- if (bp->getDistance() < minDistance) {
- priQ.push(bp);
- } else {
- delete bp;
+ std::auto_ptr<BoundablePair> bp(new BoundablePair(child, bndOther, itemDistance));
+ if (minDistance == std::numeric_limits<double>::infinity() || bp->getDistance() < minDistance) {
+ priQ.push(bp.release());
}
}
}
Modified: trunk/src/index/strtree/STRtree.cpp
===================================================================
--- trunk/src/index/strtree/STRtree.cpp 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/src/index/strtree/STRtree.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -27,6 +27,7 @@
#include <algorithm> // std::sort
#include <iostream> // for debugging
#include <limits>
+#include <geos/util/GEOSException.h>
using namespace std;
using namespace geos::geom;
@@ -191,14 +192,19 @@
ItemBoundable bnd = ItemBoundable(env, (void*) item);
BoundablePair bp(getRoot(), &bnd, itemDist);
- return nearestNeighbour(&bp);
+ return nearestNeighbour(&bp).first;
}
-const void* STRtree::nearestNeighbour(BoundablePair* initBndPair) {
- return nearestNeighbour(initBndPair, DoubleInfinity);
+std::pair<const void*, const void*> STRtree::nearestNeighbour(BoundablePair* initBndPair) {
+ return nearestNeighbour(initBndPair, std::numeric_limits<double>::infinity());
}
-const void* STRtree::nearestNeighbour(BoundablePair* initBndPair, double maxDistance) {
+std::pair<const void*, const void*> STRtree::nearestNeighbour(ItemDistance * itemDist) {
+ BoundablePair bp(this->getRoot(), this->getRoot(), itemDist);
+ return nearestNeighbour(&bp);
+}
+
+std::pair<const void*, const void*> STRtree::nearestNeighbour(BoundablePair* initBndPair, double maxDistance) {
double distanceLowerBound = maxDistance;
BoundablePair* minPair = NULL;
@@ -216,7 +222,7 @@
* So the current minDistance must be the true minimum,
* and we are done.
*/
- if (currentDistance >= distanceLowerBound)
+ if (minPair && currentDistance >= distanceLowerBound)
break;
priQ.pop();
@@ -246,16 +252,21 @@
/* Free any remaining BoundablePairs in the queue */
while(!priQ.empty()) {
- BoundablePair* bp = priQ.top();
+ BoundablePair* bndPair = priQ.top();
priQ.pop();
- delete bp;
+ if (bndPair != initBndPair)
+ delete bndPair;
}
- const void* retval = dynamic_cast<const ItemBoundable*>(minPair->getBoundable(0))->getItem();
+ if (!minPair)
+ throw util::GEOSException("Error computing nearest neighbor");
+
+ const void* item0 = dynamic_cast<const ItemBoundable*>(minPair->getBoundable(0))->getItem();
+ const void* item1 = dynamic_cast<const ItemBoundable*>(minPair->getBoundable(1))->getItem();
if (minPair != initBndPair)
delete minPair;
- return retval;
+ return std::pair<const void*, const void*>(item0, item1);
}
class STRAbstractNode: public AbstractNode{
Added: trunk/src/operation/distance/FacetSequence.cpp
===================================================================
--- trunk/src/operation/distance/FacetSequence.cpp (rev 0)
+++ trunk/src/operation/distance/FacetSequence.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -0,0 +1,117 @@
+/**********************************************************************
+ *
+ * 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)
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/CGAlgorithms.h>
+#include <geos/operation/distance/FacetSequence.h>
+
+using namespace geos::operation::distance;
+using namespace geos::algorithm;
+
+FacetSequence::FacetSequence(const CoordinateSequence* pts, size_t start, size_t end) :
+ pts(pts),
+ start(start),
+ end(end) {
+ computeEnvelope();
+}
+
+size_t FacetSequence::size() const {
+ return end - start;
+}
+
+bool FacetSequence::isPoint() const {
+ return end - start == 1;
+}
+
+double FacetSequence::distance(const FacetSequence & facetSeq) {
+ bool isPointThis = isPoint();
+ bool isPointOther = facetSeq.isPoint();
+
+ if (isPointThis && isPointOther) {
+ Coordinate pt = pts->getAt(start);
+ Coordinate seqPt = facetSeq.pts->getAt(start);
+ return pt.distance(seqPt);
+
+ } else if (isPointThis) {
+ Coordinate pt = pts->getAt(start);
+ return computePointLineDistance(pt, facetSeq);
+ } else if (isPointOther) {
+ Coordinate seqPt = facetSeq.pts->getAt(start);
+ return computePointLineDistance(seqPt, *this);
+ }
+
+ return computeLineLineDistance(facetSeq);
+}
+
+double FacetSequence::computePointLineDistance(const Coordinate & pt, const FacetSequence & facetSeq) const {
+ double minDistance = std::numeric_limits<double>::infinity();
+ double dist;
+ Coordinate q0;
+ Coordinate q1;
+
+ for (size_t i = facetSeq.start; i < facetSeq.end - 1; i++) {
+ facetSeq.pts->getAt(i, q0);
+ facetSeq.pts->getAt(i + 1, q1);
+ dist = CGAlgorithms::distancePointLine(pt, q0, q1);
+ if (dist == 0.0)
+ return dist;
+ if (dist < minDistance)
+ minDistance = dist;
+ }
+
+ return minDistance;
+}
+
+double FacetSequence::computeLineLineDistance(const FacetSequence & facetSeq) const {
+ double minDistance = std::numeric_limits<double>::infinity();
+ double dist;
+ Coordinate p0, p1, q0, q1;
+
+ for (size_t i = start; i < end - 1; i++) {
+ pts->getAt(i, p0);
+ pts->getAt(i + 1, p1);
+
+ for (size_t j = facetSeq.start; j < facetSeq.end - 1; j++) {
+ facetSeq.pts->getAt(j, q0);
+ facetSeq.pts->getAt(j + 1, q1);
+
+ dist = CGAlgorithms::distanceLineLine(p0, p1, q0, q1);
+ if (dist == 0.0)
+ return dist;
+ if (dist < minDistance)
+ minDistance = dist;
+ }
+ }
+
+ return minDistance;
+}
+
+void FacetSequence::computeEnvelope() {
+ env = Envelope();
+ for (size_t i = start; i < end; i++) {
+ env.expandToInclude(pts->getX(i), pts->getY(i));
+ }
+}
+
+const Envelope * FacetSequence::getEnvelope() const {
+ return &env;
+}
+
+const Coordinate * FacetSequence::getCoordinate(size_t index) const {
+ return &(pts->getAt(start + index));
+}
+
Added: trunk/src/operation/distance/FacetSequenceTreeBuilder.cpp
===================================================================
--- trunk/src/operation/distance/FacetSequenceTreeBuilder.cpp (rev 0)
+++ trunk/src/operation/distance/FacetSequenceTreeBuilder.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -0,0 +1,76 @@
+/**********************************************************************
+ *
+ * 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)
+ *
+ **********************************************************************/
+
+#include <geos/operation/distance/FacetSequenceTreeBuilder.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/Point.h>
+
+STRtree* FacetSequenceTreeBuilder::build(const Geometry* g) {
+ std::auto_ptr<STRtree> tree(new STRtree(STR_TREE_NODE_CAPACITY));
+ std::auto_ptr<std::vector<FacetSequence*> > sections(computeFacetSequences(g));
+ for (std::vector<FacetSequence*>::iterator it = sections->begin(); it != sections->end(); ++it) {
+ FacetSequence* section = *it;
+ tree->insert(section->getEnvelope(), section);
+ }
+
+ tree->build();
+ return tree.release();
+}
+
+std::vector<FacetSequence*> * FacetSequenceTreeBuilder::computeFacetSequences(const Geometry* g) {
+ std::auto_ptr<std::vector<FacetSequence*> > sections(new std::vector<FacetSequence*>());
+
+ class FacetSequenceAdder;
+ class FacetSequenceAdder : public geom::GeometryComponentFilter {
+ std::vector<FacetSequence*>* m_sections;
+
+ public :
+ FacetSequenceAdder(std::vector<FacetSequence*> * p_sections) :
+ m_sections(p_sections) {}
+ void filter_ro(const Geometry* geom) {
+ if (const LineString* ls = dynamic_cast<const LineString*>(geom)) {
+ const CoordinateSequence* seq = ls->getCoordinatesRO();
+ addFacetSequences(seq, *m_sections);
+ } else if (const Point* pt = dynamic_cast<const Point*>(geom)) {
+ const CoordinateSequence* seq = pt->getCoordinatesRO();
+ addFacetSequences(seq, *m_sections);
+ }
+ }
+ };
+
+ FacetSequenceAdder facetSequenceAdder(sections.get());
+ g->apply_ro(&facetSequenceAdder);
+
+ return sections.release();
+}
+
+void FacetSequenceTreeBuilder::addFacetSequences(const CoordinateSequence* pts, std::vector<FacetSequence*> & sections) {
+ size_t i = 0;
+ size_t size = pts->size();
+
+ while (i <= size - 1) {
+ size_t end = i + FACET_SEQUENCE_SIZE + 1;
+ // if only one point remains after this section, include it in this
+ // section
+ if (end >= size - 1)
+ end = size;
+ FacetSequence* sect = new FacetSequence(pts, i, end);
+ sections.push_back(sect);
+ i += FACET_SEQUENCE_SIZE;
+ }
+}
Modified: trunk/src/operation/distance/Makefile.am
===================================================================
--- trunk/src/operation/distance/Makefile.am 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/src/operation/distance/Makefile.am 2016-04-25 18:06:50 UTC (rev 4203)
@@ -11,6 +11,8 @@
ConnectedElementLocationFilter.cpp \
ConnectedElementPointFilter.cpp \
DistanceOp.cpp \
+ FacetSequence.cpp \
+ FacetSequenceTreeBuilder.cpp \
GeometryLocation.cpp
libopdistance_la_LIBADD =
Modified: trunk/src/precision/Makefile.am
===================================================================
--- trunk/src/precision/Makefile.am 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/src/precision/Makefile.am 2016-04-25 18:06:50 UTC (rev 4203)
@@ -13,6 +13,7 @@
CommonBitsRemover.cpp \
EnhancedPrecisionOp.cpp \
GeometryPrecisionReducer.cpp \
+ MinimumClearance.cpp \
PrecisionReducerCoordinateOperation.cpp \
SimpleGeometryPrecisionReducer.cpp
Added: trunk/src/precision/MinimumClearance.cpp
===================================================================
--- trunk/src/precision/MinimumClearance.cpp (rev 0)
+++ trunk/src/precision/MinimumClearance.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -0,0 +1,186 @@
+/**********************************************************************
+ *
+ * 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)
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/CGAlgorithms.h>
+#include <geos/precision/MinimumClearance.h>
+#include <geos/index/strtree/STRtree.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/CoordinateArraySequenceFactory.h>
+#include <geos/operation/distance/FacetSequenceTreeBuilder.h>
+#include <geos/geom/LineSegment.h>
+#include <geos/index/ItemVisitor.h>
+
+namespace geos {
+namespace precision {
+
+MinimumClearance::MinimumClearance(const Geometry* g) : inputGeom(g) {}
+
+double MinimumClearance::getDistance() {
+ compute();
+ return minClearance;
+}
+
+std::auto_ptr<LineString> MinimumClearance::getLine() {
+ compute();
+
+ // return empty line string if no min pts were found
+ if (minClearance == std::numeric_limits<double>::infinity())
+ return std::auto_ptr<LineString>(inputGeom->getFactory()->createLineString());
+
+ return std::auto_ptr<LineString>(inputGeom->getFactory()->createLineString(minClearancePts->clone()));
+}
+
+void MinimumClearance::compute() {
+ class MinClearanceDistance : public ItemDistance {
+ private:
+ double minDist;
+ std::vector<Coordinate> minPts;
+
+ void updatePts(const Coordinate & p, const Coordinate & seg0, const Coordinate & seg1) {
+ LineSegment seg(seg0, seg1);
+
+ minPts[0] = p;
+ seg.closestPoint(p, minPts[1]);
+ }
+
+ public:
+ MinClearanceDistance() :
+ minDist(std::numeric_limits<double>::infinity()),
+ minPts(std::vector<Coordinate>(2))
+ {}
+
+ const std::vector<Coordinate> * getCoordinates() {
+ return &minPts;
+ }
+
+ double distance(const ItemBoundable* b1, const ItemBoundable* b2) {
+ FacetSequence* fs1 = static_cast<FacetSequence*>(b1->getItem());
+ FacetSequence* fs2 = static_cast<FacetSequence*>(b2->getItem());
+
+ minDist = std::numeric_limits<double>::infinity();
+
+ return distance(fs1, fs2);
+ }
+
+ double distance(const FacetSequence* fs1, const FacetSequence* fs2) {
+ // Compute MinClearance distance metric
+
+ vertexDistance(fs1, fs2);
+ if (fs1->size() == 1 && fs2->size() == 1)
+ return minDist;
+ if (minDist <= 0.0)
+ return minDist;
+
+ segmentDistance(fs1, fs2);
+ if (minDist <= 0.0)
+ return minDist;
+
+ segmentDistance(fs2, fs1);
+ return minDist;
+ }
+
+ double vertexDistance(const FacetSequence* fs1, const FacetSequence* fs2) {
+ for (size_t i1 = 0; i1 < fs1->size(); i1++) {
+ for (size_t i2 = 0; i2 < fs2->size(); i2++) {
+ const Coordinate* p1 = fs1->getCoordinate(i1);
+ const Coordinate* p2 = fs2->getCoordinate(i2);
+ if (!p1->equals2D(*p2)) {
+ double d = p1->distance(*p2);
+ if (d < minDist) {
+ minDist = d;
+ minPts[0] = *p1;
+ minPts[1] = *p2;
+ if (d == 0.0)
+ return d;
+ }
+ }
+ }
+ }
+ return minDist;
+ }
+
+ double segmentDistance(const FacetSequence* fs1, const FacetSequence* fs2) {
+ for (size_t i1 = 0; i1 < fs1->size(); i1++) {
+ for (size_t i2 = 1; i2 < fs2->size(); i2++) {
+ const Coordinate* p = fs1->getCoordinate(i1);
+
+ const Coordinate* seg0 = fs2->getCoordinate(i2 - 1);
+ const Coordinate* seg1 = fs2->getCoordinate(i2);
+
+ if (! (p->equals2D(*seg0) || p->equals2D(*seg1))) {
+ double d = geos::algorithm::CGAlgorithms::distancePointLine(*p, *seg0, *seg1);
+ if (d < minDist) {
+ minDist = d;
+ updatePts(*p, *seg0, *seg1);
+ if (d == 0.0)
+ return d;
+ }
+ }
+ }
+ }
+ return minDist;
+ }
+ };
+
+ struct ItemDeleter : public index::ItemVisitor {
+ void visitItem(void * item) {
+ delete static_cast<FacetSequence*>(item);
+ }
+ };
+
+ struct ManagedResourceSTRtree {
+ STRtree* m_tree;
+
+ ManagedResourceSTRtree(STRtree * p_tree) : m_tree(p_tree) {}
+
+ ~ManagedResourceSTRtree() {
+ ItemDeleter id;
+ m_tree->iterate(id);
+
+ delete m_tree;
+ }
+ };
+
+ // already computed
+ if (minClearancePts.get() != NULL)
+ return;
+
+ // initialize to "No Distance Exists" state
+ minClearancePts = std::auto_ptr<CoordinateSequence>(inputGeom->getFactory()->getCoordinateSequenceFactory()->create(2, 2));
+ minClearance = std::numeric_limits<double>::infinity();
+
+ // handle empty geometries
+ if (inputGeom->isEmpty())
+ return;
+
+ ManagedResourceSTRtree tree(FacetSequenceTreeBuilder::build(inputGeom));
+ MinClearanceDistance mcd;
+ std::pair<const void *, const void *> nearest = tree.m_tree->nearestNeighbour(&mcd);
+
+ minClearance = mcd.distance(
+ static_cast<const FacetSequence *>(nearest.first),
+ static_cast<const FacetSequence *>(nearest.second));
+
+ const std::vector<Coordinate>* minClearancePtsVec = mcd.getCoordinates();
+ minClearancePts->setAt((*minClearancePtsVec)[0], 0);
+ minClearancePts->setAt((*minClearancePtsVec)[1], 1);
+}
+
+
+}
+}
Modified: trunk/tests/unit/Makefile.am
===================================================================
--- trunk/tests/unit/Makefile.am 2016-04-24 17:27:57 UTC (rev 4202)
+++ trunk/tests/unit/Makefile.am 2016-04-25 18:06:50 UTC (rev 4203)
@@ -132,6 +132,7 @@
capi/GEOSInterruptTest.cpp \
capi/GEOSIntersectsTest.cpp \
capi/GEOSIntersectionTest.cpp \
+ capi/GEOSMinimumClearanceTest.cpp \
capi/GEOSMinimumRectangleTest.cpp \
capi/GEOSMinimumWidthTest.cpp \
capi/GEOSNearestPointsTest.cpp \
Added: trunk/tests/unit/capi/GEOSMinimumClearanceTest.cpp
===================================================================
--- trunk/tests/unit/capi/GEOSMinimumClearanceTest.cpp (rev 0)
+++ trunk/tests/unit/capi/GEOSMinimumClearanceTest.cpp 2016-04-25 18:06:50 UTC (rev 4203)
@@ -0,0 +1,132 @@
+//
+// Test Suite for C-API GEOSMinimumClearance
+#include <tut.hpp>
+// geos
+#include <geos_c.h>
+// std
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+
+namespace tut
+{
+ //
+ // Test Group
+ //
+
+ // Common data used in test cases.
+ struct test_capigeosminimumclearance_data
+ {
+ static void notice(const char *fmt, ...)
+ {
+ std::fprintf( stdout, "NOTICE: ");
+
+ va_list ap;
+ va_start(ap, fmt);
+ std::vfprintf(stdout, fmt, ap);
+ va_end(ap);
+
+ std::fprintf(stdout, "\n");
+ }
+
+ test_capigeosminimumclearance_data()
+ {
+ initGEOS(notice, notice);
+ }
+
+ ~test_capigeosminimumclearance_data()
+ {
+ finishGEOS();
+ }
+
+ void testClearance(const std::string & wkx_input,
+ const std::string & wkx_expected,
+ double clearance) {
+
+ GEOSGeometry* input;
+ GEOSGeometry* expected_result;
+ if (wkx_input[0] == '0')
+ input = GEOSGeomFromHEX_buf((const unsigned char*) wkx_input.c_str(), wkx_input.length());
+ else
+ input = GEOSGeomFromWKT(wkx_input.c_str());
+
+ if (wkx_expected[0] == '0')
+ expected_result = GEOSGeomFromHEX_buf((const unsigned char*) wkx_expected.c_str(), wkx_expected.length());
+ else
+ expected_result = GEOSGeomFromWKT(wkx_expected.c_str());
+
+ double d;
+ int error = GEOSMinimumClearance(input, &d);
+
+ ensure(!error);
+ if (clearance == std::numeric_limits<double>::infinity()) {
+ ensure_equals(d, clearance);
+ } else {
+ ensure_equals("clearance", d, clearance, 1e-12);
+ }
+
+ GEOSGeometry* result = GEOSMinimumClearanceLine(input);
+ ensure(result != NULL);
+ ensure_equals(1, GEOSEquals(result, expected_result));
+
+ GEOSGeom_destroy(input);
+ GEOSGeom_destroy(expected_result);
+ GEOSGeom_destroy(result);
+ }
+
+ };
+
+ typedef test_group<test_capigeosminimumclearance_data> group;
+ typedef group::object object;
+
+ group test_capigeosminimumclearance_group("capi::GEOSMinimumClearance");
+
+ //
+ // Test Cases
+ //
+ template<>
+ template<>
+ void object::test<1>()
+ {
+ testClearance("MULTIPOINT ((100 100), (100 100))",
+ "LINESTRING EMPTY",
+ std::numeric_limits<double>::infinity());
+ }
+
+ template<>
+ template<>
+ void object::test<2>()
+ {
+ testClearance("MULTIPOINT ((100 100), (10 100), (30 100))",
+ "LINESTRING (30 100, 10 100)",
+ 20);
+ }
+
+ template<>
+ template<>
+ void object::test<3>()
+ {
+ testClearance("POLYGON ((100 100, 300 100, 200 200, 100 100))",
+ "LINESTRING (200 200, 200 100)",
+ 100);
+ }
+
+ template<>
+ template<>
+ void object::test<4>()
+ {
+ testClearance("0106000000010000000103000000010000001a00000035d42824992d5cc01b834e081dca404073b9c150872d5cc03465a71fd4c940400ec00644882d5cc03b8a73d4d1c94040376dc669882d5cc0bf9cd9aed0c940401363997e892d5cc002f4fbfecdc94040ca4e3fa88b2d5cc0a487a1d5c9c940408f1ce90c8c2d5cc0698995d1c8c94040fab836548c2d5cc0bd175fb4c7c940409f1f46088f2d5cc0962023a0c2c940407b15191d902d5cc068041bd7bfc940400397c79a912d5cc0287d21e4bcc940403201bf46922d5cc065e3c116bbc940409d9d0c8e922d5cc0060fd3beb9c940400ef7915b932d5cc09012bbb6b7c940404fe61f7d932d5cc0e4a08499b6c94040fc71fbe5932d5cc0ea9106b7b5c94040eaec6470942d5cc0c2323674b3c94040601dc70f952d5cc043588d25acc94040aea06989952d5cc03ecf9f36aac94040307f85cc952d5cc0e5eb32fca7c94040dd0a6135962d5cc01b615111a7c9404048a7ae7c962d5cc00a2aaa7ea5c94040f4328ae5962d5cc05eb87361a4c94040c49448a2972d5cc04d81cccea2c940407c80eecb992d5cc06745d4449fc9404035d42824992d5cc01b834e081dca4040",
+ "LINESTRING (-112.712119 33.575919, -112.712127 33.575885)",
+ 3.49284983912134e-05);
+ }
+
+ template<>
+ template<>
+ void object::test<5>()
+ {
+ testClearance("POLYGON EMPTY",
+ "LINESTRING EMPTY",
+ std::numeric_limits<double>::infinity());
+ }
+
+} // namespace tut
More information about the geos-commits
mailing list