geos
/
geos
5
1
Fork 16
Code Issues Pull Requests Releases Activity

Fix BufferOp to avoid artifacts in certain polygon buffers
continuous-integration/drone/push Build is passing Details

gh-autolink
mdavis 2021-04-13 11:16:57 -07:00
parent 70b1662e9b
commit 733bddf948
3 changed files with 172 additions and 24 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

37
include/geos/operation/buffer/OffsetCurveSetBuilder.h
View File

@ -139,6 +139,42 @@ private:
double offsetDistance, int side, geom::Location cwLeftLoc,
geom::Location cwRightLoc);
/**
* Tests whether the offset curve for a ring is fully inverted.
* An inverted ("inside-out") curve occurs in some specific situations
* involving a buffer distance which should result in a fully-eroded (empty) buffer.
* It can happen that the sides of a small, convex polygon
* produce offset segments which all cross one another to form
* a curve with inverted orientation.
* This happens at buffer distances slightly greater than the distance at
* which the buffer should disappear.
* The inverted curve will produce an incorrect non-empty buffer (for a shell)
* or an incorrect hole (for a hole).
* It must be discarded from the set of offset curves used in the buffer.
* Heuristics are used to reduce the number of cases which area checked,
* for efficiency and correctness.
* <p>
* See https://github.com/locationtech/jts/issues/472
*
* @param inputPts the input ring
* @param distance the buffer distance
* @param curvePts the generated offset curve
* @return true if the offset curve is inverted
*/
static bool isRingCurveInverted(
const geom::CoordinateSequence* inputPts, double dist,
const geom::CoordinateSequence* curvePts);
/**
* Computes the maximum distance out of a set of points to a linestring.
*
* @param pts the points
* @param line the linestring vertices
* @return the maximum distance
*/
static double maxDistance(
const geom::CoordinateSequence* pts, const geom::CoordinateSequence* line);
/**
* The ringCoord is assumed to contain no repeated points.
* It may be degenerate (i.e. contain only 1, 2, or 3 points).
@ -216,4 +252,3 @@ public:
#endif
#endif // ndef GEOS_OP_BUFFER_OFFSETCURVESETBUILDER_H

91
src/operation/buffer/OffsetCurveSetBuilder.cpp
View File

@ -332,9 +332,78 @@ OffsetCurveSetBuilder::addRingSide(const CoordinateSequence* coord,
}
std::vector<CoordinateSequence*> lineList;
curveBuilder.getRingCurve(coord, side, offsetDistance, lineList);
// ASSERT: lineList contains exactly 1 curve (this is teh JTS semantics)
if (lineList.size() > 0) {
const CoordinateSequence* curve = lineList[0];
/**
* If the offset curve has inverted completely it will produce
* an unwanted artifact in the result, so skip it.
*/
if (isRingCurveInverted(coord, offsetDistance, curve)) {
return;
}
}
addCurves(lineList, leftLoc, rightLoc);
}
static const int MAX_INVERTED_RING_SIZE = 9;
static const double NEARNESS_FACTOR = 0.99;
/* private static*/
bool
OffsetCurveSetBuilder::isRingCurveInverted(
const CoordinateSequence* inputPts, double dist,
const CoordinateSequence* curvePts)
{
if (dist == 0.0) return false;
/**
* Only proper rings can invert.
*/
if (inputPts->size() <= 3) return false;
/**
* Heuristic based on low chance that a ring with many vertices will invert.
* This low limit ensures this test is fairly efficient.
*/
if (inputPts->size() >= MAX_INVERTED_RING_SIZE) return false;
/**
* An inverted curve has no more points than the input ring.
* This also eliminates concave inputs (which will produce fillet arcs)
*/
if (curvePts->size() > inputPts->size()) return false;
/**
* Check if the curve vertices are all closer to the input ring
* than the buffer distance.
* If so, the curve is NOT a valid buffer curve.
*/
double distTol = NEARNESS_FACTOR * fabs(dist);
double maxDist = maxDistance(curvePts, inputPts);
bool isCurveTooClose = maxDist < distTol;
return isCurveTooClose;
}
/**
* Computes the maximum distance out of a set of points to a linestring.
*
* @param pts the points
* @param line the linestring vertices
* @return the maximum distance
*/
/* private static */
double
OffsetCurveSetBuilder::maxDistance(const CoordinateSequence* pts, const CoordinateSequence* line) {
double maxDistance = 0;
for (size_t i = 0; i < pts->size(); i++) {
Coordinate p = pts->getAt(i);
double dist = Distance::pointToSegmentString(p, line);
if (dist > maxDistance) {
maxDistance = dist;
}
}
return maxDistance;
}
/*private*/
bool
OffsetCurveSetBuilder::isErodedCompletely(const LinearRing* ring,
@ -358,28 +427,6 @@ OffsetCurveSetBuilder::isErodedCompletely(const LinearRing* ring,
if(bufferDistance < 0.0 && 2 * std::abs(bufferDistance) > envMinDimension) {
return true;
}
/*
* The following is a heuristic test to determine whether an
* inside buffer will be eroded completely->
* It is based on the fact that the minimum diameter of the ring
* pointset
* provides an upper bound on the buffer distance which would erode the
* ring->
* If the buffer distance is less than the minimum diameter, the ring
* may still be eroded, but this will be determined by
* a full topological computation->
*
*/
/* MD 7 Feb 2005 - there's an unknown bug in the MD code,
so disable this for now */
#if 0
MinimumDiameter md(ring); //=new MinimumDiameter(ring);
double minDiam = md.getLength();
return minDiam < (2 * std::fabs(bufferDistance));
#endif
return false;
}

68
tests/unit/operation/buffer/BufferOpTest.cpp
View File

@ -11,6 +11,7 @@
#include <geos/geom/Coordinate.h>
#include <geos/geom/GeometryFactory.h>
#include <geos/geom/Geometry.h>
#include <geos/geom/Polygon.h>
#include <geos/algorithm/PointLocator.h>
#include <geos/io/WKTReader.h>
#include <geos/io/WKTWriter.h>
@ -464,6 +465,71 @@ void object::test<15>
ensure_equals_geometry(gresult.get(), gexpected.get());
}
// Test for #1101 - Non-empty negative buffer of 4-pt convex polygon
template<>
template<>
void object::test<16>
()
{
std::string wkt0("POLYGON ((666360.09 429614.71, 666344.4 429597.12, 666358.47 429584.52, 666374.5 429602.33, 666360.09 429614.71))");
GeomPtr g0(wktreader.read(wkt0));
ensure_not( g0->buffer( -9 )->isEmpty() );
ensure( g0->buffer( -10 )->isEmpty() );
ensure( g0->buffer( -15 )->isEmpty() );
ensure( g0->buffer( -18 )->isEmpty() );
}
// Test for #1101 - Non-empty negative buffer of 5-pt convex polygon
template<>
template<>
void object::test<17>
()
{
std::string wkt0("POLYGON ((6 20, 16 20, 21 9, 9 0, 0 10, 6 20))");
GeomPtr g0(wktreader.read(wkt0));
ensure_not( g0->buffer( -8 )->isEmpty() );
ensure( g0->buffer( -8.6 )->isEmpty() );
ensure( g0->buffer( -9.6 )->isEmpty() );
ensure( g0->buffer( -11 )->isEmpty() );
}
// Test for #1101 - Buffer of Polygon with hole with hole eroded
template<>
template<>
void object::test<18>
()
{
std::string wkt0("POLYGON ((-6 26, 29 26, 29 -5, -6 -5, -6 26), (6 20, 16 20, 21 9, 9 0, 0 10, 6 20))");
GeomPtr g0(wktreader.read(wkt0));
GeomPtr result1 = g0->buffer( -8 );
ensure( 0 == dynamic_cast<const geos::geom::Polygon*>(result1.get())->getNumInteriorRing() );
GeomPtr result2 = g0->buffer( -8.6 );
ensure( 0 == dynamic_cast<const geos::geom::Polygon*>(result2.get())->getNumInteriorRing() );
GeomPtr result3 = g0->buffer( -9.6 );
ensure( 0 == dynamic_cast<const geos::geom::Polygon*>(result2.get())->getNumInteriorRing() );
GeomPtr result4 = g0->buffer( -11 );
ensure( 0 == dynamic_cast<const geos::geom::Polygon*>(result2.get())->getNumInteriorRing() );
}
// Test for #1101 - Non-empty negative buffer of 5-pt convex polygon
template<>
template<>
void object::test<19>
()
{
std::string wkt0("MULTIPOLYGON (((30 18, 14 0, 0 13, 16 30, 30 18)), ((180 210, 60 50, 154 6, 270 40, 290 130, 250 190, 180 210)))");
GeomPtr g0(wktreader.read(wkt0));
ensure( 2 == g0->buffer( -9 )->getNumGeometries() );
ensure( 1 == g0->buffer( -10 )->getNumGeometries() );
ensure( 1 == g0->buffer( -15 )->getNumGeometries() );
ensure( 1 == g0->buffer( -18 )->getNumGeometries() );
}
} // namespace tut