[geos-commits] [SCM] GEOS branch main updated. c11ff90e18cc14395977640154168623fa75c55a
git at osgeo.org
git at osgeo.org
Fri Apr 28 14:37:40 PDT 2023
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, main has been updated
via c11ff90e18cc14395977640154168623fa75c55a (commit)
from 23ae165cf1e6a6fc0f98a1360f986fafa0932923 (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 c11ff90e18cc14395977640154168623fa75c55a
Author: Martin Davis <mtnclimb at gmail.com>
Date: Fri Apr 28 14:36:36 2023 -0700
Fix and improve MIC and LEC (#883)
diff --git a/include/geos/algorithm/construct/MaximumInscribedCircle.h b/include/geos/algorithm/construct/MaximumInscribedCircle.h
index fa1119151..1da12c038 100644
--- a/include/geos/algorithm/construct/MaximumInscribedCircle.h
+++ b/include/geos/algorithm/construct/MaximumInscribedCircle.h
@@ -107,6 +107,22 @@ public:
*/
static std::unique_ptr<geom::LineString> getRadiusLine(const geom::Geometry* polygonal, double tolerance);
+ /**
+ * Computes the maximum number of iterations allowed.
+ * Uses a heuristic based on the area of the input geometry
+ * and the tolerance distance.
+ * The number of tolerance-sized cells that cover the input geometry area
+ * is computed, times a safety factor.
+ * This prevents massive numbers of iterations and created cells
+ * for casees where the input geometry has extremely small area
+ * (e.g. is very thin).
+ *
+ * @param geom the input geometry
+ * @param toleranceDist the tolerance distance
+ * @return the maximum number of iterations allowed
+ */
+ static std::size_t computeMaximumIterations(const geom::Geometry* geom, double toleranceDist);
+
private:
/* private members */
@@ -170,22 +186,32 @@ private:
{
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;
}
+
+ /**
+ * The Cell priority queue is sorted by the natural order of maxDistance.
+ * std::priority_queue sorts with largest first,
+ * which is what is needed for this algorithm.
+ */
+ using CellQueue = std::priority_queue<Cell>;
};
- void createInitialGrid(const geom::Envelope* env, std::priority_queue<Cell>& cellQueue);
- Cell createCentroidCell(const geom::Geometry* geom);
+ void createInitialGrid(const geom::Envelope* env, Cell::CellQueue& cellQueue);
+ Cell createInteriorPointCell(const geom::Geometry* geom);
};
@@ -193,4 +219,3 @@ private:
} // geos::algorithm::construct
} // geos::algorithm
} // geos
-
diff --git a/src/algorithm/construct/LargestEmptyCircle.cpp b/src/algorithm/construct/LargestEmptyCircle.cpp
index e1c2fe5ab..68f50223a 100644
--- a/src/algorithm/construct/LargestEmptyCircle.cpp
+++ b/src/algorithm/construct/LargestEmptyCircle.cpp
@@ -18,6 +18,7 @@
**********************************************************************/
#include <geos/algorithm/construct/LargestEmptyCircle.h>
+#include <geos/algorithm/construct/MaximumInscribedCircle.h>
#include <geos/geom/Coordinate.h>
#include <geos/geom/CoordinateSequence.h>
#include <geos/geom/Envelope.h>
@@ -28,6 +29,7 @@
#include <geos/geom/MultiPolygon.h>
#include <geos/algorithm/locate/IndexedPointInAreaLocator.h>
#include <geos/operation/distance/IndexedFacetDistance.h>
+#include <geos/util/Interrupt.h>
#include <typeinfo> // for dynamic_cast
#include <cassert>
@@ -110,21 +112,20 @@ LargestEmptyCircle::getRadiusLine()
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));
- }
+ if (!env->isFinite()) {
+ throw util::GEOSException("Non-finite envelope encountered.");
}
+
+ double cellSize = std::max(env->getWidth(), env->getHeight());
+ double hSide = cellSize / 2.0;
+
+ // Collapsed geometries just end up using the centroid
+ // as the answer and skip all the other machinery
+ if (cellSize == 0) return;
+
+ CoordinateXY c;
+ env->centre(c);
+ cellQueue.emplace(c.x, c.y, hSide, distanceToConstraints(c.x, c.y));
}
/* private */
@@ -231,12 +232,18 @@ LargestEmptyCircle::compute()
* Carry out the branch-and-bound search
* of the cell space
*/
- while (!cellQueue.empty()) {
+ std::size_t maxIter = MaximumInscribedCircle::computeMaximumIterations(boundary.get(), tolerance);
+ std::size_t iterationCount = 0;
+ while (!cellQueue.empty() && iterationCount < maxIter) {
// pick the most promising cell from the queue
Cell cell = cellQueue.top();
cellQueue.pop();
+ if ((iterationCount++ % 1000) == 0) {
+ GEOS_CHECK_FOR_INTERRUPTS();
+ }
+
// update the center cell if the candidate is further from the constraints
if (cell.getDistance() > farthestCell.getDistance()) {
farthestCell = cell;
diff --git a/src/algorithm/construct/MaximumInscribedCircle.cpp b/src/algorithm/construct/MaximumInscribedCircle.cpp
index 0175788d9..b20d14a08 100644
--- a/src/algorithm/construct/MaximumInscribedCircle.cpp
+++ b/src/algorithm/construct/MaximumInscribedCircle.cpp
@@ -78,6 +78,18 @@ MaximumInscribedCircle::getRadiusLine(const Geometry* polygonal, double toleranc
return mic.getRadiusLine();
}
+/* public static */
+std::size_t
+MaximumInscribedCircle::computeMaximumIterations(const Geometry* geom, double toleranceDist)
+{
+ double diam = geom->getEnvelopeInternal()->getDiameter();
+ double ncells = diam / toleranceDist;
+ //-- Using log of ncells allows control over number of iterations
+ std::size_t factor = (std::size_t) std::log(ncells);
+ if (factor < 1) factor = 1;
+ return 2000 + 2000 * factor;
+}
+
/* public */
std::unique_ptr<Point>
MaximumInscribedCircle::getCenter()
@@ -106,33 +118,26 @@ MaximumInscribedCircle::getRadiusLine()
return factory->createLineString(std::move(cl));
}
+int INITIAL_GRID_SIDE = 25;
+
/* private */
void
-MaximumInscribedCircle::createInitialGrid(const Envelope* env, std::priority_queue<Cell>& cellQueue)
+MaximumInscribedCircle::createInitialGrid(const Envelope* env, Cell::CellQueue& cellQueue)
{
if (!env->isFinite()) {
throw util::GEOSException("Non-finite envelope encountered.");
}
- 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;
+ double cellSize = std::max(env->getWidth(), env->getHeight());
+ double hSide = cellSize / 2.0;
// Collapsed geometries just end up using the centroid
// as the answer and skip all the other machinery
if (cellSize == 0) return;
- // 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));
- }
- }
+ CoordinateXY c;
+ env->centre(c);
+ cellQueue.emplace(c.x, c.y, hSide, distanceToBoundary(c.x, c.y));
}
/* private */
@@ -157,11 +162,11 @@ MaximumInscribedCircle::distanceToBoundary(double x, double y)
/* private */
MaximumInscribedCircle::Cell
-MaximumInscribedCircle::createCentroidCell(const Geometry* geom)
+MaximumInscribedCircle::createInteriorPointCell(const Geometry* geom)
{
Coordinate c;
- geom->getCentroid(c);
- Cell cell(c.x, c.y, 0, distanceToBoundary(c));
+ std::unique_ptr<Point> p = geom->getInteriorPoint();
+ Cell cell(p->getX(), p->getY(), 0, distanceToBoundary(c));
return cell;
}
@@ -174,27 +179,33 @@ MaximumInscribedCircle::compute()
if (done) return;
// Priority queue of cells, ordered by maximum distance from boundary
- std::priority_queue<Cell> cellQueue;
+ Cell::CellQueue cellQueue;
createInitialGrid(inputGeom->getEnvelopeInternal(), cellQueue);
// use the area centroid as the initial candidate center point
- Cell farthestCell = createCentroidCell(inputGeom);
+ Cell farthestCell = createInteriorPointCell(inputGeom);
/**
* Carry out the branch-and-bound search
* of the cell space
*/
+ std::size_t maxIter = computeMaximumIterations(inputGeom, tolerance);
std::size_t iterationCount = 0;
- while (!cellQueue.empty()) {
+ while (!cellQueue.empty() && iterationCount < maxIter) {
// pick the most promising cell from the queue
Cell cell = cellQueue.top();
cellQueue.pop();
+// std::cout << iterationCount << "] Dist: " << cell.getDistance() << " size: " << cell.getHSize() << std::endl;
if ((iterationCount++ % 1000) == 0) {
GEOS_CHECK_FOR_INTERRUPTS();
}
+ //-- if cell must be closer than furthest, terminate since all remaining cells in queue are even closer.
+ if (cell.getMaxDistance() < farthestCell.getDistance())
+ break;
+
// update the center cell if the candidate is further from the boundary
if (cell.getDistance() > farthestCell.getDistance()) {
farthestCell = cell;
@@ -237,4 +248,3 @@ MaximumInscribedCircle::compute()
} // namespace geos.algorithm.construct
} // namespace geos.algorithm
} // namespace geos
-
diff --git a/tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp b/tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp
index 8711f0c61..959beb4a0 100644
--- a/tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp
+++ b/tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp
@@ -61,6 +61,24 @@ struct test_lec_data {
ensure_equals("y coordinate does not match", lhs.y, rhs.y, tolerance);
}
+ /**
+ * A coarse distance check, mainly testing
+ * that there is not a huge number of iterations.
+ * (This will be revealed by CI taking a very long time!)
+ */
+ void
+ checkCircle(std::string wktObstacles, double tolerance)
+ {
+ std::unique_ptr<Geometry> geom(reader_.read(wktObstacles));
+ LargestEmptyCircle lec(geom.get(), tolerance);
+ std::unique_ptr<Point> centerPoint = lec.getCenter();
+ double dist = geom->distance(centerPoint.get());
+ //std::cout << dist << std::endl;
+ std::unique_ptr<LineString> radiusLine = lec.getRadiusLine();
+ double actualRadius = radiusLine->getLength();
+ ensure(std::abs(actualRadius - dist) < 2 * tolerance);
+ }
+
void
checkCircle(const Geometry *obstacles, const Geometry *boundary, double build_tolerance, double x, double y, double expectedRadius)
{
@@ -202,8 +220,8 @@ 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 );
+ checkCircle("MULTILINESTRING ((100 100, 200 150, 100 200, 250 250, 100 300, 300 350, 100 400), (70 380, 0 350, 50 300, 0 250, 50 200, 0 150, 50 120))",
+ 0.01, 77.52, 249.99, 54.81 );
}
//
@@ -243,11 +261,21 @@ void object::test<9>
0.01 );
}
-// testBoundaryEmpty
+// testThinExtent
template<>
template<>
void object::test<10>
()
+{
+ checkCircle("MULTIPOINT ((100 100), (300 100), (200 100.1))",
+ 0.01 );
+}
+
+// testBoundaryEmpty
+template<>
+template<>
+void object::test<11>
+()
{
checkCircle("MULTIPOINT ((2 2), (8 8), (7 5))",
"POLYGON EMPTY",
@@ -257,7 +285,7 @@ void object::test<10>
// testBoundarySquare
template<>
template<>
-void object::test<11>
+void object::test<12>
()
{
checkCircle("MULTIPOINT ((2 2), (6 4), (8 8))",
@@ -268,7 +296,7 @@ void object::test<11>
//testBoundarySquareObstaclesOutside
template<>
template<>
-void object::test<12>
+void object::test<13>
()
{
checkCircle("MULTIPOINT ((10 10), (10 0))",
@@ -279,7 +307,7 @@ void object::test<12>
// testBoundaryMultiSquares
template<>
template<>
-void object::test<13>
+void object::test<14>
()
{
checkCircle("MULTIPOINT ((10 10), (10 0), (5 5))",
@@ -290,7 +318,7 @@ void object::test<13>
// testBoundaryAsObstacle
template<>
template<>
-void object::test<14>
+void object::test<15>
()
{
checkCircle("GEOMETRYCOLLECTION (POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9)), POINT (4 3), POINT (7 6))",
diff --git a/tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp b/tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp
index a64e4def8..7a7add925 100644
--- a/tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp
+++ b/tests/unit/algorithm/construct/MaximumInscribedCircleTest.cpp
@@ -59,10 +59,26 @@ struct test_mic_data {
ensure_equals("y coordinate does not match", lhs.y, rhs.y, tolerance);
}
+ /**
+ * A coarse distance check, mainly testing
+ * that there is not a huge number of iterations.
+ * (This will be revealed by CI taking a very long time!)
+ */
void
- checkCircle(const Geometry *geom, double build_tolerance, double x, double y, double expectedRadius)
+ checkCircle(std::string wkt, double tolerance)
+ {
+ std::unique_ptr<Geometry> geom(reader_.read(wkt));
+ MaximumInscribedCircle mic(geom.get(), tolerance);
+ std::unique_ptr<Point> centerPoint = mic.getCenter();
+ std::unique_ptr<Geometry> bdy = geom->getBoundary();
+ double dist = bdy->distance(centerPoint.get());
+ //std::cout << dist << std::endl;
+ ensure(dist < 2 * tolerance);
+ }
+
+ void
+ checkCircle(const Geometry *geom, double tolerance, double x, double y, double expectedRadius)
{
- double tolerance = 2*build_tolerance;
MaximumInscribedCircle mic(geom, tolerance);
std::unique_ptr<Point> centerPoint = mic.getCenter();
Coordinate centerPt(*centerPoint->getCoordinate());
@@ -203,7 +219,30 @@ void object::test<8>
} catch (const util::GEOSException & e) {}
}
+ /**
+ * Tests that a nearly flat geometry doesn't make the initial cell grid huge.
+ *
+ * See https://github.com/libgeos/geos/issues/875
+ */
+// testNearlyFlat
+template<>
+template<>
+void object::test<9>
+()
+{
+ checkCircle("POLYGON ((59.3 100.00000000000001, 99.7 100.00000000000001, 99.7 100, 59.3 100, 59.3 100.00000000000001))",
+ 0.01 );
+}
+
+// testVeryThin
+template<>
+template<>
+void object::test<10>
+()
+{
+ checkCircle("POLYGON ((100 100, 200 300, 300 100, 450 250, 300 99.999999, 200 299.99999, 100 100))",
+ 0.01 );
+}
} // namespace tut
-
-----------------------------------------------------------------------
Summary of changes:
.../algorithm/construct/MaximumInscribedCircle.h | 31 +++++++++++--
src/algorithm/construct/LargestEmptyCircle.cpp | 37 +++++++++------
src/algorithm/construct/MaximumInscribedCircle.cpp | 54 +++++++++++++---------
.../algorithm/construct/LargestEmptyCircleTest.cpp | 42 ++++++++++++++---
.../construct/MaximumInscribedCircleTest.cpp | 45 ++++++++++++++++--
5 files changed, 159 insertions(+), 50 deletions(-)
hooks/post-receive
--
GEOS
More information about the geos-commits
mailing list