Fix MinimumDiameter getMinimumRectangle for flat input

include/geos/algorithm/MinimumDiameter.h

@@ -28,6 +28,7 @@
 // Forward declarations
 namespace geos {
 namespace geom {
+class GeometryFactory;
 class Geometry;
 class LineString;
 class CoordinateSequence;
@@ -92,6 +93,10 @@ class GEOS_DLL MinimumDiameter {
 
     static geom::LineSegment computeSegmentForLine(double a, double b, double c);
 
+    static std::unique_ptr<geom::Geometry> computeMaximumLine(
+                        const geom::CoordinateSequence* pts,
+                        const geom::GeometryFactory* factory);
+
 public:
     ~MinimumDiameter() = default;
 
@@ -174,4 +179,3 @@ class GEOS_DLL MinimumDiameter {
 
 } // namespace geos::algorithm
 } // namespace geos
-

src/algorithm/MinimumDiameter.cpp

@@ -240,6 +240,8 @@ MinimumDiameter::findMaxPerpDistance(const CoordinateSequence* pts,
         maxPerpDistance = nextPerpDistance;
         maxIndex = nextIndex;
         nextIndex = getNextIndex(pts, maxIndex);
+        if (nextIndex == startIndex)
+            break;
         nextPerpDistance = seg->distancePerpendicular(pts->getAt(nextIndex));
     }
 
@@ -278,7 +280,8 @@ MinimumDiameter::getMinimumRectangle()
         if(minBaseSeg.p0.equals2D(minBaseSeg.p1)) {
             return std::unique_ptr<Geometry>(inputGeom->getFactory()->createPoint(minBaseSeg.p0));
         }
-        return minBaseSeg.toGeometry(*inputGeom->getFactory());
+      //-- Min rectangle is a line. Use the diagonal of the extent
+      return computeMaximumLine(convexHullPts.get(), inputGeom->getFactory());
     }
 
     // deltas for the base segment of the minimum diameter
@@ -337,6 +340,42 @@ MinimumDiameter::getMinimumRectangle()
     return inputGeom->getFactory()->createPolygon(std::move(shell));
 }
 
+// private static
+std::unique_ptr<Geometry>
+MinimumDiameter::computeMaximumLine(const geom::CoordinateSequence* pts,
+        const GeometryFactory* factory)
+{
+    //-- find max and min pts for X and Y
+    Coordinate ptMinX = pts->getAt(0);
+    Coordinate ptMaxX = pts->getAt(0);
+    Coordinate ptMinY = pts->getAt(0);
+    Coordinate ptMaxY = pts->getAt(0);
+
+    std::size_t const n = pts->getSize();
+    for(std::size_t i = 1; i < n; ++i) {
+      const Coordinate& p = pts->getAt(i);
+      if (p.x < ptMinX.x) ptMinX = p;
+      if (p.x > ptMaxX.x) ptMaxX = p;
+      if (p.y < ptMinY.y) ptMinY = p;
+      if (p.y > ptMaxY.y) ptMaxY = p;
+    }
+    Coordinate p0 = ptMinX;
+    Coordinate p1 = ptMaxX;
+    //-- line is vertical - use Y pts
+    if (p0.x == p1.x) {
+      p0 = ptMinY;
+      p1 = ptMaxY;
+    }
+
+    const CoordinateSequenceFactory* csf =
+        factory->getCoordinateSequenceFactory();
+    auto seq = csf->create(2, 2);
+    seq->setAt(p0, 0);
+    seq->setAt(p1, 1);
+
+    return factory->createLineString(std::move(seq));
+}
+
 double
 MinimumDiameter::computeC(double a, double b, const Coordinate& p)
 {
@@ -381,4 +420,3 @@ MinimumDiameter::getMinimumDiameter(Geometry* geom)
 
 } // namespace geos.algorithm
 } // namespace geos
-

tests/unit/capi/GEOSMinimumRotatedRectangleTest.cpp

@@ -0,0 +1,120 @@
+//
+// Test Suite for C-API GEOSMinimumRotatedRectangle
+
+#include <tut/tut.hpp>
+// geos
+#include <geos_c.h>
+// std
+#include <cstdarg>
+#include <cstdio>
+#include <cstdlib>
+
+#include "capi_test_utils.h"
+
+namespace tut {
+//
+// Test Group
+//
+
+// Common data used in test cases.
+struct test_minimumrotatedrectangle_data : public capitest::utility {
+    test_minimumrotatedrectangle_data() {
+        GEOSWKTWriter_setTrim(wktw_, 1);
+        GEOSWKTWriter_setRoundingPrecision(wktw_, 8);
+    }
+
+    void checkMinRectangle(const char* wkt, const char* expected)
+    {
+        // input
+        geom1_ = GEOSGeomFromWKT(wkt);
+        ensure(nullptr != geom1_);
+
+        // result
+        geom2_ = GEOSMinimumRotatedRectangle(geom1_);
+        ensure(nullptr != geom2_);
+
+        // expected
+        if (expected) {
+            geom3_ = GEOSGeomFromWKT(expected);
+            ensure(nullptr != geom3_);
+            ensure_geometry_equals(geom2_, geom3_, 0.0001);
+        }
+    }
+};
+
+typedef test_group<test_minimumrotatedrectangle_data> group;
+typedef group::object object;
+
+group test_capigeosminimumrotatedrectangle_group("capi::GEOSMinimumRotatedRectangle");
+
+//
+// Test Cases
+//
+
+template<>
+template<>
+void object::test<1>
+()
+{
+    input_ = GEOSGeomFromWKT("POLYGON ((1 6, 6 11, 11 6, 6 1, 1 6))");
+    ensure(nullptr != input_);
+
+    GEOSGeometry* output = GEOSMinimumRotatedRectangle(input_);
+    ensure(nullptr != output);
+    ensure(0 == GEOSisEmpty(output));
+
+    wkt_ = GEOSWKTWriter_write(wktw_, output);
+    ensure_equals(std::string(wkt_), std::string("POLYGON ((6 1, 11 6, 6 11, 1 6, 6 1))"));
+
+    GEOSGeom_destroy(output);
+}
+
+// zero-length
+template<>
+template<>
+void object::test<2>
+()
+{
+    checkMinRectangle("LINESTRING (1 1, 1 1)", "POINT (1 1)");
+}
+
+// Horizontal
+template<>
+template<>
+void object::test<3>
+()
+{
+    checkMinRectangle("LINESTRING (1 1, 3 1, 5 1, 7 1)", "LINESTRING (1 1, 7 1)");
+}
+
+// Vertical
+template<>
+template<>
+void object::test<4>
+()
+{
+    checkMinRectangle("LINESTRING (1 1, 1 4, 1 7, 1 9)", "LINESTRING (1 1, 1 9)");
+}
+
+// Bent Line
+template<>
+template<>
+void object::test<5>
+()
+{
+    checkMinRectangle("LINESTRING (1 2, 3 8, 9 6)", "POLYGON ((9 6, 7 10, -1 6, 1 2, 9 6))");
+}
+
+// Failure case from https://trac.osgeo.org/postgis/ticket/5163
+template<>
+template<>
+void object::test<6>
+()
+{
+    checkMinRectangle(
+        "LINESTRING(-99.48710639268086 34.79029839231914,-99.48370699999998 34.78689899963806,-99.48152167568102 34.784713675318976)",
+        "LINESTRING (-99.48710639268086 34.79029839231914, -99.48152167568102 34.784713675318976)"
+        );
+}
+
+} // namespace tut