Buffer of linestring sometimes becomes multipolygon

[Komzpa]

Based on https://trac.osgeo.org/postgis/ticket/2181

Expected that buffer of linestring is a polygon by @dr-jts in https://lists.osgeo.org/pipermail/postgis-users/2013-January/036342.html.

Examples:

22:25:25 [kom] > SELECT 'ticket2181-linestring', GeometryType(ST_Buffer(                                                                                                                                                     
        'LINESTRING(640770.332537465 216785.425146015,640770.356311913 216785.560172686,640770.475488952 216786.226162019,640770.647396495 216787.800152365,640770.658978918 216789.268376902,640770.447752096 216792.399898366,640769.98018796 216795.624688816,640769.110287021 216798.527305023,640767.648768098 216801.891389348,640766.02448565 216804.805892778,640764.74005588 216807.134569434,640763.337363474 216809.560477487,640761.650901703 216812.246990061,640759.955600712 216815.327305019,640758.482499366 216817.868732504,640757.383693167 216819.548488661,640756.238862077 216821.226720821,640754.01503683 216824.434442434,640751.366929134 216828.435255234,640749.204368809 216832.450088865,640747.25516891 216836.227787618,640746.014020829 216838.801219168,640745.393751588 216840.360274284,640744.315976633 216842.854457673,640743.400050801 216844.65704848,640742.582575565 216846.277673322,640741.52491745 216848.478638524,640740.503835408 216850.705816576,640739.787553975 216852.270967708,640739.053289307 216853.851968469,640738.173025147 216855.546355058,640737.55488951 216856.568961104,640736.817272035 216857.822910812,640735.910795021 216859.35697228,640734.77510795 216861.236982439,640733.832969266 216862.938074642,640732.814325629 216865.039674844,640731.225095251 216869.012141189,640729.979984761 216871.879095724,640729.445974092 216873.02148841,640729.002794006 216873.679857725,640728.952197105 216873.745389857,640728.676962154 216874.089814544)',
        100));
┌───────────────────────┬──────────────┐
│       ?column?        │ geometrytype │
├───────────────────────┼──────────────┤
│ ticket2181-linestring │ MULTIPOLYGON │
└───────────────────────┴──────────────┘
(1 row)

Time: 0,691 ms
22:25:27 [kom] > SELECT ST_GeometryType(ST_AsText(ST_Buffer(ST_GeomFromText('POLYGON((19.4994 56.0304, 19.2074 56.0271, 18.9318 56.0205, 18.8957 57.9268, 20.5001 57.9498, 20.0079 56.0369, 19.7389 56.0304, 19.4994 56.0304))'), 42.0)));
┌─────────────────┐
│ st_geometrytype │
├─────────────────┤
│ ST_MultiPolygon │
└─────────────────┘
(1 row)

Time: 0,351 ms
22:25:25 [kom] > SELECT 'ticket2181-linestring', GeometryType(ST_Buffer(                                                                                                                                                     
        'LINESTRING(640770.332537465 216785.425146015,640770.356311913 216785.560172686,640770.475488952 216786.226162019,640770.647396495 216787.800152365,640770.658978918 216789.268376902,640770.447752096 216792.399898366,640769.98018796 216795.624688816,640769.110287021 216798.527305023,640767.648768098 216801.891389348,640766.02448565 216804.805892778,640764.74005588 216807.134569434,640763.337363474 216809.560477487,640761.650901703 216812.246990061,640759.955600712 216815.327305019,640758.482499366 216817.868732504,640757.383693167 216819.548488661,640756.238862077 216821.226720821,640754.01503683 216824.434442434,640751.366929134 216828.435255234,640749.204368809 216832.450088865,640747.25516891 216836.227787618,640746.014020829 216838.801219168,640745.393751588 216840.360274284,640744.315976633 216842.854457673,640743.400050801 216844.65704848,640742.582575565 216846.277673322,640741.52491745 216848.478638524,640740.503835408 216850.705816576,640739.787553975 216852.270967708,640739.053289307 216853.851968469,640738.173025147 216855.546355058,640737.55488951 216856.568961104,640736.817272035 216857.822910812,640735.910795021 216859.35697228,640734.77510795 216861.236982439,640733.832969266 216862.938074642,640732.814325629 216865.039674844,640731.225095251 216869.012141189,640729.979984761 216871.879095724,640729.445974092 216873.02148841,640729.002794006 216873.679857725,640728.952197105 216873.745389857,640728.676962154 216874.089814544)',
        100));
┌───────────────────────┬──────────────┐
│       ?column?        │ geometrytype │
├───────────────────────┼──────────────┤
│ ticket2181-linestring │ MULTIPOLYGON │
└───────────────────────┴──────────────┘
(1 row)

Time: 0,691 ms
22:25:27 [kom] > SELECT ST_GeometryType(ST_AsText(ST_Buffer(ST_GeomFromText('POLYGON((19.4994 56.0304, 19.2074 56.0271, 18.9318 56.0205, 18.8957 57.9268, 20.5001 57.9498, 20.0079 56.0369, 19.7389 56.0304, 19.4994 56.0304))'), 42.0)));
┌─────────────────┐
│ st_geometrytype │
├─────────────────┤
│ ST_MultiPolygon │
└─────────────────┘
(1 row)

Time: 0,351 ms

[dbaston]

The "extra" part of the MultiPolygon is a little sliver on the upper-right:

[dr-jts]

This occurs sometimes due to a complex interaction between noding and some heuristics used in generating a robust buffer line.

For the case of single-element geometries perhaps a simple heuristic to keep only the largest result polygon would suffice.

[Komzpa]

Letting AI loose on test case produces this after a while. Does this help with actual fix?

diff --git a/include/geos/operation/buffer/OffsetSegmentString.h b/include/geos/operation/buffer/OffsetSegmentString.h
index a3fdf1943b249476ac955afed3912f0c78e9bbbc..29ef229dd7f04c27c03cdea941fc16a7d056397e 100644
--- a/include/geos/operation/buffer/OffsetSegmentString.h
+++ b/include/geos/operation/buffer/OffsetSegmentString.h
@@ -145,50 +145,55 @@ public:
             for(std::size_t i = 0, n = pts.size(); i < n; ++i) {
                 addPt(pts[i]);
             }
         }
         else {
             for(std::size_t i = pts.size(); i > 0; --i) {
                 addPt(pts[i - 1]);
             }
         }
     }
 
     /// Check that points are a ring
     ///
     /// add the startpoint again if they are not
     void
     closeRing()
     {
         if(ptList->size() < 1) {
             return;
         }
         const geom::Coordinate& startPt = ptList->front();
         const geom::Coordinate& lastPt = ptList->back();
         if(startPt.equals(lastPt)) {
             return;
         }
+        if(startPt.distance(lastPt) < minimumVertexDistance) {
+            // drop the nearly coincident closing point to avoid creating
+            // tiny spikes at the seam
+            ptList->pop_back();
+        }
         // we ask to allow repeated as we checked this ourself
         ptList->add(startPt, true);
     }
 
     /// Get coordinates by taking ownership of them
     ///
     /// After this call, the coordinates reference in
     /// this object are dropped. Calling twice will
     /// segfault...
     ///
     /// FIXME: refactor memory management of this
     ///
     geom::CoordinateSequence*
     getCoordinates()
     {
         closeRing();
         geom::CoordinateSequence* ret = ptList;
         ptList = nullptr;
         return ret;
     }
 
     inline size_t
     size() const
     {
         return ptList ? ptList->size() : 0 ;
diff --git a/src/operation/buffer/OffsetCurveBuilder.cpp b/src/operation/buffer/OffsetCurveBuilder.cpp
index 1628ef9bb2d3dd154a23b558f8761f267747b835..d9d7a655b6423f9d5a158b6e62016526758123ce 100644
--- a/src/operation/buffer/OffsetCurveBuilder.cpp
+++ b/src/operation/buffer/OffsetCurveBuilder.cpp
@@ -332,51 +332,59 @@ OffsetCurveBuilder::getRingCurve(const CoordinateSequence* inputPts,
 }
 
 
 /* public */
 std::unique_ptr<CoordinateSequence>
 OffsetCurveBuilder::getRingCurve(const CoordinateSequence* inputPts, int side, double pDistance)
 {
     distance = pDistance;
     if (inputPts->size() <= 2)
         return getLineCurve(inputPts, distance);
 
     // optimize creating ring for for zero distance
     if (distance == 0.0) {
         return inputPts->clone();
     }
     OffsetSegmentGenerator segGen(precisionModel, bufParams, distance);
     computeRingBufferCurve(*inputPts, side, segGen);
     return segGen.getCoordinates();
 }
 
 
 /* private */
 double
 OffsetCurveBuilder::simplifyTolerance(double bufDistance)
 {
-    return bufDistance / SIMPLIFY_FACTOR;
+    double distTol = bufDistance / SIMPLIFY_FACTOR;
+    if(bufDistance > 1.0) {
+        distTol = 1.0 / SIMPLIFY_FACTOR;
+    }
+    int quadSegs = bufParams.getQuadrantSegments();
+    if(quadSegs > 1) {
+        distTol /= quadSegs;
+    }
+    return distTol;
 }
 
 /*private*/
 void
 OffsetCurveBuilder::computeLineBufferCurve(const CoordinateSequence& inputPts,
         OffsetSegmentGenerator& segGen)
 {
     double distTol = simplifyTolerance(distance);
 
     //--------- compute points for left side of line
     // Simplify the appropriate side of the line before generating
     std::unique_ptr<CoordinateSequence> simp1_ =
         BufferInputLineSimplifier::simplify(inputPts, distTol);
     const CoordinateSequence& simp1 = *simp1_;
 
 
     auto n1 = simp1.size() - 1;
     segGen.initSideSegments(simp1[0], simp1[1], Position::LEFT);
     for(std::size_t i = 2; i <= n1; ++i) {
         segGen.addNextSegment(simp1[i], true);
     }
     segGen.addLastSegment();
     // add line cap for end of line
     segGen.addLineEndCap(simp1[n1 - 1], simp1[n1]);
 
diff --git a/src/operation/buffer/OffsetSegmentGenerator.cpp b/src/operation/buffer/OffsetSegmentGenerator.cpp
index 6c7c2c3b088903f87375539731af3671cba1a6e4..f35fdb3a685aa0b0b3b6079eff2b8c65d1389feb 100644
--- a/src/operation/buffer/OffsetSegmentGenerator.cpp
+++ b/src/operation/buffer/OffsetSegmentGenerator.cpp
@@ -22,51 +22,51 @@
 
 #include <geos/algorithm/Angle.h>
 #include <geos/algorithm/Distance.h>
 #include <geos/algorithm/Orientation.h>
 #include <geos/operation/buffer/OffsetSegmentGenerator.h>
 #include <geos/operation/buffer/OffsetSegmentString.h>
 #include <geos/operation/buffer/BufferInputLineSimplifier.h>
 #include <geos/operation/buffer/BufferParameters.h>
 #include <geos/geom/Position.h>
 #include <geos/geom/CoordinateSequence.h>
 #include <geos/geom/Coordinate.h>
 #include <geos/geom/PrecisionModel.h>
 #include <geos/algorithm/NotRepresentableException.h>
 #include <geos/algorithm/CGAlgorithmsDD.h>
 #include <geos/util.h>
 
 
 using namespace geos::algorithm;
 using namespace geos::geom;
 
 namespace geos {
 namespace operation { // geos.operation
 namespace buffer { // geos.operation.buffer
 
 /*private data*/
-const double OffsetSegmentGenerator::CURVE_VERTEX_SNAP_DISTANCE_FACTOR = 1.0E-4;
+const double OffsetSegmentGenerator::CURVE_VERTEX_SNAP_DISTANCE_FACTOR = 5.0E-4;
 const double OffsetSegmentGenerator::OFFSET_SEGMENT_SEPARATION_FACTOR = 1.0E-3;
 const double OffsetSegmentGenerator::INSIDE_TURN_VERTEX_SNAP_DISTANCE_FACTOR = 1.0E-3;
 const double OffsetSegmentGenerator::SIMPLIFY_FACTOR = 100.0;
 
 /*public*/
 OffsetSegmentGenerator::OffsetSegmentGenerator(
     const PrecisionModel* newPrecisionModel,
     const BufferParameters& nBufParams,
     double dist)
     :
     maxCurveSegmentError(0.0),
     closingSegLengthFactor(1),
     segList(),
     distance(dist),
     precisionModel(newPrecisionModel),
     bufParams(nBufParams),
     li(),
     s0(),
     s1(),
     s2(),
     seg0(),
     seg1(),
     offset0(),
     offset1(),
     side(0),
diff --git a/tests/unit/operation/buffer/BufferOpTest.cpp b/tests/unit/operation/buffer/BufferOpTest.cpp
index b51914ca8ddc1dacbdf7c92f24e074b64901c074..ba0dba758d716f17f7ed8b957399ebac9334c3a9 100644
--- a/tests/unit/operation/buffer/BufferOpTest.cpp
+++ b/tests/unit/operation/buffer/BufferOpTest.cpp
@@ -598,29 +598,75 @@ void object::test<26>
 }
 
 // testElementErodedEx
 // Checks that a skinny element polygon is eroded with no internal predicision reduction due to topo exes
 // see https://github.com/libgeos/geos/issues/1182
 template<>
 template<>
 void object::test<27>
 ()
 {
     std::string wkt = "MULTIPOLYGON (((48268.99938 -49048.29324, 44429.1 -55700.232847, 44429.1 -55107.317582, 44506.1 -54974, 44429.1 -54840, 44429.1 -51569.2, 42170.10515 -49316.27944, 48268.99938 -49048.29324)), ((43433.08324 -51823.15037, 42480.09977 -55494.96132, 42477.638798 -55504.400121, 42480.20715 -55494.547587, 42482.247919 -55485.931009, 42482.431666 -55485.976608, 43433.08324 -51823.15037)))";
     checkBuffer(wkt, -1, 0.01,
 "POLYGON ((48267.2218198241 -49049.37231112561, 44430.1 -55696.500291383236, 44430.1 -55107.58561879013, 44506.96594366424 -54974.50014155032, 44507.05088958367 -54974.309530288774, 44507.09442572395 -54974.10543945913, 44507.09465615313 -54973.896756903174, 44507.05157083636 -54973.69257042533, 44506.96704607295 -54973.50177203271, 44430.1 -54839.73314639927, 44430.1 -51569.2, 44430.08071956347 -51569.00457958696, 44430.02362172434 -51568.81669475566, 44429.93090822513 -51568.64359050924, 44429.80615417933 -51568.49194189856, 42172.42282165639 -49317.178566110095, 48267.2218198241 -49049.37231112561))");    
 }
 
 // testBufferByZeroKeepsAllElements
 // Test a buffer-by-zero case which revealed a bug in the DeptSegment comparator.
 // See https://github.com/locationtech/jts/issues/1131
 template<>
 template<>
 void object::test<28>
 ()
 {
     std::string wkt = "MULTIPOLYGON (((24 95.239, 24 96, 24 99, 24.816 99, 24 95.239)), ((3 90, 3 93, 3 96, 3 99, 21 99, 21 96, 21 93, 21 90, 3 90)))";
     checkBuffer(wkt, 0, 0.01,
-"MULTIPOLYGON (((24 95.239, 24 96, 24 99, 24.816 99, 24 95.239)), ((3 90, 3 93, 3 96, 3 99, 21 99, 21 96, 21 93, 21 90, 3 90)))");    
+"MULTIPOLYGON (((24 95.239, 24 96, 24 99, 24.816 99, 24 95.239)), ((3 90, 3 93, 3 96, 3 99, 21 99, 21 96, 21 93, 21 90, 3 90)))");
+}
+
+// Buffer of complex linestring should be a single polygon without slivers
+template<>
+template<>
+void object::test<29>
+()
+{
+    const std::string wktLine =
+        "LINESTRING(640770.332537465 216785.425146015,640770.356311913 216785.560172686,640770.475488952 216786.226162019,"
+        "640770.647396495 216787.800152365,640770.658978918 216789.268376902,640770.447752096 216792.399898366,640769.98018796 "
+        "216795.624688816,640769.110287021 216798.527305023,640767.648768098 216801.891389348,640766.02448565 216804.805892778,"
+        "640764.74005588 216807.134569434,640763.337363474 216809.560477487,640761.650901703 216812.246990061,640759.955600712 "
+        "216815.327305019,640758.482499366 216817.868732504,640757.383693167 216819.548488661,640756.238862077 216821.226720821,"
+        "640754.01503683 216824.434442434,640751.366929134 216828.435255234,640749.204368809 216832.450088865,640747.25516891 "
+        "216836.227787618,640746.014020829 216838.801219168,640745.393751588 216840.360274284,640744.315976633 216842.854457673,"
+        "640743.400050801 216844.65704848,640742.582575565 216846.277673322,640741.52491745 216848.478638524,640740.503835408 "
+        "216850.705816576,640739.787553975 216852.270967708,640739.053289307 216853.851968469,640738.173025147 216855.546355058,"
+        "640737.55488951 216856.568961104,640736.817272035 216857.822910812,640735.910795021 216859.35697228,640734.77510795 "
+        "216861.236982439,640733.832969266 216862.938074642,640732.814325629 216865.039674844,640731.225095251 216869.012141189,"
+        "640729.979984761 216871.879095724,640729.445974092 216873.02148841,640729.002794006 216873.679857725,640728.952197105 "
+        "216873.745389857,640728.676962154 216874.089814544)";
+
+    auto buffered = buffer(wktLine, 100.0);
+    ensure_equals(buffered->getGeometryTypeId(), geos::geom::GEOS_POLYGON);
+    ensure_equals(buffered->getNumGeometries(), std::size_t(1));
+    ensure_equals(buffered->getArea(), 51020.066327663415, 1e-3);
+    ensure_equals(buffered->getLength(), 825.03211479112, 1e-6);
+}
+
+// Buffer of polygon should remain polygon with expected area/perimeter
+template<>
+template<>
+void object::test<30>
+()
+{
+    const std::string wktPoly =
+        "POLYGON((19.4994 56.0304,19.2074 56.0271,18.9318 56.0205,18.8957 57.9268,20.5001 57.9498,20.0079 56.0369,19.7389 56.0304,"
+        "19.4994 56.0304))";
+
+    auto gPoly = wktreader.read(wktPoly);
+    auto buffered = BufferOp::bufferOp(gPoly.get(), 42.0);
+    ensure_equals(buffered->getGeometryTypeId(), geos::geom::GEOS_POLYGON);
+    ensure_equals(buffered->getNumGeometries(), std::size_t(1));
+    ensure_equals(buffered->getArea(), 5784.289705520652, 1e-3);
+    ensure_equals(buffered->getLength(), 270.0312290469457, 1e-6);
 }
 
 } // namespace tut

[dr-jts]

@dbaston see locationtech/jts#1161 for a slightly better code location for this fix.