Port IndexedFacetDistance from JTS

NEWS

@@ -1,6 +1,9 @@
 Changes in 3.7.0
 YYYY-MM-DD
 
+- New things:
+  - CAPI: GEOSDistanceIndexed (#795, Dan Baston)
+
 Changes in 3.6.0
 2016-10-25
 

capi/geos_c.cpp

@@ -278,6 +278,12 @@ GEOSDistance(const Geometry *g1, const Geometry *g2, double *dist)
     return GEOSDistance_r( handle, g1, g2, dist );
 }
 
+int
+GEOSDistanceIndexed(const Geometry *g1, const Geometry *g2, double *dist)
+{
+    return GEOSDistanceIndexed_r( handle, g1, g2, dist );
+}
+
 int
 GEOSHausdorffDistance(const Geometry *g1, const Geometry *g2, double *dist)
 {

capi/geos_c.h.in

@@ -1100,6 +1100,9 @@ extern int GEOS_DLL GEOSLength_r(GEOSContextHandle_t handle,
 extern int GEOS_DLL GEOSDistance_r(GEOSContextHandle_t handle,
                                    const GEOSGeometry* g1,
                                    const GEOSGeometry* g2, double *dist);
+extern int GEOS_DLL GEOSDistanceIndexed_r(GEOSContextHandle_t handle,
+                                   const GEOSGeometry* g1,
+                                   const GEOSGeometry* g2, double *dist);
 extern int GEOS_DLL GEOSHausdorffDistance_r(GEOSContextHandle_t handle,
                                    const GEOSGeometry *g1,
                                    const GEOSGeometry *g2,
@@ -1963,6 +1966,8 @@ extern int GEOS_DLL GEOSArea(const GEOSGeometry* g, double *area);
 extern int GEOS_DLL GEOSLength(const GEOSGeometry* g, double *length);
 extern int GEOS_DLL GEOSDistance(const GEOSGeometry* g1, const GEOSGeometry* g2,
     double *dist);
+extern int GEOS_DLL GEOSDistanceIndexed(const GEOSGeometry* g1, const GEOSGeometry* g2,
+    double *dist);
 extern int GEOS_DLL GEOSHausdorffDistance(const GEOSGeometry *g1,
         const GEOSGeometry *g2, double *dist);
 extern int GEOS_DLL GEOSHausdorffDistanceDensify(const GEOSGeometry *g1,

capi/geos_ts_c.cpp

@@ -53,6 +53,7 @@
 #include <geos/operation/buffer/BufferOp.h>
 #include <geos/operation/buffer/BufferParameters.h>
 #include <geos/operation/distance/DistanceOp.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
 #include <geos/operation/linemerge/LineMerger.h>
 #include <geos/operation/overlay/OverlayOp.h>
 #include <geos/operation/overlay/snap/GeometrySnapper.h>
@@ -134,6 +135,7 @@ using geos::io::WKBWriter;
 using geos::operation::overlay::OverlayOp;
 using geos::operation::overlay::overlayOp;
 using geos::operation::geounion::CascadedPolygonUnion;
+using geos::operation::distance::IndexedFacetDistance;
 using geos::operation::buffer::BufferParameters;
 using geos::operation::buffer::BufferBuilder;
 using geos::precision::GeometryPrecisionReducer;
@@ -1157,6 +1159,40 @@ GEOSDistance_r(GEOSContextHandle_t extHandle, const Geometry *g1, const Geometry
     return 0;
 }
 
+int
+GEOSDistanceIndexed_r(GEOSContextHandle_t extHandle, const Geometry *g1, const Geometry *g2, double *dist)
+{
+    assert(0 != dist);
+
+    if ( 0 == extHandle )
+    {
+        return 0;
+    }
+
+    GEOSContextHandleInternal_t *handle = 0;
+    handle = reinterpret_cast<GEOSContextHandleInternal_t*>(extHandle);
+    if ( 0 == handle->initialized )
+    {
+        return 0;
+    }
+
+    try
+    {
+        *dist = IndexedFacetDistance::distance(g1, g2);
+        return 1;
+    }
+    catch (const std::exception &e)
+    {
+        handle->ERROR_MESSAGE("%s", e.what());
+    }
+    catch (...)
+    {
+        handle->ERROR_MESSAGE("Unknown exception thrown");
+    }
+
+    return 0;
+}
+
 int
 GEOSHausdorffDistance_r(GEOSContextHandle_t extHandle, const Geometry *g1, const Geometry *g2, double *dist)
 {
@@ -6074,6 +6110,7 @@ GEOSSTRtree_nearest_generic_r(GEOSContextHandle_t extHandle,
                               GEOSDistanceCallback distancefn,
                               void* userdata)
 {
+    using namespace geos::index::strtree;
 
     GEOSContextHandleInternal_t *handle = 0;
 

include/geos/index/strtree/BoundablePair.h

@@ -36,13 +36,10 @@
  * @author Martin Davis
  *
  */
-
-using namespace geos::index::strtree;
-
-
 namespace geos {
 namespace index {
 namespace strtree {
+
 class BoundablePair {
 	private:
 		const Boundable* boundable1;

include/geos/index/strtree/STRtree.h

@@ -143,10 +143,12 @@ using AbstractSTRtree::query;
 	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);
+	std::pair<const void*, const void*> nearestNeighbour(STRtree *tree, ItemDistance *itemDist);
 
 	bool remove(const geom::Envelope *itemEnv, void* item) {
 		return AbstractSTRtree::remove(itemEnv, item);
 	}
+
 };
 
 } // namespace geos::index::strtree

include/geos/operation/distance/FacetSequence.h

@@ -23,40 +23,38 @@
 #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 geom::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;
+                geom::Envelope env;
 
                 double computeLineLineDistance(const FacetSequence & facetSeq) const;
 
-                double computePointLineDistance(const Coordinate & pt, const FacetSequence & facetSeq) const;
+                double computePointLineDistance(const geom::Coordinate & pt, const FacetSequence & facetSeq) const;
 
                 void computeEnvelope();
 
             public:
-                const Envelope * getEnvelope() const;
+                const geom::Envelope * getEnvelope() const;
 
-                const Coordinate * getCoordinate(size_t index) const;
+                const geom::Coordinate * getCoordinate(size_t index) const;
 
                 size_t size() const;
 
                 bool isPoint() const;
 
-                double distance(const FacetSequence & facetSeq);
+                double distance(const FacetSequence & facetSeq) const;
 
-                FacetSequence(const CoordinateSequence *pts, size_t start, size_t end);
+                FacetSequence(const geom::CoordinateSequence *pts, size_t start, size_t end);
             };
 
         }

include/geos/operation/distance/FacetSequenceTreeBuilder.h

@@ -24,10 +24,6 @@
 #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 {
@@ -39,11 +35,11 @@ namespace geos {
                 // 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);
+                static void addFacetSequences(const geom::CoordinateSequence* pts, std::vector<FacetSequence*> & sections);
+                static std::vector<FacetSequence*> * computeFacetSequences(const geom::Geometry* g);
 
             public:
-                static STRtree* build(const Geometry* g);
+                static geos::index::strtree::STRtree* build(const geom::Geometry* g);
             };
         }
     }

include/geos/operation/distance/IndexedFacetDistance.h

@@ -0,0 +1,45 @@
+/**********************************************************************
+ *
+ * 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/IndexedFacetDistance.java (f6187ee2 JTS-1.14)
+ *
+ **********************************************************************/
+
+#ifndef GEOS_INDEXEDFACETDISTANCE_H
+#define GEOS_INDEXEDFACETDISTANCE_H
+
+#include <geos/operation/distance/FacetSequenceTreeBuilder.h>
+
+namespace geos {
+    namespace operation {
+        namespace distance {
+            class GEOS_DLL IndexedFacetDistance {
+            public:
+                IndexedFacetDistance(const geom::Geometry * g) :
+                        cachedTree(FacetSequenceTreeBuilder::build(g))
+                        {}
+
+                static double distance(const geom::Geometry * g1, const geom::Geometry * g2);
+
+                double getDistance(const geom::Geometry * g) const;
+
+            private:
+                std::auto_ptr<geos::index::strtree::STRtree> cachedTree;
+
+            };
+        }
+    }
+}
+
+#endif //GEOS_INDEXEDFACETDISTANCE_H

include/geos/operation/distance/Makefile.am

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

src/index/strtree/BoundablePair.cpp

@@ -21,6 +21,10 @@
 #include <geos/index/strtree/AbstractNode.h>
 #include <geos/util/IllegalArgumentException.h>
 
+namespace geos {
+namespace index {
+namespace strtree {
+
 BoundablePair::BoundablePair(const Boundable* boundable1, const Boundable* boundable2, ItemDistance* itemDistance) :
 	boundable1(boundable1),
 	boundable2(boundable2),
@@ -107,3 +111,7 @@ void BoundablePair::expand(const Boundable* bndComposite, const Boundable* bndOt
 	}
 }
 
+}
+}
+}
+

src/index/strtree/STRtree.cpp

@@ -204,6 +204,11 @@ std::pair<const void*, const void*> STRtree::nearestNeighbour(ItemDistance * ite
 	return nearestNeighbour(&bp);
 }
 
+std::pair<const void*, const void*> STRtree::nearestNeighbour(STRtree* tree, ItemDistance* itemDist) {
+	BoundablePair bp(getRoot(), tree->getRoot(), itemDist);
+	return nearestNeighbour(&bp);
+}
+
 std::pair<const void*, const void*> STRtree::nearestNeighbour(BoundablePair* initBndPair, double maxDistance) {
 	double distanceLowerBound = maxDistance;
 	BoundablePair* minPair = NULL;

src/operation/distance/FacetSequence.cpp

@@ -19,6 +19,7 @@
 #include <geos/algorithm/CGAlgorithms.h>
 #include <geos/operation/distance/FacetSequence.h>
 
+using namespace geos::geom;
 using namespace geos::operation::distance;
 using namespace geos::algorithm;
 
@@ -37,7 +38,7 @@ bool FacetSequence::isPoint() const {
     return end - start == 1;
 }
 
-double FacetSequence::distance(const FacetSequence & facetSeq) {
+double FacetSequence::distance(const FacetSequence & facetSeq) const {
     bool isPointThis = isPoint();
     bool isPointOther = facetSeq.isPoint();
 

src/operation/distance/FacetSequenceTreeBuilder.cpp

@@ -20,6 +20,13 @@
 #include <geos/geom/LineString.h>
 #include <geos/geom/Point.h>
 
+using namespace geos::geom;
+using namespace geos::index::strtree;
+
+namespace geos {
+namespace operation {
+namespace distance {
+
 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));
@@ -74,3 +81,8 @@ void FacetSequenceTreeBuilder::addFacetSequences(const CoordinateSequence* pts,
         i += FACET_SEQUENCE_SIZE;
     }
 }
+
+}
+}
+}
+

src/operation/distance/IndexedFacetDistance.cpp

@@ -0,0 +1,48 @@
+/**********************************************************************
+ *
+ * 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/IndexedFacetDistance.java (f6187ee2 JTS-1.14)
+ *
+ **********************************************************************/
+
+#include <geos/index/strtree/STRtree.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
+
+using namespace geos::geom;
+using namespace geos::index::strtree;
+
+namespace geos {
+    namespace operation {
+        namespace distance {
+            double IndexedFacetDistance::distance(const Geometry * g1, const Geometry * g2) {
+                IndexedFacetDistance ifd(g1);
+                return ifd.getDistance(g2);
+            }
+
+            double IndexedFacetDistance::getDistance(const Geometry * g) const {
+                struct : public ItemDistance {
+                    double distance(const ItemBoundable* item1, const ItemBoundable* item2) {
+                        return static_cast<const FacetSequence*>(item1->getItem())->distance(*static_cast<const FacetSequence*>(item2->getItem()));
+                    }
+                } itemDistance;
+
+                std::auto_ptr<STRtree> tree2(FacetSequenceTreeBuilder::build(g));
+
+                std::pair<const void*, const void*> obj = cachedTree->nearestNeighbour(tree2.get(), dynamic_cast<ItemDistance*>(&itemDistance));
+
+                return static_cast<const FacetSequence*>(obj.first)->distance(*static_cast<const FacetSequence*>(obj.second));
+            }
+        }
+    }
+}

src/operation/distance/Makefile.am

@@ -11,6 +11,7 @@ libopdistance_la_SOURCES = \
     ConnectedElementLocationFilter.cpp \
     ConnectedElementPointFilter.cpp \
     DistanceOp.cpp \
+    IndexedFacetDistance.cpp \
     FacetSequence.cpp \
     FacetSequenceTreeBuilder.cpp \
     GeometryLocation.cpp