[geos-commits] [SCM] GEOS branch master updated. 66908f9e268661a25cc38bd9d17a4bc566b516e9
git at osgeo.org
git at osgeo.org
Tue Apr 28 13:40:40 PDT 2020
This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GEOS".
The branch, master has been updated
via 66908f9e268661a25cc38bd9d17a4bc566b516e9 (commit)
via 0909042df5f8d4b9ae1832b3fd67713e206057fb (commit)
from 5ac945758d82a134c1c9bcde0c5c3f723b29f4c3 (commit)
Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.
- Log -----------------------------------------------------------------
commit 66908f9e268661a25cc38bd9d17a4bc566b516e9
Merge: 0909042 5ac9457
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Tue Apr 28 13:40:34 2020 -0700
Merge branch 'master' of https://git.osgeo.org/gitea/geos/geos
commit 0909042df5f8d4b9ae1832b3fd67713e206057fb
Author: Paul Ramsey <pramsey at cleverelephant.ca>
Date: Tue Apr 28 13:39:11 2020 -0700
Port MaximumInscribedCircle and LargestEmptyCircle
From https://github.com/locationtech/jts/pull/530
Closes #1029
diff --git a/.azure-pipelines.yml b/.azure-pipelines.yml
index 3c16549..2d4d4e9 100644
--- a/.azure-pipelines.yml
+++ b/.azure-pipelines.yml
@@ -158,11 +158,10 @@ stages:
matrix:
#Disabled until MSVC 16.4 is released
#(see https://developercommunity.visualstudio.com/content/problem/726778/this-snippet-compiles-totally-fine-on-godbolt-msvc.html?childToView=744771#comment-744771)
- #VS 2019 C++20 Strict:
- # TOOLSET: msvc-14.2
- # CXXSTD: 20
- # CXXFLAGS: -permissive-
- # VM_IMAGE: 'windows-2019'
+ VS 2019 C++20 Strict:
+ TOOLSET: msvc-14.2
+ CXXSTD: 20
+ VM_IMAGE: 'windows-2019'
VS 2017 C++17:
TOOLSET: msvc-14.1
CXXSTD: 17
@@ -171,10 +170,6 @@ stages:
TOOLSET: msvc-14.1
CXXSTD: 14 # default
VM_IMAGE: 'vs2017-win2016'
- VS 2015 C++14:
- TOOLSET: msvc-14.0
- CXXSTD: 14 # default
- VM_IMAGE: 'vs2015-win2012r2'
pool:
vmImage: $(VM_IMAGE)
steps:
@@ -202,7 +197,7 @@ stages:
- job: 'macOS'
pool:
- vmImage: 'macOS-10.13'
+ vmImage: 'macOS-10.14'
strategy:
matrix:
Xcode 10.1:
@@ -215,30 +210,6 @@ stages:
Xcode 9.4.1:
CXXSTD: 11, 14, 17
XCODE_APP: /Applications/Xcode_9.4.1.app
- Xcode 9.4:
- CXXSTD: 11, 14, 17
- XCODE_APP: /Applications/Xcode_9.4.app
- Xcode 9.3.1:
- CXXSTD: 11, 14, 17
- XCODE_APP: /Applications/Xcode_9.3.1.app
- Xcode 9.3:
- CXXSTD: 11, 14
- XCODE_APP: /Applications/Xcode_9.3.app
- Xcode 9.2:
- CXXSTD: 11, 14
- XCODE_APP: /Applications/Xcode_9.2.app
- Xcode 9.1:
- CXXSTD: 11
- XCODE_APP: /Applications/Xcode_9.1.app
- Xcode 9.0.1:
- CXXSTD: 11
- XCODE_APP: /Applications/Xcode_9.0.1.app
- Xcode 9.0:
- CXXSTD: 11
- XCODE_APP: /Applications/Xcode_9.app
- Xcode 8.3.3:
- CXXSTD: 11
- XCODE_APP: /Applications/Xcode_8.3.3.app
steps:
- script: |
set -e
diff --git a/.gitignore b/.gitignore
index 57712f8..697cb79 100644
--- a/.gitignore
+++ b/.gitignore
@@ -8,6 +8,7 @@ Makefile
.vimrc
.dirstamp
.tm_properties
+*~
*.o
*.orig
*.lo
diff --git a/NEWS b/NEWS
index db44ee1..f04c50f 100644
--- a/NEWS
+++ b/NEWS
@@ -2,13 +2,13 @@ Changes in 3.9.0
2020-xx-xx
- New things:
- -
+ - MaximumInscribedCircle and LargestEmptyCircle (JTS-530, Paul Ramsey)
- Improvements:
- Stack allocate segments in OverlapUnion (Paul Ramsey)
- Improve performance of GEOSisValid (Dan Baston)
- - Update geos-config tool for consistency
- and escape paths (https://git.osgeo.org/gitea/geos/geos/pulls/99)
+ - Update geos-config tool for consistency
+ and escape paths (https://git.osgeo.org/gitea/geos/geos/pulls/99)
changes mostly affect CMake MSVC builds (#1015, Mike Taves)
- Testing on Rasberry Pi 32-bit (berrie) (#1017, Bruce Rindahl, Regina Obe)
- Replace ttmath with JTS DD double-double implementation (Paul Ramsey)
diff --git a/capi/geos_c.cpp b/capi/geos_c.cpp
index acf5cb3..8e4547d 100644
--- a/capi/geos_c.cpp
+++ b/capi/geos_c.cpp
@@ -458,6 +458,18 @@ extern "C" {
}
Geometry*
+ GEOSMaximumInscribedCircle(const Geometry* g, double tolerance)
+ {
+ return GEOSMaximumInscribedCircle_r(handle, g, tolerance);
+ }
+
+ Geometry*
+ GEOSLargestEmptyCircle(const Geometry* g, double tolerance)
+ {
+ return GEOSLargestEmptyCircle_r(handle, g, tolerance);
+ }
+
+ Geometry*
GEOSMinimumWidth(const Geometry* g)
{
return GEOSMinimumWidth_r(handle, g);
diff --git a/capi/geos_c.h.in b/capi/geos_c.h.in
index 0de0be1..e411f79 100644
--- a/capi/geos_c.h.in
+++ b/capi/geos_c.h.in
@@ -561,6 +561,9 @@ extern GEOSGeometry GEOS_DLL *GEOSConvexHull_r(GEOSContextHandle_t handle,
extern GEOSGeometry GEOS_DLL *GEOSMinimumRotatedRectangle_r(GEOSContextHandle_t handle,
const GEOSGeometry* g);
+extern GEOSGeometry GEOS_DLL *GEOSMaximumInscribedCircle_r(GEOSContextHandle_t handle, const GEOSGeometry* g, double tolerance);
+extern GEOSGeometry GEOS_DLL *GEOSLargestEmptyCircle_r(GEOSContextHandle_t handle, const GEOSGeometry* g, double tolerance);
+
/* Returns a LINESTRING geometry which represents the minimum diameter of the geometry.
* The minimum diameter is defined to be the width of the smallest band that
* contains the geometry, where a band is a strip of the plane defined
@@ -1599,6 +1602,38 @@ extern GEOSGeometry GEOS_DLL *GEOSConvexHull(const GEOSGeometry* g);
*/
extern GEOSGeometry GEOS_DLL *GEOSMinimumRotatedRectangle(const GEOSGeometry* g);
+/* Constructs the Maximum Inscribed Circle for a polygonal geometry, up to a specified tolerance.
+ * The Maximum Inscribed Circle is determined by a point in the interior of the area
+ * which has the farthest distance from the area boundary, along with a boundary point at that distance.
+ * In the context of geography the center of the Maximum Inscribed Circle is known as the
+ * Pole of Inaccessibility. A cartographic use case is to determine a suitable point
+ * to place a map label within a polygon.
+ * The radius length of the Maximum Inscribed Circle is a measure of how "narrow" a polygon is. It is the
+ * distance at which the negative buffer becomes empty.
+ * The class supports polygons with holes and multipolygons.
+ * The implementation uses a successive-approximation technique over a grid of square cells covering the area geometry.
+ * The grid is refined using a branch-and-bound algorithm. Point containment and distance are computed in a performant
+ * way by using spatial indexes.
+ * Returns a two-point linestring, with one point at the center of the inscribed circle and the other
+ * on the boundary of the inscribed circle.
+*/
+extern GEOSGeometry GEOS_DLL *GEOSMaximumInscribedCircle(const GEOSGeometry* g, double tolerance);
+
+/* Constructs the Largest Empty Circle for a set of obstacle geometries, up to a
+ * specified tolerance. The obstacles are point and line geometries.
+ * The Largest Empty Circle is the largest circle which has its center in the convex hull of the
+ * obstacles (the boundary), and whose interior does not intersect with any obstacle.
+ * The circle center is the point in the interior of the boundary which has the farthest distance from
+ * the obstacles (up to tolerance). The circle is determined by the center point and a point lying on an
+ * obstacle indicating the circle radius.
+ * The implementation uses a successive-approximation technique over a grid of square cells covering the obstacles and boundary.
+ * The grid is refined using a branch-and-bound algorithm. Point containment and distance are computed in a performant
+ * way by using spatial indexes.
+ * Returns a two-point linestring, with one point at the center of the inscribed circle and the other
+ * on the boundary of the inscribed circle.
+ */
+extern GEOSGeometry GEOS_DLL *GEOSLargestEmptyCircle(const GEOSGeometry* g, double tolerance);
+
/* Returns a LINESTRING geometry which represents the minimum diameter of the geometry.
* The minimum diameter is defined to be the width of the smallest band that
* contains the geometry, where a band is a strip of the plane defined
diff --git a/capi/geos_ts_c.cpp b/capi/geos_ts_c.cpp
index d46a0af..a4b2692 100644
--- a/capi/geos_ts_c.cpp
+++ b/capi/geos_ts_c.cpp
@@ -48,6 +48,8 @@
#include <geos/algorithm/MinimumBoundingCircle.h>
#include <geos/algorithm/MinimumDiameter.h>
#include <geos/algorithm/Orientation.h>
+#include <geos/algorithm/construct/MaximumInscribedCircle.h>
+#include <geos/algorithm/construct/LargestEmptyCircle.h>
#include <geos/algorithm/distance/DiscreteHausdorffDistance.h>
#include <geos/algorithm/distance/DiscreteFrechetDistance.h>
#include <geos/simplify/DouglasPeuckerSimplifier.h>
@@ -1159,6 +1161,28 @@ extern "C" {
}
Geometry*
+ GEOSMaximumInscribedCircle_r(GEOSContextHandle_t extHandle, const Geometry* g, double tolerance)
+ {
+ return execute(extHandle, [&]() {
+ geos::algorithm::construct::MaximumInscribedCircle mic(g, tolerance);
+ auto g3 = mic.getRadiusLine();
+ g3->setSRID(g->getSRID());
+ return g3.release();
+ });
+ }
+
+ Geometry*
+ GEOSLargestEmptyCircle_r(GEOSContextHandle_t extHandle, const Geometry* g, double tolerance)
+ {
+ return execute(extHandle, [&]() {
+ geos::algorithm::construct::LargestEmptyCircle lec(g, tolerance);
+ auto g3 = lec.getRadiusLine();
+ g3->setSRID(g->getSRID());
+ return g3.release();
+ });
+ }
+
+ Geometry*
GEOSMinimumWidth_r(GEOSContextHandle_t extHandle, const Geometry* g)
{
return execute(extHandle, [&]() {
diff --git a/configure.ac b/configure.ac
index 1eade62..1c149aa 100644
--- a/configure.ac
+++ b/configure.ac
@@ -364,8 +364,9 @@ AC_OUTPUT([
macros/Makefile
src/Makefile
src/algorithm/Makefile
- src/algorithm/locate/Makefile
+ src/algorithm/construct/Makefile
src/algorithm/distance/Makefile
+ src/algorithm/locate/Makefile
src/geom/Makefile
src/geom/prep/Makefile
src/geom/util/Makefile
@@ -376,6 +377,7 @@ AC_OUTPUT([
include/geos/algorithm/Makefile
include/geos/algorithm/locate/Makefile
include/geos/algorithm/distance/Makefile
+ include/geos/algorithm/construct/Makefile
include/geos/geom/Makefile
include/geos/geom/prep/Makefile
include/geos/geom/util/Makefile
diff --git a/include/geos/algorithm/Makefile.am b/include/geos/algorithm/Makefile.am
index dbd5a61..67fc727 100644
--- a/include/geos/algorithm/Makefile.am
+++ b/include/geos/algorithm/Makefile.am
@@ -3,7 +3,8 @@
#
SUBDIRS = \
locate \
- distance
+ distance \
+ construct
EXTRA_DIST =
diff --git a/include/geos/algorithm/construct/LargestEmptyCircle.h b/include/geos/algorithm/construct/LargestEmptyCircle.h
new file mode 100644
index 0000000..fb8c8de
--- /dev/null
+++ b/include/geos/algorithm/construct/LargestEmptyCircle.h
@@ -0,0 +1,207 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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: algorithm/construct/LargestEmptyCircle.java
+ * https://github.com/locationtech/jts/commit/98274a7ea9b40651e9de6323dc10fb2cac17a245
+ *
+ **********************************************************************/
+
+#ifndef GEOS_ALGORITHM_CONSTRUCT_LARGESTCIRCLE_H
+#define GEOS_ALGORITHM_CONSTRUCT_LARGESTCIRCLE_H
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/Envelope.h>
+#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
+
+#include <memory>
+#include <queue>
+
+
+
+namespace geos {
+namespace geom {
+class Coordinate;
+class Envelope;
+class Geometry;
+class GeometryFactory;
+class LineString;
+class Point;
+}
+namespace operation {
+namespace distance {
+class IndexedFacetDistance;
+}
+}
+}
+
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+namespace construct { // geos::algorithm::construct
+
+/**
+ * Computes the Euclidean distance (L2 metric) from a Point to a Geometry.
+ *
+ * Also computes two points which are separated by the distance.
+ */
+class GEOS_DLL LargestEmptyCircle {
+
+public:
+
+ /**
+ * Creates a new instance of a Largest Empty Circle construction.
+ *
+ * @param obstacles a geometry representing the obstacles (points and lines)
+ * @param tolerance the distance tolerance for computing the circle center point
+ */
+ LargestEmptyCircle(const geom::Geometry* polygonal, double tolerance);
+ ~LargestEmptyCircle() = default;
+
+ /**
+ * Computes the center point of the Largest Empty Circle
+ * `within a set of obstacles, up to a given tolerance distance.
+ *
+ * @param obstacles a geometry representing the obstacles (points and lines)
+ * @param tolerance the distance tolerance for computing the center point
+ * @return the center point of the Largest Empty Circle
+ */
+ static std::unique_ptr<geom::Point> getCenter(const geom::Geometry* polygonal, double tolerance);
+
+ /**
+ * Computes a radius line of the Largest Empty Circle
+ * within a set of obstacles, up to a given distance tolerance.
+ *
+ * @param obstacles a geometry representing the obstacles (points and lines)
+ * @param tolerance the distance tolerance for computing the center point
+ * @return a line from the center of the circle to a point on the edge
+ */
+ static std::unique_ptr<geom::LineString> getRadiusLine(const geom::Geometry* polygonal, double tolerance);
+
+ std::unique_ptr<geom::Point> getCenter();
+ std::unique_ptr<geom::Point> getRadiusPoint();
+ std::unique_ptr<geom::LineString> getRadiusLine();
+
+
+private:
+
+ /* private members */
+ double tolerance;
+ const geom::Geometry* obstacles;
+ const geom::GeometryFactory* factory;
+ std::unique_ptr<geom::Geometry> boundary; // convexhull(obstacles)
+ operation::distance::IndexedFacetDistance obstacleDistance;
+ bool done;
+ std::unique_ptr<algorithm::locate::IndexedPointInAreaLocator> ptLocator;
+ std::unique_ptr<operation::distance::IndexedFacetDistance> boundaryDistance;
+ geom::Coordinate centerPt;
+ geom::Coordinate radiusPt;
+
+ /* private methods */
+ void setBoundary(const geom::Geometry* obstacles);
+
+ /**
+ * Computes the signed distance from a point to the constraints
+ * (obstacles and boundary).
+ * Points outside the boundary polygon are assigned a negative distance.
+ * Their containing cells will be last in the priority queue
+ * (but will still end up being tested since they may be refined).
+ *
+ * @param p the point to compute the distance for
+ * @return the signed distance to the constraints (negative indicates outside the boundary)
+ */
+ double distanceToConstraints(const geom::Coordinate& c);
+ double distanceToConstraints(double x, double y);
+ void compute();
+
+ /* private class */
+ class Cell {
+ private:
+ static constexpr double SQRT2 = 1.4142135623730951;
+ double x;
+ double y;
+ double hSize;
+ double distance;
+ double maxDist;
+
+ public:
+ Cell(double p_x, double p_y, double p_hSize, double p_distanceToConstraints)
+ : x(p_x)
+ , y(p_y)
+ , hSize(p_hSize)
+ , distance(p_distanceToConstraints)
+ , maxDist(p_distanceToConstraints + (p_hSize*SQRT2))
+ {};
+
+ geom::Envelope getEnvelope() const
+ {
+ geom::Envelope env(x-hSize, x+hSize, y-hSize, y+hSize);
+ return env;
+ }
+
+ bool isFullyOutside() const
+ {
+ return maxDist < 0.0;
+ }
+ bool isOutside() const
+ {
+ return distance < 0.0;
+ }
+ double getMaxDistance() const
+ {
+ return maxDist;
+ }
+ double getDistance() const
+ {
+ return distance;
+ }
+ double getHSize() const
+ {
+ return hSize;
+ }
+ double getX() const
+ {
+ return x;
+ }
+ double getY() const
+ {
+ return y;
+ }
+ bool operator< (const Cell& rhs) const
+ {
+ return maxDist < rhs.maxDist;
+ }
+ bool operator> (const Cell& rhs) const
+ {
+ return maxDist > rhs.maxDist;
+ }
+ bool operator==(const Cell& rhs) const
+ {
+ return maxDist == rhs.maxDist;
+ }
+ };
+
+ bool mayContainCircleCenter(const Cell& cell, const Cell& farthestCell);
+ void createInitialGrid(const geom::Envelope* env, std::priority_queue<Cell>& cellQueue);
+ Cell createCentroidCell(const geom::Geometry* geom);
+
+};
+
+
+} // geos::algorithm::construct
+} // geos::algorithm
+} // geos
+
+#endif // GEOS_ALGORITHM_CONSTRUCT_LARGESTCIRCLE_H
diff --git a/include/geos/algorithm/construct/Makefile.am b/include/geos/algorithm/construct/Makefile.am
new file mode 100644
index 0000000..5986e03
--- /dev/null
+++ b/include/geos/algorithm/construct/Makefile.am
@@ -0,0 +1,12 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/)
+#
+SUBDIRS =
+
+EXTRA_DIST =
+
+geosdir = $(includedir)/geos/algorithm/construct
+
+geos_HEADERS = \
+ LargestEmptyCircle.h \
+ MaximumInscribedCircle.h
diff --git a/include/geos/algorithm/construct/MaximumInscribedCircle.h b/include/geos/algorithm/construct/MaximumInscribedCircle.h
new file mode 100644
index 0000000..fef9ccf
--- /dev/null
+++ b/include/geos/algorithm/construct/MaximumInscribedCircle.h
@@ -0,0 +1,198 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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: algorithm/construct/MaximumInscribedCircle.java
+ * https://github.com/locationtech/jts/commit/98274a7ea9b40651e9de6323dc10fb2cac17a245
+ *
+ **********************************************************************/
+
+#ifndef GEOS_ALGORITHM_CONSTRUCT_MAXIMUMCIRCLE_H
+#define GEOS_ALGORITHM_CONSTRUCT_MAXIMUMCIRCLE_H
+
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/Point.h>
+#include <geos/geom/Envelope.h>
+#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
+
+#include <memory>
+#include <queue>
+
+
+
+namespace geos {
+namespace geom {
+class Coordinate;
+class Envelope;
+class Geometry;
+class GeometryFactory;
+class LineString;
+class Point;
+}
+}
+
+using geos::algorithm::locate::IndexedPointInAreaLocator;
+using geos::operation::distance::IndexedFacetDistance;
+
+namespace geos {
+namespace algorithm { // geos::algorithm
+namespace construct { // geos::algorithm::construct
+
+/**
+ * Computes the Euclidean distance (L2 metric) from a Point to a Geometry.
+ *
+ * Also computes two points which are separated by the distance.
+ */
+class GEOS_DLL MaximumInscribedCircle {
+
+public:
+
+ MaximumInscribedCircle(const geom::Geometry* polygonal, double tolerance);
+ ~MaximumInscribedCircle() = default;
+
+ /**
+ * Gets the center point of the maximum inscribed circle
+ * (up to the tolerance distance).
+ *
+ * @return the center point of the maximum inscribed circle
+ */
+ std::unique_ptr<geom::Point> getCenter();
+
+ /**
+ * Gets a point defining the radius of the Maximum Inscribed Circle.
+ * This is a point on the boundary which is
+ * nearest to the computed center of the Maximum Inscribed Circle.
+ * The line segment from the center to this point
+ * is a radius of the constructed circle, and this point
+ * lies on the boundary of the circle.
+ *
+ * @return a point defining the radius of the Maximum Inscribed Circle
+ */
+ std::unique_ptr<geom::Point> getRadiusPoint();
+
+ /**
+ * Gets a line representing a radius of the Largest Empty Circle.
+ *
+ * @return a line from the center of the circle to a point on the edge
+ */
+ std::unique_ptr<geom::LineString> getRadiusLine();
+
+ /**
+ * Computes the center point of the Maximum Inscribed Circle
+ * of a polygonal geometry, up to a given tolerance distance.
+ *
+ * @param polygonal a polygonal geometry
+ * @param tolerance the distance tolerance for computing the center point
+ * @return the center point of the maximum inscribed circle
+ */
+ static std::unique_ptr<geom::Point> getCenter(const geom::Geometry* polygonal, double tolerance);
+
+ /**
+ * Computes a radius line of the Maximum Inscribed Circle
+ * of a polygonal geometry, up to a given tolerance distance.
+ *
+ * @param polygonal a polygonal geometry
+ * @param tolerance the distance tolerance for computing the center point
+ * @return a line from the center to a point on the circle
+ */
+ static std::unique_ptr<geom::LineString> getRadiusLine(const geom::Geometry* polygonal, double tolerance);
+
+private:
+
+ /* private members */
+ const geom::Geometry* inputGeom;
+ std::unique_ptr<geom::Geometry> inputGeomBoundary;
+ double tolerance;
+ IndexedFacetDistance indexedDistance;
+ IndexedPointInAreaLocator ptLocator;
+ const geom::GeometryFactory* factory;
+ bool done;
+ geom::Coordinate centerPt;
+ geom::Coordinate radiusPt;
+
+ /* private methods */
+ double distanceToBoundary(const geom::Coordinate& c);
+ double distanceToBoundary(double x, double y);
+ void compute();
+
+ /* private class */
+ class Cell {
+ private:
+ static constexpr double SQRT2 = 1.4142135623730951;
+ double x;
+ double y;
+ double hSize;
+ double distance;
+ double maxDist;
+
+ public:
+ Cell(double p_x, double p_y, double p_hSize, double p_distanceToBoundary)
+ : x(p_x)
+ , y(p_y)
+ , hSize(p_hSize)
+ , distance(p_distanceToBoundary)
+ , maxDist(p_distanceToBoundary+(p_hSize*SQRT2))
+ {};
+
+ geom::Envelope getEnvelope() const
+ {
+ geom::Envelope env(x-hSize, x+hSize, y-hSize, y+hSize);
+ return env;
+ }
+
+ double getMaxDistance() const
+ {
+ return maxDist;
+ }
+ double getDistance() const
+ {
+ return distance;
+ }
+ double getHSize() const
+ {
+ return hSize;
+ }
+ double getX() const
+ {
+ return x;
+ }
+ double getY() const
+ {
+ return y;
+ }
+ bool operator< (const Cell& rhs) const
+ {
+ return maxDist < rhs.maxDist;
+ }
+ bool operator> (const Cell& rhs) const
+ {
+ return maxDist > rhs.maxDist;
+ }
+ bool operator==(const Cell& rhs) const
+ {
+ return maxDist == rhs.maxDist;
+ }
+ };
+
+ void createInitialGrid(const geom::Envelope* env, std::priority_queue<Cell>& cellQueue);
+ Cell createCentroidCell(const geom::Geometry* geom);
+
+};
+
+
+} // geos::algorithm::construct
+} // geos::algorithm
+} // geos
+
+#endif // GEOS_ALGORITHM_CONSTRUCT_MAXIMUMCIRCLE_H
diff --git a/include/geos/geom/Envelope.h b/include/geos/geom/Envelope.h
index 7f7da87..9ecca80 100644
--- a/include/geos/geom/Envelope.h
+++ b/include/geos/geom/Envelope.h
@@ -453,6 +453,28 @@ public:
*/
double distance(const Envelope* env) const;
+ /** \brief
+ * Computes the distance between one Coordinate and an Envelope
+ * defined by two other Coordinates. The order of the Coordinates
+ * used to define the envelope is not significant.
+ *
+ * @param c the coordinate to from which distance should be found
+ * @param p1 first coordinate defining an envelope
+ * @param p2 second coordinate defining an envelope.
+ */
+ static double distanceToCoordinate(const Coordinate & c, const Coordinate & p1, const Coordinate & p2);
+
+ /** \brief
+ * Computes the squared distance between one Coordinate and an Envelope
+ * defined by two other Coordinates. The order of the Coordinates
+ * used to define the envelope is not significant.
+ *
+ * @param c the coordinate to from which distance should be found
+ * @param p1 first coordinate defining an envelope
+ * @param p2 second coordinate defining an envelope.
+ */
+ static double distanceSquaredToCoordinate(const Coordinate & c, const Coordinate & p1, const Coordinate & p2);
+
size_t hashCode() const;
private:
diff --git a/include/geos/geom/Envelope.inl b/include/geos/geom/Envelope.inl
index 30599f0..7942cb6 100644
--- a/include/geos/geom/Envelope.inl
+++ b/include/geos/geom/Envelope.inl
@@ -19,7 +19,9 @@
#ifndef GEOS_GEOM_ENVELOPE_INL
#define GEOS_GEOM_ENVELOPE_INL
+#include <algorithm> // std::min, std::max
#include <cassert>
+#include <numeric> // std::signbit
#include <geos/geom/Coordinate.h>
#include <geos/geom/Envelope.h>
@@ -262,6 +264,25 @@ Envelope::setToNull()
maxy = -1;
}
+INLINE double
+Envelope::distanceToCoordinate(const Coordinate & c, const Coordinate & p0, const Coordinate & p1) {
+ return std::sqrt(distanceSquaredToCoordinate(c, p0, p1));
+}
+
+INLINE double
+Envelope::distanceSquaredToCoordinate(const Coordinate & c, const Coordinate & p0, const Coordinate & p1) {
+ double xa = c.x - p0.x;
+ double xb = c.x - p1.x;
+ double ya = c.y - p0.y;
+ double yb = c.y - p1.y;
+
+ // If sign of a and b are not the same, then Envelope spans c and distance is zero.
+ double dx = (std::signbit(xa) == std::signbit(xb)) * std::min(std::abs(xa), std::abs(xb));
+ double dy = (std::signbit(ya) == std::signbit(yb)) * std::min(std::abs(ya), std::abs(yb));
+
+ return dx*dx + dy*dy;
+}
+
} // namespace geos::geom
} // namespace geos
diff --git a/include/geos/index/strtree/Boundable.h b/include/geos/index/strtree/Boundable.h
index e8c34fb..290ee1f 100644
--- a/include/geos/index/strtree/Boundable.h
+++ b/include/geos/index/strtree/Boundable.h
@@ -40,8 +40,7 @@ public:
virtual const void* getBounds() const = 0;
virtual bool isLeaf() const = 0;
- virtual
- ~Boundable() {}
+ virtual ~Boundable() {}
};
diff --git a/src/algorithm/Makefile.am b/src/algorithm/Makefile.am
index 75dfb6c..a4a8180 100644
--- a/src/algorithm/Makefile.am
+++ b/src/algorithm/Makefile.am
@@ -3,7 +3,8 @@
#
SUBDIRS = \
locate \
- distance
+ distance \
+ construct
noinst_LTLIBRARIES = libalgorithm.la
@@ -36,4 +37,5 @@ libalgorithm_la_SOURCES = \
libalgorithm_la_LIBADD = \
locate/liblocation.la \
- distance/libdistance.la
+ distance/libdistance.la \
+ construct/libconstruct.la
diff --git a/src/algorithm/construct/LargestEmptyCircle.cpp b/src/algorithm/construct/LargestEmptyCircle.cpp
new file mode 100644
index 0000000..6506442
--- /dev/null
+++ b/src/algorithm/construct/LargestEmptyCircle.cpp
@@ -0,0 +1,272 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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: algorithm/construct/LargestEmptyCircle.java
+ * https://github.com/locationtech/jts/commit/98274a7ea9b40651e9de6323dc10fb2cac17a245
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/construct/LargestEmptyCircle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CoordinateSequenceFactory.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/MultiPolygon.h>
+#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
+
+#include <typeinfo> // for dynamic_cast
+#include <cassert>
+
+using namespace geos::geom;
+
+
+namespace geos {
+namespace algorithm { // geos.algorithm
+namespace construct { // geos.algorithm.construct
+
+
+
+LargestEmptyCircle::LargestEmptyCircle(const Geometry* p_obstacles, double p_tolerance)
+ : tolerance(p_tolerance)
+ , obstacles(p_obstacles)
+ , factory(p_obstacles->getFactory())
+ , boundary(p_obstacles->convexHull())
+ , obstacleDistance(p_obstacles)
+ , done(false)
+{
+ if (p_obstacles->isEmpty()) {
+ throw util::IllegalArgumentException("Empty obstacles geometry is not supported");
+ }
+
+ // if boundary does not enclose an area cannot create a ptLocator
+ if (boundary->getDimension() >= 2) {
+ ptLocator.reset(new algorithm::locate::IndexedPointInAreaLocator(*(boundary.get())));
+ boundaryDistance.reset(new operation::distance::IndexedFacetDistance(boundary.get()));
+ }
+}
+
+
+/* public static */
+std::unique_ptr<Point>
+LargestEmptyCircle::getCenter(const Geometry* polygonal, double tolerance)
+{
+ LargestEmptyCircle lec(polygonal, tolerance);
+ return lec.getCenter();
+}
+
+/* public static */
+std::unique_ptr<LineString>
+LargestEmptyCircle::getRadiusLine(const Geometry* polygonal, double tolerance)
+{
+ LargestEmptyCircle lec(polygonal, tolerance);
+ return lec.getRadiusLine();
+}
+
+/* public */
+std::unique_ptr<Point>
+LargestEmptyCircle::getCenter()
+{
+ compute();
+ return std::unique_ptr<Point>(factory->createPoint(centerPt));
+}
+
+/* public */
+std::unique_ptr<Point>
+LargestEmptyCircle::getRadiusPoint()
+{
+ compute();
+ return std::unique_ptr<Point>(factory->createPoint(radiusPt));
+}
+
+/* public */
+std::unique_ptr<LineString>
+LargestEmptyCircle::getRadiusLine()
+{
+ compute();
+
+ auto cl = factory->getCoordinateSequenceFactory()->create(2);
+ cl->setAt(centerPt, 0);
+ cl->setAt(radiusPt, 1);
+ return factory->createLineString(std::move(cl));
+}
+
+
+/* private */
+void
+LargestEmptyCircle::createInitialGrid(const Envelope* env, std::priority_queue<Cell>& cellQueue)
+{
+ double minX = env->getMinX();
+ double maxX = env->getMaxX();
+ double minY = env->getMinY();
+ double maxY = env->getMaxY();
+ double width = env->getWidth();
+ double height = env->getHeight();
+ double cellSize = std::min(width, height);
+ double hSize = cellSize / 2.0;
+
+ // compute initial grid of cells to cover area
+ for (double x = minX; x < maxX; x += cellSize) {
+ for (double y = minY; y < maxY; y += cellSize) {
+ cellQueue.emplace(x+hSize, y+hSize, hSize, distanceToConstraints(x+hSize, y+hSize));
+ }
+ }
+}
+
+/* private */
+bool
+LargestEmptyCircle::mayContainCircleCenter(const Cell& cell, const Cell& farthestCell)
+{
+ /**
+ * Every point in the cell lies outside the boundary,
+ * so they cannot be the center point
+ */
+ if (cell.isFullyOutside())
+ return false;
+
+ /**
+ * The cell is outside, but overlaps the boundary
+ * so it may contain a point which should be checked.
+ * This is only the case if the potential overlap distance
+ * is larger than the tolerance.
+ */
+ if (cell.isOutside()) {
+ bool isOverlapSignificant = cell.getMaxDistance() > tolerance;
+ return isOverlapSignificant;
+ }
+
+ /**
+ * Cell is inside the boundary. It may contain the center
+ * if the maximum possible distance is greater than the current distance
+ * (up to tolerance).
+ */
+ double potentialIncrease = cell.getMaxDistance() - farthestCell.getDistance();
+ return potentialIncrease > tolerance;
+}
+
+
+/* private */
+double
+LargestEmptyCircle::distanceToConstraints(const Coordinate& c)
+{
+ bool isOutside = ptLocator && (Location::EXTERIOR == ptLocator->locate(&c));
+ std::unique_ptr<Point> pt(factory->createPoint(c));
+ if (isOutside) {
+ double boundaryDist = boundaryDistance->distance(pt.get());
+ return -boundaryDist;
+
+ }
+ double dist = obstacleDistance.distance(pt.get());
+ return dist;
+}
+
+/* private */
+double
+LargestEmptyCircle::distanceToConstraints(double x, double y)
+{
+ Coordinate coord(x, y);
+ return distanceToConstraints(coord);
+}
+
+/* private */
+LargestEmptyCircle::Cell
+LargestEmptyCircle::createCentroidCell(const Geometry* geom)
+{
+ Coordinate c;
+ geom->getCentroid(c);
+ Cell cell(c.x, c.y, 0, distanceToConstraints(c));
+ return cell;
+}
+
+
+/* private */
+void
+LargestEmptyCircle::compute()
+{
+
+ // check if already computed
+ if (done) return;
+
+ // if ptLocator is not present then result is degenerate (represented as zero-radius circle)
+ if (!ptLocator) {
+ const Coordinate* pt = obstacles->getCoordinate();
+ centerPt = *pt;
+ radiusPt = *pt;
+ done = true;
+ return;
+ }
+
+ // Priority queue of cells, ordered by decreasing distance from constraints
+ std::priority_queue<Cell> cellQueue;
+ createInitialGrid(obstacles->getEnvelopeInternal(), cellQueue);
+
+ Cell farthestCell = createCentroidCell(obstacles);
+
+ /**
+ * Carry out the branch-and-bound search
+ * of the cell space
+ */
+ while (!cellQueue.empty()) {
+
+ // pick the most promising cell from the queue
+ Cell cell = cellQueue.top();
+ cellQueue.pop();
+
+ // update the center cell if the candidate is further from the constraints
+ if (cell.getDistance() > farthestCell.getDistance()) {
+ farthestCell = cell;
+ }
+
+ /**
+ * If this cell may contain a better approximation to the center
+ * of the empty circle, then refine it (partition into subcells
+ * which are added into the queue for further processing).
+ * Otherwise the cell is pruned (not investigated further),
+ * since no point in it can be further than the current farthest distance.
+ */
+ if (mayContainCircleCenter(cell, farthestCell)) {
+ // split the cell into four sub-cells
+ double h2 = cell.getHSize() / 2;
+ cellQueue.emplace(cell.getX()-h2, cell.getY()-h2, h2, distanceToConstraints(cell.getX()-h2, cell.getY()-h2));
+ cellQueue.emplace(cell.getX()+h2, cell.getY()-h2, h2, distanceToConstraints(cell.getX()+h2, cell.getY()-h2));
+ cellQueue.emplace(cell.getX()-h2, cell.getY()+h2, h2, distanceToConstraints(cell.getX()-h2, cell.getY()+h2));
+ cellQueue.emplace(cell.getX()+h2, cell.getY()+h2, h2, distanceToConstraints(cell.getX()+h2, cell.getY()+h2));
+ }
+ }
+
+ // the farthest cell is the best approximation to the MIC center
+ Cell centerCell = farthestCell;
+ centerPt.x = centerCell.getX();
+ centerPt.y = centerCell.getY();
+
+ // compute radius point
+ std::unique_ptr<Point> centerPoint(factory->createPoint(centerPt));
+ std::vector<geom::Coordinate> nearestPts = obstacleDistance.nearestPoints(centerPoint.get());
+ radiusPt = nearestPts[0];
+
+ // flag computation
+ done = true;
+}
+
+
+
+} // namespace geos.algorithm.construct
+} // namespace geos.algorithm
+} // namespace geos
+
+
diff --git a/src/algorithm/construct/Makefile.am b/src/algorithm/construct/Makefile.am
new file mode 100644
index 0000000..d9e47aa
--- /dev/null
+++ b/src/algorithm/construct/Makefile.am
@@ -0,0 +1,14 @@
+#
+# This file is part of project GEOS (http://trac.osgeo.org/geos/)
+#
+SUBDIRS =
+
+noinst_LTLIBRARIES = libconstruct.la
+
+AM_CPPFLAGS = -I$(top_srcdir)/include
+
+libconstruct_la_SOURCES = \
+ LargestEmptyCircle.cpp \
+ MaximumInscribedCircle.cpp
+
+libconstruct_la_LIBADD =
diff --git a/src/algorithm/construct/MaximumInscribedCircle.cpp b/src/algorithm/construct/MaximumInscribedCircle.cpp
new file mode 100644
index 0000000..b68c942
--- /dev/null
+++ b/src/algorithm/construct/MaximumInscribedCircle.cpp
@@ -0,0 +1,232 @@
+/**********************************************************************
+ *
+ * GEOS - Geometry Engine Open Source
+ * http://geos.osgeo.org
+ *
+ * Copyright (C) 2020 Paul Ramsey <pramsey at cleverelephant.ca>
+ *
+ * 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: algorithm/construct/MaximumInscribedCircle.java
+ * https://github.com/locationtech/jts/commit/98274a7ea9b40651e9de6323dc10fb2cac17a245
+ *
+ **********************************************************************/
+
+#include <geos/algorithm/construct/MaximumInscribedCircle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateSequence.h>
+#include <geos/geom/CoordinateSequenceFactory.h>
+#include <geos/geom/Envelope.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/Polygon.h>
+#include <geos/geom/MultiPolygon.h>
+#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
+#include <geos/operation/distance/IndexedFacetDistance.h>
+
+#include <typeinfo> // for dynamic_cast
+#include <cassert>
+
+using namespace geos::geom;
+
+
+namespace geos {
+namespace algorithm { // geos.algorithm
+namespace construct { // geos.algorithm.construct
+
+
+/* public */
+MaximumInscribedCircle::MaximumInscribedCircle(const Geometry* polygonal, double p_tolerance)
+ : inputGeom(polygonal)
+ , inputGeomBoundary(polygonal->getBoundary())
+ , tolerance(p_tolerance)
+ , indexedDistance(inputGeomBoundary.get())
+ , ptLocator(*polygonal)
+ , factory(polygonal->getFactory())
+ , done(false)
+{
+ if (!(typeid(*polygonal) == typeid(Polygon) ||
+ typeid(*polygonal) == typeid(MultiPolygon))) {
+ throw util::IllegalArgumentException("Input geometry must be a Polygon or MultiPolygon");
+ }
+
+ if (polygonal->isEmpty()) {
+ throw util::IllegalArgumentException("Empty input geometry is not supported");
+ }
+}
+
+
+/* public static */
+std::unique_ptr<Point>
+MaximumInscribedCircle::getCenter(const Geometry* polygonal, double tolerance)
+{
+ MaximumInscribedCircle mic(polygonal, tolerance);
+ return mic.getCenter();
+}
+
+/* public static */
+std::unique_ptr<LineString>
+MaximumInscribedCircle::getRadiusLine(const Geometry* polygonal, double tolerance)
+{
+ MaximumInscribedCircle mic(polygonal, tolerance);
+ return mic.getRadiusLine();
+}
+
+/* public */
+std::unique_ptr<Point>
+MaximumInscribedCircle::getCenter()
+{
+ compute();
+ auto pt = factory->createPoint(centerPt);
+ return std::unique_ptr<Point>(pt);
+}
+
+/* public */
+std::unique_ptr<Point>
+MaximumInscribedCircle::getRadiusPoint()
+{
+ compute();
+ auto pt = factory->createPoint(radiusPt);
+ return std::unique_ptr<Point>(pt);
+}
+
+/* public */
+std::unique_ptr<LineString>
+MaximumInscribedCircle::getRadiusLine()
+{
+ compute();
+
+ auto cl = factory->getCoordinateSequenceFactory()->create(2);
+ cl->setAt(centerPt, 0);
+ cl->setAt(radiusPt, 1);
+ return factory->createLineString(std::move(cl));
+}
+
+/* private */
+void
+MaximumInscribedCircle::createInitialGrid(const Envelope* env, std::priority_queue<Cell>& cellQueue)
+{
+ double minX = env->getMinX();
+ double maxX = env->getMaxX();
+ double minY = env->getMinY();
+ double maxY = env->getMaxY();
+ double width = env->getWidth();
+ double height = env->getHeight();
+ double cellSize = std::min(width, height);
+ double hSize = cellSize / 2.0;
+
+ // compute initial grid of cells to cover area
+ for (double x = minX; x < maxX; x += cellSize) {
+ for (double y = minY; y < maxY; y += cellSize) {
+ cellQueue.emplace(x+hSize, y+hSize, hSize, distanceToBoundary(x+hSize, y+hSize));
+ }
+ }
+}
+
+/* private */
+double
+MaximumInscribedCircle::distanceToBoundary(const Coordinate& c)
+{
+ std::unique_ptr<Point> pt(factory->createPoint(c));
+ double dist = indexedDistance.distance(pt.get());
+ // double dist = inputGeomBoundary->distance(pt.get());
+ bool isOutside = (Location::EXTERIOR == ptLocator.locate(&c));
+ if (isOutside) return -dist;
+ return dist;
+}
+
+/* private */
+double
+MaximumInscribedCircle::distanceToBoundary(double x, double y)
+{
+ Coordinate coord(x, y);
+ return distanceToBoundary(coord);
+}
+
+/* private */
+MaximumInscribedCircle::Cell
+MaximumInscribedCircle::createCentroidCell(const Geometry* geom)
+{
+ Coordinate c;
+ geom->getCentroid(c);
+ Cell cell(c.x, c.y, 0, distanceToBoundary(c));
+ return cell;
+}
+
+/* private */
+void
+MaximumInscribedCircle::compute()
+{
+
+ // check if already computed
+ if (done) return;
+
+ // Priority queue of cells, ordered by maximum distance from boundary
+ std::priority_queue<Cell> cellQueue;
+
+ createInitialGrid(inputGeom->getEnvelopeInternal(), cellQueue);
+
+ // use the area centroid as the initial candidate center point
+ Cell farthestCell = createCentroidCell(inputGeom);
+
+ /**
+ * Carry out the branch-and-bound search
+ * of the cell space
+ */
+ while (!cellQueue.empty()) {
+ // pick the most promising cell from the queue
+ Cell cell = cellQueue.top();
+ cellQueue.pop();
+
+ // std::cout << i << ": (" << cell.getX() << ", " << cell.getY() << ") " << cell.getHSize() << " dist = " << cell.getDistance() << std::endl;
+
+ // update the center cell if the candidate is further from the boundary
+ if (cell.getDistance() > farthestCell.getDistance()) {
+ farthestCell = cell;
+ }
+ /**
+ * Refine this cell if the potential distance improvement
+ * is greater than the required tolerance.
+ * Otherwise the cell is pruned (not investigated further),
+ * since no point in it is further than
+ * the current farthest distance.
+ */
+ double potentialIncrease = cell.getMaxDistance() - farthestCell.getDistance();
+ if (potentialIncrease > tolerance) {
+ // split the cell into four sub-cells
+ double h2 = cell.getHSize() / 2;
+ cellQueue.emplace(cell.getX()-h2, cell.getY()-h2, h2, distanceToBoundary(cell.getX()-h2, cell.getY()-h2));
+ cellQueue.emplace(cell.getX()+h2, cell.getY()-h2, h2, distanceToBoundary(cell.getX()+h2, cell.getY()-h2));
+ cellQueue.emplace(cell.getX()-h2, cell.getY()+h2, h2, distanceToBoundary(cell.getX()-h2, cell.getY()+h2));
+ cellQueue.emplace(cell.getX()+h2, cell.getY()+h2, h2, distanceToBoundary(cell.getX()+h2, cell.getY()+h2));
+ }
+ }
+ // std::cout << "number of iterations: " << i << std::endl;
+
+ // the farthest cell is the best approximation to the MIC center
+ Cell centerCell = farthestCell;
+ centerPt.x = centerCell.getX();
+ centerPt.y = centerCell.getY();
+
+ // compute radius point
+ std::unique_ptr<Point> centerPoint(factory->createPoint(centerPt));
+ std::vector<geom::Coordinate> nearestPts = indexedDistance.nearestPoints(centerPoint.get());
+ radiusPt = nearestPts[0];
+
+ // flag computation
+ done = true;
+}
+
+
+
+
+} // namespace geos.algorithm.construct
+} // namespace geos.algorithm
+} // namespace geos
+
diff --git a/src/index/strtree/BoundablePair.cpp b/src/index/strtree/BoundablePair.cpp
index 60409b7..aa7820c 100644
--- a/src/index/strtree/BoundablePair.cpp
+++ b/src/index/strtree/BoundablePair.cpp
@@ -77,7 +77,8 @@ BoundablePair::isLeaves() const
bool
BoundablePair::isComposite(const Boundable* item)
{
- return dynamic_cast<const AbstractNode*>(item) != nullptr;
+ return !(item->isLeaf());
+ // return dynamic_cast<const AbstractNode*>(item) != nullptr;
}
double
diff --git a/src/operation/distance/FacetSequence.cpp b/src/operation/distance/FacetSequence.cpp
index 6633e81..e1589c7 100644
--- a/src/operation/distance/FacetSequence.cpp
+++ b/src/operation/distance/FacetSequence.cpp
@@ -95,6 +95,9 @@ FacetSequence::nearestLocations(const FacetSequence& facetSeq) const
const Coordinate& seqPt = facetSeq.pts->getAt(facetSeq.start);
GeometryLocation gl1(geom, start, pt);
GeometryLocation gl2(facetSeq.geom, facetSeq.start, seqPt);
+ locs.clear();
+ locs.push_back(gl1);
+ locs.push_back(gl2);
}
else if (isPointThis) {
const Coordinate& pt = pts->getAt(start);
@@ -120,15 +123,11 @@ FacetSequence::computeDistancePointLine(const Coordinate& pt,
std::vector<GeometryLocation> *locs) const
{
double minDistance = std::numeric_limits<double>::infinity();
- double dist;
for(size_t i = facetSeq.start; i < facetSeq.end - 1; i++) {
const Coordinate& q0 = facetSeq.pts->getAt(i);
const Coordinate& q1 = facetSeq.pts->getAt(i + 1);
- dist = Distance::pointToSegment(pt, q0, q1);
- if(dist == 0.0) {
- return dist;
- }
+ double dist = Distance::pointToSegment(pt, q0, q1);
if(dist < minDistance) {
minDistance = dist;
if (locs != nullptr) {
diff --git a/src/operation/overlay/snap/LineStringSnapper.cpp b/src/operation/overlay/snap/LineStringSnapper.cpp
index aadef3b..9b89580 100644
--- a/src/operation/overlay/snap/LineStringSnapper.cpp
+++ b/src/operation/overlay/snap/LineStringSnapper.cpp
@@ -24,6 +24,7 @@
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateList.h>
+#include <geos/geom/Envelope.h>
#include <geos/util/UniqueCoordinateArrayFilter.h>
#include <geos/geom/LineSegment.h>
#include <geos/util/Interrupt.h>
@@ -431,6 +432,10 @@ LineStringSnapper::findSegmentToSnap(
}
}
+ if (Envelope::distanceSquaredToCoordinate(snapPt, seg.p0, seg.p1) >= minDist*minDist) {
+ continue;
+ }
+
double dist = seg.distance(snapPt);
if(dist >= minDist) {
#if GEOS_DEBUG
diff --git a/tests/unit/Makefile.am b/tests/unit/Makefile.am
index 3c836fc..cef4374 100644
--- a/tests/unit/Makefile.am
+++ b/tests/unit/Makefile.am
@@ -46,6 +46,8 @@ geos_unit_SOURCES = \
algorithm/CGAlgorithms/isPointInRingTest.cpp \
algorithm/CGAlgorithms/signedAreaTest.cpp \
algorithm/ConvexHullTest.cpp \
+ algorithm/construct/LargestEmptyCircleTest.cpp \
+ algorithm/construct/MaximumInscribedCircleTest.cpp \
algorithm/distance/DiscreteFrechetDistanceTest.cpp \
algorithm/distance/DiscreteHausdorffDistanceTest.cpp \
algorithm/InteriorPointAreaTest.cpp \
@@ -87,6 +89,7 @@ geos_unit_SOURCES = \
capi/GEOSisValidDetailTest.cpp \
capi/GEOSLineString_PointTest.cpp \
capi/GEOSMakeValidTest.cpp \
+ capi/GEOSMaximumInscribedCircleTest.cpp \
capi/GEOSMinimumBoundingCircleTest.cpp \
capi/GEOSMinimumClearanceTest.cpp \
capi/GEOSMinimumRectangleTest.cpp \
diff --git a/tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp b/tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp
new file mode 100644
index 0000000..a424772
--- /dev/null
+++ b/tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp
@@ -0,0 +1,243 @@
+//
+// Test Suite for geos::algorithm::construct::LargestEmptyCircle
+
+
+#include <tut/tut.hpp>
+// geos
+#include <geos/operation/distance/IndexedFacetDistance.h>
+#include <geos/algorithm/construct/LargestEmptyCircle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateArraySequence.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/PrecisionModel.h>
+#include <geos/io/WKTReader.h>
+#include <geos/io/WKTWriter.h>
+// std
+#include <sstream>
+#include <string>
+#include <memory>
+
+
+
+using namespace geos;
+using namespace geos::geom;
+
+using geos::algorithm::construct::LargestEmptyCircle;
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_lec_data {
+ geos::geom::Geometry* geom_;
+ geos::geom::PrecisionModel pm_;
+ geos::geom::GeometryFactory::Ptr factory_;
+ geos::io::WKTReader reader_;
+ geos::io::WKTWriter writer_;
+
+ test_lec_data():
+ geom_(nullptr),
+ pm_(geos::geom::PrecisionModel::FLOATING),
+ factory_(GeometryFactory::create(&pm_, 0)),
+ reader_(factory_.get())
+ {}
+
+ ~test_lec_data()
+ {
+ factory_->destroyGeometry(geom_);
+ geom_ = nullptr;
+ }
+
+ void
+ ensure_equals_coordinate(const geos::geom::Coordinate &lhs,
+ const geos::geom::Coordinate &rhs, double tolerance)
+ {
+ ensure_equals("x coordinate does not match", lhs.x, rhs.x, tolerance);
+ ensure_equals("y coordinate does not match", lhs.y, rhs.y, tolerance);
+ }
+
+ void
+ checkCircle(const Geometry *geom, double build_tolerance, double x, double y, double expectedRadius)
+ {
+ double tolerance = 2*build_tolerance;
+ LargestEmptyCircle lec(geom, tolerance);
+ std::unique_ptr<Point> centerPoint = lec.getCenter();
+ const Coordinate* centerPt = centerPoint->getCoordinate();
+ Coordinate expectedCenter(x, y);
+ ensure_equals_coordinate(*centerPt, expectedCenter, tolerance);
+ std::unique_ptr<LineString> radiusLine = lec.getRadiusLine();
+ double actualRadius = radiusLine->getLength();
+ ensure_equals("radius", actualRadius, expectedRadius, tolerance);
+ const Coordinate& linePt0 = radiusLine->getCoordinateN(0);
+ const Coordinate& linePt1 = radiusLine->getCoordinateN(1);
+
+ // std::cout << std::endl;
+ // std::cout << writer_.write(geom) << std::endl;
+ // std::cout << writer_.write(radiusLine.get()) << std::endl;
+
+ ensure_equals_coordinate(*centerPt, linePt0, tolerance);
+ const Coordinate* radiusPt = lec.getRadiusPoint()->getCoordinate();
+ ensure_equals_coordinate(*radiusPt, linePt1, tolerance);
+ }
+
+ void
+ checkCircleZeroRadius(const Geometry *geom, double tolerance)
+ {
+ LargestEmptyCircle lec(geom, tolerance);
+ std::unique_ptr<Point> centerPoint = lec.getCenter();
+ const Coordinate* centerPt = centerPoint->getCoordinate();
+ std::unique_ptr<LineString> radiusLine = lec.getRadiusLine();
+ double actualRadius = radiusLine->getLength();
+ ensure_equals("radius", actualRadius, 0.0, tolerance);
+
+ const Coordinate& linePt0 = radiusLine->getCoordinateN(0);
+ const Coordinate& linePt1 = radiusLine->getCoordinateN(1);
+
+ ensure_equals_coordinate(*centerPt, linePt0, tolerance);
+ const Coordinate* radiusPt = lec.getRadiusPoint()->getCoordinate();
+ ensure_equals_coordinate(*radiusPt, linePt1, tolerance);
+
+
+ }
+
+ void
+ checkCircleZeroRadius(std::string wkt, double tolerance)
+ {
+ std::unique_ptr<Geometry> geom(reader_.read(wkt));
+ checkCircleZeroRadius(geom.get(), tolerance);
+ }
+
+ void
+ checkCircle(std::string wkt, double tolerance, double x, double y, double expectedRadius)
+ {
+ std::unique_ptr<Geometry> geom(reader_.read(wkt));
+ checkCircle(geom.get(), tolerance, x, y, expectedRadius);
+ }
+
+};
+
+typedef test_group<test_lec_data> group;
+typedef group::object object;
+
+group test_lec_group("geos::algorithm::construct::LargestEmptyCircle");
+
+//
+// testPointsSquare
+//
+template<>
+template<>
+void object::test<1>
+()
+{
+
+ checkCircle("MULTIPOINT ((100 100), (100 200), (200 200), (200 100))",
+ 0.01, 150, 150, 70.71 );}
+
+//
+// testPointsTriangleOnHull
+//
+template<>
+template<>
+void object::test<2>
+()
+{
+ checkCircle("MULTIPOINT ((100 100), (300 100), (150 50))",
+ 0.01, 216.66, 99.99, 83.33 );
+}
+
+//
+// testPointsTriangleInterior
+//
+template<>
+template<>
+void object::test<3>
+()
+{
+ checkCircle("MULTIPOINT ((100 100), (300 100), (200 250))",
+ 0.01, 200.00, 141.66, 108.33 );
+}
+
+//
+// testLinesOpenDiamond
+//
+template<>
+template<>
+void object::test<4>
+()
+{
+
+ checkCircle("MULTILINESTRING ((50 100, 150 50), (250 50, 350 100), (350 150, 250 200), (50 150, 150 200))",
+ 0.01, 200, 125, 90.13 );}
+
+//
+// testLinesCrossed
+//
+template<>
+template<>
+void object::test<5>
+()
+{
+ checkCircle("MULTILINESTRING ((100 100, 300 300), (100 200, 300 0))",
+ 0.01, 299.99, 150.00, 106.05 );
+}
+
+
+//
+// testLinesZigzag
+//
+template<>
+template<>
+void object::test<6>
+()
+{
+
+ checkCircle("MULTILINESTRING ((100 100, 200 150, 100 200, 250 250, 100 300, 300 350, 100 400), (50 400, 0 350, 50 300, 0 250, 50 200, 0 150, 50 100))",
+ 0.01, 77.52, 349.99, 54.81 );
+}
+
+//
+// testPointsLinesTriangle
+//
+template<>
+template<>
+void object::test<7>
+()
+{
+checkCircle("GEOMETRYCOLLECTION (LINESTRING (100 100, 300 100), POINT (250 200))",
+ 0.01, 196.49, 164.31, 64.31 );
+}
+
+
+//
+// testPointsLinesTriangle
+//
+template<>
+template<>
+void object::test<8>
+()
+{
+checkCircleZeroRadius("POINT (100 100)",
+ 0.01 );
+}
+
+//
+// testLineFlat
+//
+template<>
+template<>
+void object::test<9>
+()
+{
+ checkCircleZeroRadius("LINESTRING (0 0, 50 50)",
+ 0.01 );
+}
+
+
+} // namespace tut
+
+
diff --git a/tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp b/tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp
new file mode 100644
index 0000000..59c674d
--- /dev/null
+++ b/tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp
@@ -0,0 +1,180 @@
+//
+// Test Suite for geos::algorithm::construct::MaximumInscribedCircle
+
+#include <tut/tut.hpp>
+// geos
+#include <geos/operation/distance/IndexedFacetDistance.h>
+#include <geos/algorithm/construct/MaximumInscribedCircle.h>
+#include <geos/geom/Coordinate.h>
+#include <geos/geom/CoordinateArraySequence.h>
+#include <geos/geom/Dimension.h>
+#include <geos/geom/Geometry.h>
+#include <geos/geom/GeometryFactory.h>
+#include <geos/geom/LineString.h>
+#include <geos/geom/PrecisionModel.h>
+#include <geos/io/WKTReader.h>
+#include <geos/io/WKTWriter.h>
+// std
+#include <sstream>
+#include <string>
+#include <memory>
+
+
+
+using namespace geos;
+using namespace geos::geom;
+
+using geos::algorithm::construct::MaximumInscribedCircle;
+
+namespace tut {
+//
+// Test Group
+//
+
+// dummy data, not used
+struct test_mic_data {
+ geos::geom::Geometry* geom_;
+ geos::geom::PrecisionModel pm_;
+ geos::geom::GeometryFactory::Ptr factory_;
+ geos::io::WKTReader reader_;
+ geos::io::WKTWriter writer_;
+
+ test_mic_data():
+ geom_(nullptr),
+ pm_(geos::geom::PrecisionModel::FLOATING),
+ factory_(GeometryFactory::create(&pm_, 0)),
+ reader_(factory_.get())
+ {}
+
+ ~test_mic_data()
+ {
+ factory_->destroyGeometry(geom_);
+ geom_ = nullptr;
+ }
+
+ void
+ ensure_equals_coordinate(const geos::geom::Coordinate &lhs,
+ const geos::geom::Coordinate &rhs, double tolerance)
+ {
+ ensure_equals("x coordinate does not match", lhs.x, rhs.x, tolerance);
+ ensure_equals("y coordinate does not match", lhs.y, rhs.y, tolerance);
+ }
+
+ void
+ checkCircle(const Geometry *geom, double build_tolerance, double x, double y, double expectedRadius)
+ {
+ double tolerance = 2*build_tolerance;
+ MaximumInscribedCircle mic(geom, tolerance);
+ std::unique_ptr<Point> centerPoint = mic.getCenter();
+ const Coordinate* centerPt = centerPoint->getCoordinate();
+ Coordinate expectedCenter(x, y);
+ ensure_equals_coordinate(*centerPt, expectedCenter, tolerance);
+ std::unique_ptr<LineString> radiusLine = mic.getRadiusLine();
+ double actualRadius = radiusLine->getLength();
+ ensure_equals("radius", actualRadius, expectedRadius, tolerance);
+ const Coordinate& linePt0 = radiusLine->getCoordinateN(0);
+ const Coordinate& linePt1 = radiusLine->getCoordinateN(1);
+
+ // std::cout << std::endl;
+ // std::cout << writer_.write(geom) << std::endl;
+ // std::cout << writer_.write(radiusLine.get()) << std::endl;
+
+ ensure_equals_coordinate(*centerPt, linePt0, tolerance);
+ const Coordinate* radiusPt = mic.getRadiusPoint()->getCoordinate();
+ ensure_equals_coordinate(*radiusPt, linePt1, tolerance);
+ }
+
+ void
+ checkCircle(std::string wkt, double tolerance, double x, double y, double expectedRadius)
+ {
+ std::unique_ptr<Geometry> geom(reader_.read(wkt));
+ checkCircle(geom.get(), tolerance, x, y, expectedRadius);
+ }
+
+};
+
+typedef test_group<test_mic_data> group;
+typedef group::object object;
+
+group test_mic_group("geos::algorithm::construct::MaximumInscribedCircle");
+
+//
+// testSquare
+//
+template<>
+template<>
+void object::test<1>
+()
+{
+ checkCircle("POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))",
+ 0.001, 150, 150, 50);
+}
+
+//
+// testDiamond
+//
+template<>
+template<>
+void object::test<2>
+()
+{
+ checkCircle("POLYGON ((150 250, 50 150, 150 50, 250 150, 150 250))",
+ 0.001, 150, 150, 70.71 );
+}
+
+//
+// testCircle
+//
+template<>
+template<>
+void object::test<3>
+()
+{
+ std::unique_ptr<Geometry> geom(reader_.read("POINT (100 100)"));
+ std::unique_ptr<Geometry> circle(geom->buffer(100, 20));
+ // MIC radius is less than 100 because buffer boundary segments lie inside circle
+ checkCircle(circle.get(), 0.01, 100, 100, 99.9229);
+}
+
+//
+// testDoubleKite
+//
+template<>
+template<>
+void object::test<4>
+()
+{
+ checkCircle("MULTIPOLYGON (((150 200, 100 150, 150 100, 250 150, 150 200)), ((400 250, 300 150, 400 50, 560 150, 400 250)))",
+ 0.001, 411.38877, 149.9996185, 78.7634662 );
+}
+
+//
+// Invalid polygon collapsed to a line
+//
+template<>
+template<>
+void object::test<5>
+()
+{
+ checkCircle("POLYGON ((100 100, 200 200, 100 100, 100 100))",
+ 0.01, 150, 150, 0 );
+}
+
+
+// //
+// // Invalid polygon collapsed to a point
+// //
+template<>
+template<>
+void object::test<6>
+()
+{
+ checkCircle("POLYGON ((100 100, 100 100, 100 100, 100 100))",
+ 0.01, 100, 100, 0 );
+}
+
+
+
+
+} // namespace tut
+
diff --git a/tests/unit/capi/GEOSMaximumInscribedCircleTest.cpp b/tests/unit/capi/GEOSMaximumInscribedCircleTest.cpp
new file mode 100644
index 0000000..210f9e8
--- /dev/null
+++ b/tests/unit/capi/GEOSMaximumInscribedCircleTest.cpp
@@ -0,0 +1,106 @@
+// $Id$
+//
+// Test Suite for C-API GEOSGetCentroid
+
+#include <tut/tut.hpp>
+// geos
+#include <geos_c.h>
+// std
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+#include <cstring>
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used in test cases.
+struct test_capimaximuminscribedcircle_data {
+ GEOSGeometry* geom1_;
+ GEOSGeometry* geom2_;
+ GEOSWKTWriter* wktw_;
+ char* wkt_;
+ double area_;
+
+ 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_capimaximuminscribedcircle_data()
+ : geom1_(nullptr), geom2_(nullptr), wkt_(nullptr)
+ {
+ initGEOS(notice, notice);
+ wktw_ = GEOSWKTWriter_create();
+ GEOSWKTWriter_setTrim(wktw_, 1);
+ GEOSWKTWriter_setRoundingPrecision(wktw_, 8);
+ }
+
+ ~test_capimaximuminscribedcircle_data()
+ {
+ GEOSGeom_destroy(geom1_);
+ GEOSGeom_destroy(geom2_);
+ GEOSWKTWriter_destroy(wktw_);
+ GEOSFree(wkt_);
+ geom1_ = nullptr;
+ geom2_ = nullptr;
+ wkt_ = nullptr;
+ finishGEOS();
+ }
+
+};
+
+typedef test_group<test_capimaximuminscribedcircle_data> group;
+typedef group::object object;
+
+group test_capimaximuminscribedcircle_group("capi::GEOSMaximumInscribedCircle");
+
+//
+// Test Cases
+//
+
+// Single point
+template<>
+template<>
+void object::test<1>
+()
+{
+ geom1_ = GEOSGeomFromWKT("POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))");
+ ensure(nullptr != geom1_);
+ geom2_ = GEOSMaximumInscribedCircle(geom1_, 0.001);
+ ensure(nullptr != geom2_);
+
+ wkt_ = GEOSWKTWriter_write(wktw_, geom2_);
+
+ ensure_equals(std::string(wkt_), std::string("LINESTRING (150 150, 150 200)"));
+}
+
+// Single point
+template<>
+template<>
+void object::test<2>
+()
+{
+ geom1_ = GEOSGeomFromWKT("MULTIPOINT ((100 100), (100 200), (200 200), (200 100))");
+ ensure(nullptr != geom1_);
+ geom2_ = GEOSLargestEmptyCircle(geom1_, 0.001);
+ ensure(nullptr != geom2_);
+
+ wkt_ = GEOSWKTWriter_write(wktw_, geom2_);
+
+ ensure_equals(std::string(wkt_), std::string("LINESTRING (150 150, 100 100)"));
+}
+
+
+} // namespace tut
+
diff --git a/tests/unit/geom/EnvelopeTest.cpp b/tests/unit/geom/EnvelopeTest.cpp
index a68c1c6..bd98f6e 100644
--- a/tests/unit/geom/EnvelopeTest.cpp
+++ b/tests/unit/geom/EnvelopeTest.cpp
@@ -242,5 +242,67 @@ void object::test<9>
ensure(empty == exemplar);
}
+// Test point-to-envelope distance
+template<>
+template<>
+void object::test<10>
+()
+{
+ using geos::geom::Coordinate;
+ using geos::geom::Envelope;
+
+ // Create a 5x5 grid of points and use them to test various
+ // spatial arrangements of the envelope and test point
+
+ // 0 1 2 3 4
+ // 5 6 7 8 9
+ // 10 11 12 13 14
+ // 15 16 17 18 19
+ // 20 21 22 23 24
+ std::vector<Coordinate> c(25);
+
+ std::cout<<std::endl;
+ for (size_t i = 0; i < c.size(); i++) {
+ c[i].x = static_cast<double>(i % 5);
+ c[i].y = static_cast<double>(5 - (i / 5));
+ std::cout<< c[i] << std::endl;
+ }
+
+ // point contained in envelope
+ ensure_equals( Envelope::distanceToCoordinate(c[18], c[22], c[9]), 0);
+ ensure_equals( Envelope::distanceToCoordinate(c[18], c[14], c[18]), 0);
+ ensure_equals( Envelope::distanceToCoordinate(c[18], c[14], c[17]), 0);
+ ensure_equals( Envelope::distanceToCoordinate(c[18], c[19], c[22]), 0);
+
+ // envelope above point
+ ensure_equals(Envelope::distanceToCoordinate(c[17], c[5], c[4]), 2);
+
+ // envelope below point
+ ensure_equals(Envelope::distanceToCoordinate(c[7], c[20], c[19]), 2);
+
+ // envelope left of point
+ ensure_equals(Envelope::distanceToCoordinate(c[13], c[20], c[11]), 2);
+
+ // envelope right of point
+ ensure_equals(Envelope::distanceToCoordinate(c[5], c[9], c[8]), 3);
+
+ // envelope upper-left of point
+ ensure_equals(Envelope::distanceToCoordinate(c[17], c[6], c[0]),
+ c[17].distance(c[6]));
+
+ // envelope upper-right of point
+ ensure_equals(Envelope::distanceToCoordinate(c[21], c[9], c[13]),
+ c[21].distance(c[13]));
+
+ // envelope lower-left of point
+ ensure_equals(Envelope::distanceToCoordinate(c[3], c[10], c[21]),
+ c[3].distance(c[11]));
+
+ // envelope lower-right of point
+ ensure_equals(Envelope::distanceToCoordinate(c[6], c[12], c[14]),
+ c[6].distance(c[12]));
+
+}
+
} // namespace tut
diff --git a/tests/unit/io/WKTReaderTest.cpp b/tests/unit/io/WKTReaderTest.cpp
index fe370d0..a50fe6f 100644
--- a/tests/unit/io/WKTReaderTest.cpp
+++ b/tests/unit/io/WKTReaderTest.cpp
@@ -126,14 +126,14 @@ void object::test<6>
try {
geom = wktreader.read("POLYGON( EMPTY, (1 1,2 2,1 2,1 1))");
- fail("Didn't get expected exception");
+ fail("Did not get expected exception");
}
catch(const geos::util::IllegalArgumentException& ex) {
ensure("Got expected exception", true);
- ex.what();
+ (void)(ex.what());
}
catch(...) {
- fail("Got unexpected execpetion.");
+ fail("Got unexpected exception");
}
}
@@ -162,7 +162,7 @@ void object::test<7>
}
catch(const geos::util::IllegalArgumentException& ex) {
ensure("Got expected exception", true);
- ex.what();
+ (void)(ex.what());
}
catch(...) {
fail("Got unexpected exception");
diff --git a/tests/unit/operation/distance/IndexedFacetDistanceTest.cpp b/tests/unit/operation/distance/IndexedFacetDistanceTest.cpp
index 3dcabb5..a0e2c8a 100644
--- a/tests/unit/operation/distance/IndexedFacetDistanceTest.cpp
+++ b/tests/unit/operation/distance/IndexedFacetDistanceTest.cpp
@@ -1,12 +1,18 @@
//
// Test Suite for geos::operation::distance::DistanceOp class.
+#define _USE_MATH_DEFINES
+
+// std
+#include <memory>
+#include <string>
+#include <vector>
+#include <cmath>
+
// tut
#include <tut/tut.hpp>
// geos
-#include <geos/algorithm/PointLocator.h>
-#include <geos/io/WKTReader.h>
-#include <geos/io/WKBReader.h>
+#include <geos/profiler.h>
#include <geos/constants.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateArraySequence.h>
@@ -16,12 +22,14 @@
#include <geos/geom/LineSegment.h>
#include <geos/geom/LineString.h>
#include <geos/geom/PrecisionModel.h>
+#include <geos/io/WKTReader.h>
+#include <geos/io/WKTWriter.h>
#include <geos/operation/distance/DistanceOp.h>
#include <geos/operation/distance/IndexedFacetDistance.h>
-// std
-#include <memory>
-#include <string>
-#include <vector>
+
+#ifndef M_PI
+#define M_PI 3.14159265358979323846
+#endif
namespace tut {
//
@@ -30,13 +38,21 @@ namespace tut {
// Common data used by tests
struct test_facetdistanceop_data {
- geos::io::WKTReader wktreader;
typedef std::unique_ptr<geos::geom::Geometry> GeomPtr;
typedef std::unique_ptr<geos::geom::CoordinateSequence> CSPtr;
+ geos::io::WKTWriter _wktwriter;
+ geos::geom::PrecisionModel _pm;
+ geos::geom::GeometryFactory::Ptr _factory;
+ geos::io::WKTReader _wktreader;
+
+
test_facetdistanceop_data()
- : wktreader()
+ : _wktwriter()
+ , _pm(geos::geom::PrecisionModel::FLOATING)
+ , _factory(geos::geom::GeometryFactory::create(&_pm, 0))
+ , _wktreader(_factory.get())
{}
void
@@ -44,8 +60,8 @@ struct test_facetdistanceop_data {
geos::geom::Coordinate& p1, geos::geom::Coordinate& p2)
{
using geos::operation::distance::IndexedFacetDistance;
- GeomPtr g1(wktreader.read(wkt1));
- GeomPtr g2(wktreader.read(wkt2));
+ GeomPtr g1(_wktreader.read(wkt1));
+ GeomPtr g2(_wktreader.read(wkt2));
std::vector<geos::geom::Coordinate> pts;
pts = IndexedFacetDistance::nearestPoints(g1.get(), g2.get());
ensure(fabs(pts[0].distance(pts[1])-distance) < 1e08);
@@ -56,6 +72,37 @@ struct test_facetdistanceop_data {
return;
}
+ int
+ angle2sincircle(double theta_deg, double radius, double amplitude, double* x, double* y)
+ {
+ // π=π
+πsin(ππ)
+ int n = 16;
+ double theta = theta_deg * M_PI / 180.0;
+ double a = radius * amplitude;
+ double r = radius + a * std::sin(n*theta);
+ if (x) *x = r * std::cos(theta);
+ if (y) *y = r * std::sin(theta);
+ return 1;
+ }
+
+
+ std::unique_ptr<geos::geom::LineString>
+ makeSinCircle(size_t nvertices, double radius, double amplitude)
+ {
+ geos::geom::CoordinateArraySequence cs;
+ std::vector<geos::geom::Coordinate> coords;
+ for (size_t i = 0; i < nvertices; i++) {
+ geos::geom::Coordinate c;
+ angle2sincircle(i*360.0/nvertices, radius, amplitude, &c.x, &c.y);
+ cs.add(c);
+ }
+
+ std::unique_ptr<geos::geom::LineString> ls(_factory->createLineString(cs));
+ return ls;
+ }
+
+
+
};
typedef test_group<test_facetdistanceop_data> group;
@@ -64,7 +111,6 @@ typedef group::object object;
group test_facetdistanceop_group("geos::operation::distance::IndexedFacetDistance");
-
//
// Test Cases
//
@@ -77,8 +123,8 @@ void object::test<1>
std::string wkt0("POINT(0 0)");
std::string wkt1("POINT(10 0)");
- GeomPtr g0(wktreader.read(wkt0));
- GeomPtr g1(wktreader.read(wkt1));
+ GeomPtr g0(_wktreader.read(wkt0));
+ GeomPtr g1(_wktreader.read(wkt1));
double d = IndexedFacetDistance::distance(g0.get(), g1.get());
ensure_equals(d, 10);
}
@@ -158,8 +204,8 @@ void object::test<6>
std::string wkt0("POLYGON((100 200, 200 200, 200 100, 100 100, 100 200))");
std::string wkt1("POINT(150 150)");
- GeomPtr g0(wktreader.read(wkt0));
- GeomPtr g1(wktreader.read(wkt1));
+ GeomPtr g0(_wktreader.read(wkt0));
+ GeomPtr g1(_wktreader.read(wkt1));
double d = IndexedFacetDistance::distance(g0.get(), g1.get());
ensure_equals(d, 50);
}
@@ -174,12 +220,121 @@ void object::test<7>
std::string wkt0("POLYGON((100 200, 200 200, 200 100, 100 100, 100 200))");
std::string wkt1("POINT(150 150)");
- GeomPtr g0(wktreader.read(wkt0));
- GeomPtr g1(wktreader.read(wkt1));
+ GeomPtr g0(_wktreader.read(wkt0));
+ GeomPtr g1(_wktreader.read(wkt1));
double d = IndexedFacetDistance::distance(g0.get(), g1.get());
ensure_equals(d, 50);
}
+template<>
+template<>
+void object::test<8>
+()
+{
+ using geos::operation::distance::IndexedFacetDistance;
+ std::string wkt0("POLYGON((100 200, 200 200, 200 100, 100 100, 100 200))");
+ std::string wkt1("POINT(150 150)");
+ GeomPtr g0(_wktreader.read(wkt0));
+ GeomPtr g1(_wktreader.read(wkt1));
+ IndexedFacetDistance ifd(g0.get());
+ double d = ifd.distance(g1.get());
+ ensure_equals(d, 50);
+}
+
+// Invalid polygon collapsed to a line
+template<>
+template<>
+void object::test<9>
+()
+{
+ using geos::operation::distance::IndexedFacetDistance;
+ std::string wkt0("POLYGON((100 100, 200 200, 100 100, 100 100))");
+ std::string wkt1("POINT(150 150)");
+ GeomPtr g0(_wktreader.read(wkt0));
+ GeomPtr g1(_wktreader.read(wkt1));
+ IndexedFacetDistance ifd(g0.get());
+ double d = ifd.distance(g1.get());
+ ensure_equals("incorrect distance", d, 0.0, 0.001);
+
+ std::vector<geos::geom::Coordinate> nearestPts = ifd.nearestPoints(g1.get());
+ ensure_equals("nearest points x", nearestPts[0].x, nearestPts[1].x, 0.00001);
+ ensure_equals("nearest points y", nearestPts[0].y, nearestPts[1].y, 0.00001);
+}
+
+// Invalid polygon collapsed to a line
+template<>
+template<>
+void object::test<10>
+()
+{
+ int npoints = 1000; // vertices in sinstar test shape
+ int ncells = 100; // number of colums/rows in test grid square
+
+ double radius = 100.0;
+ double amplitude = 0.3; // how far the sin deviates from perfect circle (0.0)
+ double width = radius * (1+amplitude); // total radius of shape
+ double cellsize = 2.0*width/ncells; // how big a cell is
+
+ // Build a sine star of the requested size and prepare
+ // the IndexedFacetDistance for it
+ using geos::operation::distance::IndexedFacetDistance;
+ std::unique_ptr<geos::geom::LineString> ls = makeSinCircle(npoints, radius, amplitude);
+ IndexedFacetDistance ifd(ls.get());
+ // std::string wkt = _wktwriter.write(ls.get());
+ // std::cout << wkt << std::endl;
+
+ // Build out the set of test points ahead of time so that
+ // point creation overhead isn't included in the test timings
+ std::vector<std::unique_ptr<geos::geom::Point>> pts;
+ for (double x = -width; x < width; x += cellsize) {
+ for (double y = -width; y < width; y += cellsize) {
+ geos::geom::Coordinate c(x, y);
+ pts.emplace(pts.end(), _factory->createPoint(c));
+ }
+ }
+
+ geos::util::Profiler prof;
+
+ enum modes {
+ TestIndexedFacetDistance = 0,
+ TestGeometryDistance = 1,
+ TestBoth = 2
+ };
+
+ std::vector<std::string> profiles = {"TestIndexedFacetDistance",
+ "TestGeometryDistance",
+ "TestBoth"};
+
+ bool perfTest = false;
+ std::vector<int> m = {TestBoth};
+ if (perfTest) {
+ m.push_back(TestIndexedFacetDistance);
+ m.push_back(TestGeometryDistance);
+ }
+
+ for (auto it = m.begin() ; it != m.end(); ++it)
+ {
+ prof.start(profiles[*it]);
+ for (size_t j = 0; j < pts.size(); j++) {
+ double dist_ifd = 0.0, dist_geom = 0.0;
+ if (*it == TestIndexedFacetDistance || *it == TestBoth)
+ dist_ifd = ifd.distance(pts[j].get());
+
+ if (*it == TestGeometryDistance || *it == TestBoth)
+ dist_geom = ls->distance(pts[j].get());
+
+ if (*it == TestBoth)
+ ensure_equals("distance", dist_ifd, dist_geom, 0.00001);
+ }
+
+ prof.stop(profiles[*it]);
+ }
+ if (perfTest) {
+ std::cout << "npoints=" << npoints << " ncells=" << ncells << std::endl;
+ std::cout << prof << std::endl;
+ }
+}
+
// TODO: finish the tests by adding:
diff --git a/tests/unit/operation/geounion/CoverageUnionTest.cpp b/tests/unit/operation/geounion/CoverageUnionTest.cpp
index 2312bfa..91de550 100644
--- a/tests/unit/operation/geounion/CoverageUnionTest.cpp
+++ b/tests/unit/operation/geounion/CoverageUnionTest.cpp
@@ -78,7 +78,9 @@ namespace tut {
try {
auto u1 = CoverageUnion::Union(coll.get());
fail();
- } catch(const geos::util::TopologyException & e) {}
+ } catch(const geos::util::TopologyException & e) {
+ (void)(e); // silence unused warning
+ }
}
template<>
diff --git a/tests/xmltester/XMLTester.cpp b/tests/xmltester/XMLTester.cpp
index b44cf31..11eae81 100644
--- a/tests/xmltester/XMLTester.cpp
+++ b/tests/xmltester/XMLTester.cpp
@@ -166,7 +166,12 @@ getIndent(unsigned int numIndents)
void
tolower(std::string& str)
{
- std::transform(str.begin(), str.end(), str.begin(), (int(*)(int))std::tolower);
+ std::transform(
+ str.begin(),
+ str.end(),
+ str.begin(),
+ [](unsigned char c){ return std::tolower(c); }
+ );
}
std::string
-----------------------------------------------------------------------
Summary of changes:
.azure-pipelines.yml | 39 +--
.gitignore | 1 +
NEWS | 6 +-
capi/geos_c.cpp | 12 +
capi/geos_c.h.in | 35 +++
capi/geos_ts_c.cpp | 24 ++
configure.ac | 4 +-
include/geos/algorithm/Makefile.am | 3 +-
.../geos/algorithm/construct/LargestEmptyCircle.h | 207 ++++++++++++++++
.../geos/{math => algorithm/construct}/Makefile.am | 5 +-
.../algorithm/construct/MaximumInscribedCircle.h | 198 +++++++++++++++
include/geos/index/strtree/Boundable.h | 3 +-
src/algorithm/Makefile.am | 6 +-
src/algorithm/construct/LargestEmptyCircle.cpp | 272 +++++++++++++++++++++
src/algorithm/construct/Makefile.am | 14 ++
src/algorithm/construct/MaximumInscribedCircle.cpp | 232 ++++++++++++++++++
src/index/strtree/BoundablePair.cpp | 3 +-
src/operation/distance/FacetSequence.cpp | 9 +-
tests/unit/Makefile.am | 3 +
.../algorithm/construct/LargestEmptyCircleTest.cpp | 243 ++++++++++++++++++
.../construct/MaximumInscribedCircleTest.cpp | 180 ++++++++++++++
tests/unit/capi/GEOSMaximumInscribedCircleTest.cpp | 106 ++++++++
tests/unit/io/WKTReaderTest.cpp | 8 +-
.../distance/IndexedFacetDistanceTest.cpp | 191 +++++++++++++--
.../unit/operation/geounion/CoverageUnionTest.cpp | 4 +-
tests/xmltester/XMLTester.cpp | 7 +-
26 files changed, 1740 insertions(+), 75 deletions(-)
create mode 100644 include/geos/algorithm/construct/LargestEmptyCircle.h
copy include/geos/{math => algorithm/construct}/Makefile.am (52%)
create mode 100644 include/geos/algorithm/construct/MaximumInscribedCircle.h
create mode 100644 src/algorithm/construct/LargestEmptyCircle.cpp
create mode 100644 src/algorithm/construct/Makefile.am
create mode 100644 src/algorithm/construct/MaximumInscribedCircle.cpp
create mode 100644 tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp
create mode 100644 tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp
create mode 100644 tests/unit/capi/GEOSMaximumInscribedCircleTest.cpp
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list