Enhance LargestEmptyCircle boundary handling (#859)

include/geos/algorithm/construct/LargestEmptyCircle.h

@@ -116,18 +116,16 @@ class GEOS_DLL LargestEmptyCircle {
     /* private members */
     double tolerance;
     const geom::Geometry* obstacles;
+    std::unique_ptr<geom::Geometry> boundary;
     const geom::GeometryFactory* factory;
-    std::unique_ptr<geom::Geometry> boundary; // convexhull(obstacles)
+    geom::Envelope gridEnv;
     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).
@@ -140,6 +138,7 @@ class GEOS_DLL LargestEmptyCircle {
     */
     double distanceToConstraints(const geom::Coordinate& c);
     double distanceToConstraints(double x, double y);
+    void initBoundary();
     void compute();
 
     /* private class */

src/algorithm/construct/LargestEmptyCircle.cpp

@@ -35,13 +35,10 @@
 
 using namespace geos::geom;
 
-
 namespace geos {
 namespace algorithm { // geos.algorithm
 namespace construct { // geos.algorithm.construct
 
-
-
 LargestEmptyCircle::LargestEmptyCircle(const Geometry* p_obstacles, double p_tolerance)
     : LargestEmptyCircle(p_obstacles, nullptr, p_tolerance)
 {
@@ -54,29 +51,14 @@ LargestEmptyCircle::LargestEmptyCircle(const Geometry* p_obstacles, const Geomet
     , obstacleDistance(p_obstacles)
     , done(false)
 {
-    if (!p_boundary)
-    {
-        boundary = p_obstacles->convexHull();
-    }
-    else
-    {
-        boundary = p_boundary->clone();
-    }
-
     if (obstacles->isEmpty()) {
         throw util::IllegalArgumentException("Empty obstacles geometry is not supported");
     }
-    if (boundary->isEmpty()) {
-        throw util::IllegalArgumentException("Empty obstacles geometry is not supported");
-    }
-    if (!boundary->covers(obstacles)) {
-        throw util::IllegalArgumentException("Boundary geometry does not cover obstacles");
+    if (! p_boundary || p_boundary->isEmpty()) {
+        boundary = obstacles->convexHull();
     }
-
-    // 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()));
+    else {
+        boundary = p_boundary->clone();
     }
 }
 
@@ -211,15 +193,26 @@ LargestEmptyCircle::createCentroidCell(const Geometry* geom)
     return cell;
 }
 
+/* private */
+void
+LargestEmptyCircle::initBoundary()
+{
+    gridEnv = *(boundary->getEnvelopeInternal());
+    // 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()));
+    }
+}
 
 /* private */
 void
 LargestEmptyCircle::compute()
 {
-
     // check if already computed
     if (done) return;
 
+    initBoundary();
     // if ptLocator is not present then result is degenerate (represented as zero-radius circle)
     if (!ptLocator) {
         const Coordinate* pt = obstacles->getCoordinate();
@@ -231,7 +224,7 @@ LargestEmptyCircle::compute()
 
     // Priority queue of cells, ordered by decreasing distance from constraints
     std::priority_queue<Cell> cellQueue;
-    createInitialGrid(obstacles->getEnvelopeInternal(), cellQueue);
+    createInitialGrid(&gridEnv, cellQueue);
 
     Cell farthestCell = createCentroidCell(obstacles);
 
@@ -286,5 +279,3 @@ LargestEmptyCircle::compute()
 } // namespace geos.algorithm.construct
 } // namespace geos.algorithm
 } // namespace geos
-
-

tests/unit/algorithm/construct/LargestEmptyCircleTest.cpp

@@ -63,10 +63,10 @@ struct test_lec_data {
     }
 
     void
-    checkCircle(const Geometry *geom, double build_tolerance, double x, double y, double expectedRadius)
+    checkCircle(const Geometry *obstacles, const Geometry *boundary, double build_tolerance, double x, double y, double expectedRadius)
     {
         double tolerance = 2*build_tolerance;
-        LargestEmptyCircle lec(geom, tolerance);
+        LargestEmptyCircle lec(obstacles, boundary, tolerance);
         std::unique_ptr<Point> centerPoint = lec.getCenter();
         const Coordinate* centerPt = centerPoint->getCoordinate();
         Coordinate expectedCenter(x, y);
@@ -115,7 +115,15 @@ struct test_lec_data {
     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);
+        checkCircle(geom.get(), nullptr, tolerance, x, y, expectedRadius);
+    }
+
+    void
+    checkCircle(std::string wktObstacles, std::string wktBoundary, double tolerance, double x, double y, double expectedRadius)
+    {
+        std::unique_ptr<Geometry> obstacles(reader_.read(wktObstacles));
+        std::unique_ptr<Geometry> boundary(reader_.read(wktBoundary));
+        checkCircle(obstacles.get(), boundary.get(), tolerance, x, y, expectedRadius);
     }
 
 };
@@ -236,7 +244,59 @@ void object::test<9>
        0.01 );
 }
 
+// testBoundaryEmpty
+template<>
+template<>
+void object::test<10>
+()
+{
+ checkCircle("MULTIPOINT ((2 2), (8 8), (7 5))",
+        "POLYGON EMPTY",
+        0.01, 4.127, 4.127, 3 );
+}
 
-} // namespace tut
+// testBoundarySquare
+template<>
+template<>
+void object::test<11>
+()
+{
+    checkCircle("MULTIPOINT ((2 2), (6 4), (8 8))",
+        "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+        0.01, 1.00390625, 8.99609375, 7.065 );
+}
+
+//testBoundarySquareObstaclesOutside
+template<>
+template<>
+void object::test<12>
+()
+{
+    checkCircle("MULTIPOINT ((10 10), (10 0))",
+        "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+        0.01, 1.0044, 4.997, 10.29 );
+}
 
+// testBoundaryMultiSquares
+template<>
+template<>
+void object::test<13>
+()
+{
+    checkCircle("MULTIPOINT ((10 10), (10 0), (5 5))",
+        "MULTIPOLYGON (((1 9, 9 9, 9 1, 1 1, 1 9)), ((15 20, 20 20, 20 15, 15 15, 15 20)))",
+        0.01, 19.995, 19.997, 14.137 );
+}
+
+// testBoundaryAsObstacle
+template<>
+template<>
+void object::test<14>
+()
+{
+    checkCircle("GEOMETRYCOLLECTION (POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9)), POINT (4 3), POINT (7 6))",
+        "POLYGON ((1 9, 9 9, 9 1, 1 1, 1 9))",
+        0.01, 4, 6, 3 );
+}
 
+} // namespace tut