Fixed MinimumBoundingCircle.getMaximumDiameter() to handle case where…

… extremalPts[2]-[0] is longest, JTS 533
klassenjs · Apr 28, 2020 · cba677c · cba677c
1 parent 46fcfd5
commit cba677c
Show file tree

Hide file tree

Showing 4 changed files with 91 additions and 20 deletions.
diff --git a/NEWS b/NEWS
@@ -15,6 +15,7 @@ Changes in 3.9.0
   - Fix bug in DistanceOp for geometries with empty components (#1026, Paul Ramsey)
   - Remove undefined behaviour in CAPI (#1021, Greg Troxel)
   - Fix buffering issue (#1022, JTS-525, Paul Ramsey)
+  - MinimumBoundingCircle.getMaximumDiameter fix (JTS-533, Paul Ramsey)
 
 
 Changes in 3.8.0

diff --git a/include/geos/algorithm/MinimumBoundingCircle.h b/include/geos/algorithm/MinimumBoundingCircle.h
@@ -56,6 +56,7 @@ class GEOS_DLL MinimumBoundingCircle {
     geom::Coordinate pointWitMinAngleWithX(std::vector<geom::Coordinate>& pts, geom::Coordinate& P);
     geom::Coordinate pointWithMinAngleWithSegment(std::vector<geom::Coordinate>& pts,
             geom::Coordinate& P, geom::Coordinate& Q);
+    std::vector<geom::Coordinate> farthestPoints(std::vector<geom::Coordinate>& pts);
 
 
 public:
@@ -93,7 +94,7 @@ class GEOS_DLL MinimumBoundingCircle {
     * @return a empty LineString if the input is empty
     * @return a Point if the input is a point
     */
-    std::unique_ptr<geom::Geometry> getFarthestPoints();
+    std::unique_ptr<geom::Geometry> getMaximumDiameter();
 
     /**
     * Gets a geometry representing the diameter of the computed Minimum Bounding

diff --git a/src/algorithm/MinimumBoundingCircle.cpp b/src/algorithm/MinimumBoundingCircle.cpp
@@ -61,22 +61,54 @@ MinimumBoundingCircle::getCircle()
 
 /*public*/
 std::unique_ptr<Geometry>
-MinimumBoundingCircle::getFarthestPoints()
+MinimumBoundingCircle::getMaximumDiameter()
 {
     compute();
+    size_t dims = input->getCoordinateDimension();
+    size_t len = 2;
     switch(extremalPts.size()) {
-    case 0:
-        return input->getFactory()->createLineString();
-    case 1:
-        return std::unique_ptr<Geometry>(input->getFactory()->createPoint(centre));
+        case 0:
+            return input->getFactory()->createLineString();
+        case 1:
+            return std::unique_ptr<Geometry>(input->getFactory()->createPoint(centre));
+        case 2: {
+            auto cs = input->getFactory()->getCoordinateSequenceFactory()->create(len, dims);
+            cs->setAt(extremalPts.front(), 0);
+            cs->setAt(extremalPts.back(), 1);
+            return input->getFactory()->createLineString(std::move(cs));
+        }
+        default: {
+            std::vector<Coordinate> fp = farthestPoints(extremalPts);
+            auto cs = input->getFactory()->getCoordinateSequenceFactory()->create(len, dims);
+            cs->setAt(fp.front(), 0);
+            cs->setAt(fp.back(), 1);
+            return input->getFactory()->createLineString(std::move(cs));
+        }
     }
 
-    size_t dims = input->getCoordinateDimension();
-    size_t len = 2;
-    auto cs = input->getFactory()->getCoordinateSequenceFactory()->create(len, dims);
-    cs->setAt(extremalPts[0], 0);
-    cs->setAt(extremalPts[extremalPts.size() - 1], 1);
-    return input->getFactory()->createLineString(std::move(cs));
+}
+
+/* private */
+std::vector<Coordinate>
+MinimumBoundingCircle::farthestPoints(std::vector<Coordinate>& pts)
+{
+    std::vector<Coordinate> fp;
+    double dist01 = pts[0].distance(pts[1]);
+    double dist12 = pts[1].distance(pts[2]);
+    double dist20 = pts[2].distance(pts[0]);
+    if (dist01 >= dist12 && dist01 >= dist20) {
+        fp.push_back(pts[0]);
+        fp.push_back(pts[1]);
+        return fp;
+    }
+    if (dist12 >= dist01 && dist12 >= dist20) {
+        fp.push_back(pts[1]);
+        fp.push_back(pts[2]);
+        return fp;
+    }
+    fp.push_back(pts[2]);
+    fp.push_back(pts[0]);
+    return fp;
 }
 
 /*public*/
@@ -100,6 +132,7 @@ MinimumBoundingCircle::getDiameter()
     return input->getFactory()->createLineString(std::move(cs));
 }
 
+
 /*public*/
 std::vector<Coordinate>
 MinimumBoundingCircle::getExtremalPoints()

diff --git a/tests/unit/algorithm/MinimumBoundingCircleTest.cpp b/tests/unit/algorithm/MinimumBoundingCircleTest.cpp
@@ -76,11 +76,15 @@ struct test_minimumboundingcircle_data {
         }
         ensure(isEqual);
 
-        if(centreOut.isNull()) {
-            ensure(centreOut.distance(actualCentre) < 0.0001);
+        if(!centreOut.isNull()) {
+            if (centreOut.distance(actualCentre) > 0.001) {
+                std::cout << "centreOut " << centreOut << std::endl;
+                std::cout << "actualCentre " << actualCentre << std::endl;
+            }
+            ensure_equals("centerOut", centreOut.distance(actualCentre), 0.0, 0.001);
         }
         if(radiusOut >= 0) {
-            ensure(fabs(radiusOut - actualRadius) < 0.0001);
+            ensure_equals("radius", actualRadius, radiusOut, 0.0001);
         }
     }
 
@@ -100,15 +104,10 @@ typedef group::object object;
 group test_minimumboundingcircle_group("geos::algorithm::MinimumBoundingCircle");
 
 
-
-
-
-
 //
 // Test Cases
 //
 
-
 template<>
 template<>
 void object::test<1>
@@ -201,7 +200,44 @@ void object::test<7>
         247.4360455914027);
 }
 
+template<>
+template<>
+void object::test<8>
+()
+{
+    Coordinate c(196.026, 159.103);
+    doMinimumBoundingCircleTest(
+        "POLYGON ((100 200, 300 150, 110 100, 100 200))",
+        "MULTIPOINT ((110 100), (300 150), (100 200))",
+        c,
+        104.372);
+}
 
+template<>
+template<>
+void object::test<9>
+()
+{
+    Coordinate c(196.026, 140.897);
+    doMinimumBoundingCircleTest(
+        "POLYGON ((110 200, 300 150, 100 100, 110 200))",
+        "MULTIPOINT ((100 100), (300 150), (110 200))",
+        c,
+        104.37204);
+}
+
+template<>
+template<>
+void object::test<10>
+()
+{
+    Coordinate c(3, 2);
+    doMinimumBoundingCircleTest(
+        "POLYGON ((0 0, 6 0, 5 5, 0 0))",
+        "MULTIPOINT ((0 0), (6 0), (5 5))",
+        c,
+        3.60555);
+}
 
 } // namespace tut