From bb5839a2a396d296fc9b8021162b4ca98a7369bc Mon Sep 17 00:00:00 2001 From: dbaston Date: Fri, 22 Apr 2016 11:09:01 -0400 Subject: [PATCH 1/2] Implement GEOSMinimumClearance and GEOSMinimumClearanceLine --- NEWS | 3 + capi/geos_c.cpp | 12 ++ capi/geos_c.h.in | 33 ++++ capi/geos_ts_c.cpp | 67 +++++++ include/geos/index/strtree/BoundablePair.h | 2 +- include/geos/index/strtree/STRtree.h | 5 +- .../geos/operation/distance/FacetSequence.h | 66 +++++++ .../distance/FacetSequenceTreeBuilder.h | 52 +++++ include/geos/operation/distance/Makefile.am | 2 + include/geos/precision/Makefile.am | 1 + include/geos/precision/MinimumClearance.h | 60 ++++++ src/index/strtree/BoundablePair.cpp | 13 +- src/index/strtree/STRtree.cpp | 29 ++- src/operation/distance/FacetSequence.cpp | 117 +++++++++++ .../distance/FacetSequenceTreeBuilder.cpp | 76 +++++++ src/operation/distance/Makefile.am | 2 + src/precision/Makefile.am | 1 + src/precision/MinimumClearance.cpp | 185 ++++++++++++++++++ tests/unit/Makefile.am | 1 + tests/unit/capi/GEOSMinimumClearanceTest.cpp | 123 ++++++++++++ 20 files changed, 831 insertions(+), 19 deletions(-) create mode 100644 include/geos/operation/distance/FacetSequence.h create mode 100644 include/geos/operation/distance/FacetSequenceTreeBuilder.h create mode 100644 include/geos/precision/MinimumClearance.h create mode 100644 src/operation/distance/FacetSequence.cpp create mode 100644 src/operation/distance/FacetSequenceTreeBuilder.cpp create mode 100644 src/precision/MinimumClearance.cpp create mode 100644 tests/unit/capi/GEOSMinimumClearanceTest.cpp diff --git a/NEWS b/NEWS index 1a6b0e84dc..51206be65f 100644 --- a/NEWS +++ b/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: diff --git a/capi/geos_c.cpp b/capi/geos_c.cpp index 2c32635bec..fec3674940 100644 --- a/capi/geos_c.cpp +++ b/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) { diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in index bc12865351..4e1b90b219 100644 --- a/capi/geos_c.h.in +++ b/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); diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp index ef67bfe290..9f6fa4eb7e 100644 --- a/capi/geos_ts_c.cpp +++ b/capi/geos_ts_c.cpp @@ -116,6 +116,7 @@ using namespace std; #undef VERBOSE_EXCEPTIONS #include +#include // 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(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(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) diff --git a/include/geos/index/strtree/BoundablePair.h b/include/geos/index/strtree/BoundablePair.h index e78038349e..0e13823ebf 100644 --- a/include/geos/index/strtree/BoundablePair.h +++ b/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 diff --git a/include/geos/index/strtree/STRtree.h b/include/geos/index/strtree/STRtree.h index 05f77df1ff..f06f97acb1 100644 --- a/include/geos/index/strtree/STRtree.h +++ b/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 nearestNeighbour(BoundablePair* initBndPair); + std::pair nearestNeighbour(ItemDistance* itemDist); + std::pair nearestNeighbour(BoundablePair* initBndPair, double maxDistance); bool remove(const geom::Envelope *itemEnv, void* item) { return AbstractSTRtree::remove(itemEnv, item); diff --git a/include/geos/operation/distance/FacetSequence.h b/include/geos/operation/distance/FacetSequence.h new file mode 100644 index 0000000000..e329c5a2b9 --- /dev/null +++ b/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 +#include +#include + +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 diff --git a/include/geos/operation/distance/FacetSequenceTreeBuilder.h b/include/geos/operation/distance/FacetSequenceTreeBuilder.h new file mode 100644 index 0000000000..0f1dad94a8 --- /dev/null +++ b/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 +#include +#include +#include + +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 & sections); + static std::vector * computeFacetSequences(const Geometry* g); + + public: + static STRtree* build(const Geometry* g); + }; + } + } +} + +#endif //GEOS_FACETSEQUENCETREEBUILDER_H diff --git a/include/geos/operation/distance/Makefile.am b/include/geos/operation/distance/Makefile.am index 18cb44385f..a93fce92b6 100644 --- a/include/geos/operation/distance/Makefile.am +++ b/include/geos/operation/distance/Makefile.am @@ -11,4 +11,6 @@ geos_HEADERS = \ ConnectedElementLocationFilter.h \ ConnectedElementPointFilter.h \ DistanceOp.h \ + FacetSequence.h \ + FacetSequenceTreeBuilder.h \ GeometryLocation.h diff --git a/include/geos/precision/Makefile.am b/include/geos/precision/Makefile.am index 5c76288327..6f886752c2 100644 --- a/include/geos/precision/Makefile.am +++ b/include/geos/precision/Makefile.am @@ -13,5 +13,6 @@ geos_HEADERS = \ CommonBitsRemover.h \ EnhancedPrecisionOp.h \ GeometryPrecisionReducer.h \ + MinimumClearance.h \ PrecisionReducerCoordinateOperation.h \ SimpleGeometryPrecisionReducer.h diff --git a/include/geos/precision/MinimumClearance.h b/include/geos/precision/MinimumClearance.h new file mode 100644 index 0000000000..b05b607a43 --- /dev/null +++ b/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 +#include +#include + +namespace geos { +namespace precision { +class MinimumClearance { + private: + const geom::Geometry* inputGeom; + double minClearance; + std::auto_ptr minClearancePts; + + void compute(); + public: + MinimumClearance(const geom::Geometry* g); + + /** + * Gets the Minimum Clearance distance. + * + * @return the value of the minimum clearance distance + * or DBL_MAX 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 LINESTRING EMPTY if no Minimum Clearance distance exists + */ + std::auto_ptr getLine(); +}; +} +} + +#endif + + diff --git a/src/index/strtree/BoundablePair.cpp b/src/index/strtree/BoundablePair.cpp index e346007025..5db6849a2e 100644 --- a/src/index/strtree/BoundablePair.cpp +++ b/src/index/strtree/BoundablePair.cpp @@ -35,10 +35,11 @@ const Boundable* BoundablePair::getBoundable(int i) const { 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 @@ void BoundablePair::expand(const Boundable* bndComposite, const Boundable* bndOt std::vector *children = ((AbstractNode*) bndComposite)->getChildBoundables(); for(std::vector::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 bp(new BoundablePair(child, bndOther, itemDistance)); + if (minDistance == std::numeric_limits::infinity() || bp->getDistance() < minDistance) { + priQ.push(bp.release()); } } } diff --git a/src/index/strtree/STRtree.cpp b/src/index/strtree/STRtree.cpp index 1d5083452a..f5cd9505fe 100644 --- a/src/index/strtree/STRtree.cpp +++ b/src/index/strtree/STRtree.cpp @@ -27,6 +27,7 @@ #include // std::sort #include // for debugging #include +#include using namespace std; using namespace geos::geom; @@ -191,14 +192,19 @@ const void* STRtree::nearestNeighbour(const Envelope* env, const void* item, Ite ItemBoundable bnd = ItemBoundable(env, (void*) item); BoundablePair bp(getRoot(), &bnd, itemDist); - return nearestNeighbour(&bp); + return nearestNeighbour(&bp).first; +} + +std::pair STRtree::nearestNeighbour(BoundablePair* initBndPair) { + return nearestNeighbour(initBndPair, std::numeric_limits::infinity()); } -const void* STRtree::nearestNeighbour(BoundablePair* initBndPair) { - return nearestNeighbour(initBndPair, DoubleInfinity); +std::pair STRtree::nearestNeighbour(ItemDistance * itemDist) { + BoundablePair bp(this->getRoot(), this->getRoot(), itemDist); + return nearestNeighbour(&bp); } -const void* STRtree::nearestNeighbour(BoundablePair* initBndPair, double maxDistance) { +std::pair STRtree::nearestNeighbour(BoundablePair* initBndPair, double maxDistance) { double distanceLowerBound = maxDistance; BoundablePair* minPair = NULL; @@ -216,7 +222,7 @@ const void* STRtree::nearestNeighbour(BoundablePair* initBndPair, double maxDist * 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 @@ const void* STRtree::nearestNeighbour(BoundablePair* initBndPair, double maxDist /* 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(minPair->getBoundable(0))->getItem(); + if (!minPair) + throw util::GEOSException("Error computing nearest neighbor"); + + const void* item0 = dynamic_cast(minPair->getBoundable(0))->getItem(); + const void* item1 = dynamic_cast(minPair->getBoundable(1))->getItem(); if (minPair != initBndPair) delete minPair; - return retval; + return std::pair(item0, item1); } class STRAbstractNode: public AbstractNode{ diff --git a/src/operation/distance/FacetSequence.cpp b/src/operation/distance/FacetSequence.cpp new file mode 100644 index 0000000000..f285d53eda --- /dev/null +++ b/src/operation/distance/FacetSequence.cpp @@ -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 +#include + +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::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::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)); +} + diff --git a/src/operation/distance/FacetSequenceTreeBuilder.cpp b/src/operation/distance/FacetSequenceTreeBuilder.cpp new file mode 100644 index 0000000000..b2990f129a --- /dev/null +++ b/src/operation/distance/FacetSequenceTreeBuilder.cpp @@ -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 +#include +#include + +STRtree* FacetSequenceTreeBuilder::build(const Geometry* g) { + std::auto_ptr tree(new STRtree(STR_TREE_NODE_CAPACITY)); + std::auto_ptr > sections(computeFacetSequences(g)); + for (std::vector::iterator it = sections->begin(); it != sections->end(); ++it) { + FacetSequence* section = *it; + tree->insert(section->getEnvelope(), section); + } + + tree->build(); + return tree.release(); +} + +std::vector * FacetSequenceTreeBuilder::computeFacetSequences(const Geometry* g) { + std::auto_ptr > sections(new std::vector()); + + class FacetSequenceAdder; + class FacetSequenceAdder : public geom::GeometryComponentFilter { + std::vector* m_sections; + + public : + FacetSequenceAdder(std::vector * p_sections) : + m_sections(p_sections) {} + void filter_ro(const Geometry* geom) { + if (const LineString* ls = dynamic_cast(geom)) { + const CoordinateSequence* seq = ls->getCoordinatesRO(); + addFacetSequences(seq, *m_sections); + } else if (const Point* pt = dynamic_cast(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 & 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; + } +} diff --git a/src/operation/distance/Makefile.am b/src/operation/distance/Makefile.am index ec9a6af56c..0321697082 100644 --- a/src/operation/distance/Makefile.am +++ b/src/operation/distance/Makefile.am @@ -11,6 +11,8 @@ libopdistance_la_SOURCES = \ ConnectedElementLocationFilter.cpp \ ConnectedElementPointFilter.cpp \ DistanceOp.cpp \ + FacetSequence.cpp \ + FacetSequenceTreeBuilder.cpp \ GeometryLocation.cpp libopdistance_la_LIBADD = diff --git a/src/precision/Makefile.am b/src/precision/Makefile.am index b9fb2b76b5..152936f5ca 100644 --- a/src/precision/Makefile.am +++ b/src/precision/Makefile.am @@ -13,6 +13,7 @@ libprecision_la_SOURCES = \ CommonBitsRemover.cpp \ EnhancedPrecisionOp.cpp \ GeometryPrecisionReducer.cpp \ + MinimumClearance.cpp \ PrecisionReducerCoordinateOperation.cpp \ SimpleGeometryPrecisionReducer.cpp diff --git a/src/precision/MinimumClearance.cpp b/src/precision/MinimumClearance.cpp new file mode 100644 index 0000000000..9e07d667e4 --- /dev/null +++ b/src/precision/MinimumClearance.cpp @@ -0,0 +1,185 @@ +/********************************************************************** + * + * 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 +#include +#include +#include +#include +#include +#include +#include + +namespace geos { +namespace precision { + +MinimumClearance::MinimumClearance(const Geometry* g) : inputGeom(g) {} + +double MinimumClearance::getDistance() { + compute(); + return minClearance; +} + +std::auto_ptr MinimumClearance::getLine() { + compute(); + + // return empty line string if no min pts were found + if (minClearance == std::numeric_limits::infinity()) + return std::auto_ptr(inputGeom->getFactory()->createLineString()); + + return std::auto_ptr(inputGeom->getFactory()->createLineString(minClearancePts->clone())); +} + +void MinimumClearance::compute() { + class MinClearanceDistance : public ItemDistance { + private: + double minDist; + std::vector 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::infinity()), + minPts(std::vector(2)) + {} + + const std::vector * getCoordinates() { + return &minPts; + } + + double distance(const ItemBoundable* b1, const ItemBoundable* b2) { + FacetSequence* fs1 = static_cast(b1->getItem()); + FacetSequence* fs2 = static_cast(b2->getItem()); + + minDist = std::numeric_limits::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(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(inputGeom->getFactory()->getCoordinateSequenceFactory()->create(2, 2)); + + // handle empty geometries + if (inputGeom->isEmpty()) + return; + + ManagedResourceSTRtree tree(FacetSequenceTreeBuilder::build(inputGeom)); + MinClearanceDistance mcd; + std::pair nearest = tree.m_tree->nearestNeighbour(&mcd); + + minClearance = mcd.distance( + static_cast(nearest.first), + static_cast(nearest.second)); + + const std::vector* minClearancePtsVec = mcd.getCoordinates(); + minClearancePts->setAt((*minClearancePtsVec)[0], 0); + minClearancePts->setAt((*minClearancePtsVec)[1], 1); +} + + +} +} diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am index 983d1832a4..cf94a82101 100644 --- a/tests/unit/Makefile.am +++ b/tests/unit/Makefile.am @@ -132,6 +132,7 @@ geos_unit_SOURCES = \ capi/GEOSInterruptTest.cpp \ capi/GEOSIntersectsTest.cpp \ capi/GEOSIntersectionTest.cpp \ + capi/GEOSMinimumClearanceTest.cpp \ capi/GEOSMinimumRectangleTest.cpp \ capi/GEOSMinimumWidthTest.cpp \ capi/GEOSNearestPointsTest.cpp \ diff --git a/tests/unit/capi/GEOSMinimumClearanceTest.cpp b/tests/unit/capi/GEOSMinimumClearanceTest.cpp new file mode 100644 index 0000000000..bfc537e8fb --- /dev/null +++ b/tests/unit/capi/GEOSMinimumClearanceTest.cpp @@ -0,0 +1,123 @@ +// +// Test Suite for C-API GEOSMinimumClearance +#include +// geos +#include +// std +#include +#include +#include + +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::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 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::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); + } + +} // namespace tut From 0868446c832bb535059b4728b3385a377e3a0c11 Mon Sep 17 00:00:00 2001 From: dbaston Date: Sun, 24 Apr 2016 17:53:57 -0400 Subject: [PATCH 2/2] add tests for empty input geometry --- src/precision/MinimumClearance.cpp | 1 + tests/unit/capi/GEOSMinimumClearanceTest.cpp | 9 +++++++++ 2 files changed, 10 insertions(+) diff --git a/src/precision/MinimumClearance.cpp b/src/precision/MinimumClearance.cpp index 9e07d667e4..f7595d095f 100644 --- a/src/precision/MinimumClearance.cpp +++ b/src/precision/MinimumClearance.cpp @@ -162,6 +162,7 @@ void MinimumClearance::compute() { // initialize to "No Distance Exists" state minClearancePts = std::auto_ptr(inputGeom->getFactory()->getCoordinateSequenceFactory()->create(2, 2)); + minClearance = std::numeric_limits::infinity(); // handle empty geometries if (inputGeom->isEmpty()) diff --git a/tests/unit/capi/GEOSMinimumClearanceTest.cpp b/tests/unit/capi/GEOSMinimumClearanceTest.cpp index bfc537e8fb..932ecaf684 100644 --- a/tests/unit/capi/GEOSMinimumClearanceTest.cpp +++ b/tests/unit/capi/GEOSMinimumClearanceTest.cpp @@ -120,4 +120,13 @@ namespace tut 3.49284983912134e-05); } + template<> + template<> + void object::test<5>() + { + testClearance("POLYGON EMPTY", + "LINESTRING EMPTY", + std::numeric_limits::infinity()); + } + } // namespace tut