Skip to content

Commit

Permalink
Fix buffer Inverted Ring Removal check (#1056)
Browse files Browse the repository at this point in the history
  • Loading branch information
@dr-jts
dr-jts authored Mar 19, 2024
1 parent bc8d092 commit 6c1a7ce
Show file tree
Hide file tree
Showing 5 changed files with 146 additions and 62 deletions.
13 changes: 9 additions & 4 deletions include/geos/geom/LineSegment.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -250,13 +250,18 @@ class GEOS_DLL LineSegment {
//
/// @param ret will be set to the midpoint of the segment
///
void midPoint(Coordinate& ret) const
CoordinateXY midPoint() const
{
ret = Coordinate(
(p0.x + p1.x) / 2,
(p0.y + p1.y) / 2);
return midPoint(p0, p1);
};

static CoordinateXY midPoint(const CoordinateXY& pt0, const CoordinateXY& pt1)
{
return CoordinateXY(
(pt0.x + pt1.x) / 2,
(pt0.y + pt1.y) / 2);
}

/// Computes the distance between this line segment and another one.
double distance(const LineSegment& ls) const
{
Expand Down
19 changes: 12 additions & 7 deletions include/geos/operation/buffer/BufferCurveSetBuilder.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -169,14 +169,19 @@ class GEOS_DLL BufferCurveSetBuilder {
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
* Tests if there are points on the raw offset curve which may
* lie on the final buffer curve
* (i.e. they are (approximately) at the buffer distance from the input ring).
* For efficiency this only tests a limited set of points on the curve.
*
* @param inputRing
* @param distance
* @param curveRing
* @return true if the curve contains points lying at the required buffer distance
*/
static double maxDistance(
const geom::CoordinateSequence* pts, const geom::CoordinateSequence* line);
static bool hasPointOnBuffer(
const CoordinateSequence* inputRing, double dist,
const CoordinateSequence* curveRing);

/**
* The ringCoord is assumed to contain no repeated points.
Expand Down
71 changes: 42 additions & 29 deletions src/operation/buffer/BufferCurveSetBuilder.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,10 +38,12 @@
#include <geos/geomgraph/Label.h>
#include <geos/noding/NodedSegmentString.h>
#include <geos/util.h>
#include <geos/io/WKTWriter.h>

#include <algorithm> // for min
#include <cmath>
#include <cassert>
#include <iomanip>
#include <memory>
#include <vector>
#include <typeinfo>
Expand Down Expand Up @@ -357,56 +359,67 @@ BufferCurveSetBuilder::addRingSide(const CoordinateSequence* coord,
/* private static*/
bool
BufferCurveSetBuilder::isRingCurveInverted(
const CoordinateSequence* inputPts, double dist,
const CoordinateSequence* curvePts)
const CoordinateSequence* inputRing, double dist,
const CoordinateSequence* curveRing)
{
if (dist == 0.0) return false;
/**
* Only proper rings can invert.
*/
if (inputPts->size() <= 3) return false;
if (inputRing->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;
if (inputRing->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() > INVERTED_CURVE_VERTEX_FACTOR * inputPts->size()) return false;
if (curveRing->size() > INVERTED_CURVE_VERTEX_FACTOR * inputRing->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.
* If curve contains points which are on the buffer,
* it is not inverted and can be included in the raw curves.
*/
double distTol = NEARNESS_FACTOR * fabs(dist);
double maxDist = maxDistance(curvePts, inputPts);
bool isCurveTooClose = maxDist < distTol;
return isCurveTooClose;
if (hasPointOnBuffer(inputRing, dist, curveRing))
return false;

//-- curve is inverted, so discard it
return true;
//std::cout << std::setprecision(10) << io::WKTWriter::toLineString(*curveRing) << std::endl;
//std::cout << "isRingCurveInverted: " << isCurveTooClose << " maxDist = " << maxDist << std::endl;
}

/**
* 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
BufferCurveSetBuilder::maxDistance(const CoordinateSequence* pts, const CoordinateSequence* line) {
double maxDistance = 0;
for (std::size_t i = 0; i < pts->size(); i++) {
const Coordinate& p = pts->getAt(i);
double dist = Distance::pointToSegmentString(p, line);
if (dist > maxDistance) {
maxDistance = dist;
/* private static*/
bool
BufferCurveSetBuilder::hasPointOnBuffer(
const CoordinateSequence* inputRing, double dist,
const CoordinateSequence* curveRing)
{
double distTol = NEARNESS_FACTOR * fabs(dist);

for (std::size_t i = 0; i < curveRing->size(); i++) {
const CoordinateXY& v = curveRing->getAt(i);

//-- check curve vertices
double distVertex = Distance::pointToSegmentString(v, inputRing);
if (distVertex > distTol) {
return true;
}

//-- check curve segment midpoints
std::size_t iNext = (i < curveRing->size() - 1) ? i + 1 : 0;
const CoordinateXY& vnext = curveRing->getAt(iNext);
CoordinateXY midPt = LineSegment::midPoint(v, vnext);

double distMid = Distance::pointToSegmentString(midPt, inputRing);
if (distMid > distTol) {
return true;
}
}
return maxDistance;
return false;
}

/*private*/
Expand Down
28 changes: 25 additions & 3 deletions tests/unit/geom/LineSegmentTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,16 +11,17 @@
// std
#include <iostream>

using geos::geom::Coordinate;
using geos::geom::CoordinateXY;
using geos::geom::LineSegment;

namespace tut {
//
// Test Group
//

struct test_lineseg_data {

typedef geos::geom::Coordinate Coordinate;
typedef geos::geom::LineSegment LineSegment;

geos::geom::Coordinate ph1;
geos::geom::Coordinate ph2;
geos::geom::Coordinate pv1;
Expand Down Expand Up @@ -119,6 +120,17 @@ struct test_lineseg_data {
ensure_equals("checkDistancePerpendicularOriented", expected, dist, 0.000001);
}

void checkMidPoint(
double x0, double y0,
double x1, double y1,
double px, double py)
{
LineSegment seg(x0, y0, x1, y1);
Coordinate expected(px, py);
Coordinate actual = Coordinate(seg.midPoint());
ensure_equals_xy(actual, expected);
}

test_lineseg_data()
: ph1(0, 2), ph2(10, 2), pv1(0, 0), pv2(0, 10), h1(ph1, ph2), v1(pv1, pv2)
{}
Expand Down Expand Up @@ -331,6 +343,16 @@ void object::test<13>()
checkDistancePerpendicularOriented(1,1, 1,1, 1,2, 1);
}

// test midpoint
template<>
template<>
void object::test<14>()
{
//-- right of line
checkMidPoint(1,1, 1,3, 1,2);
checkMidPoint(1,1, 1,1, 1,1);
checkMidPoint(1,1, 5,5, 3,3);
}


} // namespace tut
77 changes: 58 additions & 19 deletions tests/unit/operation/buffer/BufferOpTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -43,6 +43,22 @@ struct test_bufferop_data {
{
ensure_equals(default_quadrant_segments, int(8));
}

void checkBufferEmpty(const std::string& wkt, double dist, bool isEmpty)
{
std::unique_ptr<Geometry> geom = wktreader.read(wkt);
std::unique_ptr<Geometry> actual = geom->buffer(dist);
ensure_equals(actual->isEmpty(), isEmpty);
}

void checkBuffer(const std::string& wkt, double dist, double tolerance, const std::string& wktExpected)
{
std::unique_ptr<Geometry> geom = wktreader.read(wkt);
std::unique_ptr<Geometry> actual = geom->buffer(dist);
std::unique_ptr<Geometry> expected = wktreader.read(wktExpected);
ensure_equals_geometry(expected.get(), actual.get(), tolerance);
}

private:
// noncopyable
test_bufferop_data(test_bufferop_data const& other) = delete;
Expand Down Expand Up @@ -482,37 +498,33 @@ void object::test<15>
ensure_equals_geometry(gresult.get(), gexpected.get());
}

// Test for #1101 - Non-empty negative buffer of 4-pt convex polygon
// Test for https://trac.osgeo.org/geos/ticket/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( GeomPtr(g0->buffer( -9 ))->isEmpty() );
ensure( GeomPtr(g0->buffer( -10 ))->isEmpty() );
ensure( GeomPtr(g0->buffer( -15 ))->isEmpty() );
ensure( GeomPtr(g0->buffer( -18 ))->isEmpty() );
checkBufferEmpty(wkt0, -9, false);
checkBufferEmpty(wkt0, -10, true);
checkBufferEmpty(wkt0, -15, true);
checkBufferEmpty(wkt0, -18, true);
}

// Test for #1101 - Non-empty negative buffer of 5-pt convex polygon
// Test for https://trac.osgeo.org/geos/ticket/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( GeomPtr(g0->buffer( -8 ))->isEmpty() );
ensure( GeomPtr(g0->buffer( -8.6 ))->isEmpty() );
ensure( GeomPtr(g0->buffer( -9.6 ))->isEmpty() );
ensure( GeomPtr(g0->buffer( -11 ))->isEmpty() );
checkBufferEmpty(wkt0, -8, false);
checkBufferEmpty(wkt0, -8.6, true);
checkBufferEmpty(wkt0, -9.6, true);
checkBufferEmpty(wkt0, -11, true);
}

// Test for #1101 - Buffer of Polygon with hole with hole eroded
// Test for https://trac.osgeo.org/geos/ticket/1101 - Buffer of Polygon with hole with hole eroded
template<>
template<>
void object::test<18>
Expand All @@ -534,7 +546,7 @@ void object::test<18>
ensure( 0 == dynamic_cast<const geos::geom::Polygon*>(result2.get())->getNumInteriorRing() );
}

// Test for #1101 - Non-empty negative buffer of 5-pt convex polygon
// Test for https://trac.osgeo.org/geos/ticket/1101 - Non-empty negative buffer of 5-pt convex polygon
template<>
template<>
void object::test<19>
Expand All @@ -556,9 +568,6 @@ template<>
void object::test<20>
()
{
using geos::operation::buffer::BufferOp;
using geos::operation::buffer::BufferParameters;

std::string wkt0("LINESTRING (-20 0, 0 20, 20 0, 0 -20, -20 0)");
GeomPtr g0(wktreader.read(wkt0));

Expand Down Expand Up @@ -609,4 +618,34 @@ void object::test<22>
ensure_equals(result->getArea(), 200);
}

// Checks a bug in the inverted-ring-removal heuristic.
// See https://github.com/libgeos/geos/issues/984
template<>
template<>
void object::test<24>
()
{
std::string wkt("MULTIPOLYGON (((833454.7163917861 6312507.405413097, 833455.3726665961 6312510.208920742, 833456.301153878 6312514.207390314, 833492.2432584754 6312537.770332065, 833493.0901320165 6312536.098774815, 833502.6580673696 6312517.561360772, 833503.9404352929 6312515.0542803425, 833454.7163917861 6312507.405413097)))");

checkBuffer(wkt, -3.8, 0.1,
"POLYGON ((833490.79 6312532.27, 833498.15 6312518, 833459.97 6312512.07, 833490.79 6312532.27))");
checkBuffer(wkt, -7, 0.1,
"POLYGON ((833489.57 6312527.65, 833493.27 6312520.48, 833474.09 6312517.5, 833489.57 6312527.65))");
}

// Checks a bug in the inverted-ring-removal heuristic.
// See https://github.com/libgeos/geos/issues/984
template<>
template<>
void object::test<23>
()
{
std::string wkt("POLYGON ((182719.04521570954238996 224897.14115349075291306, 182807.02887436276068911 224880.64421749324537814, 182808.47314301913138479 224877.25002362736267969, 182718.38701137207681313 224740.00115247094072402, 182711.82697281913715415 224742.08599378637154587, 182717.1393717635946814 224895.61432328051887453, 182719.04521570954238996 224897.14115349075291306))");

checkBuffer(wkt, -5, 0.1,
"POLYGON ((182722 224891.5, 182802 224876.5, 182717 224747, 182722 224891.5))");
checkBuffer(wkt, -30, 0.1,
"POLYGON ((182745.98 224861.57, 182760.51 224858.84, 182745.07 224835.33, 182745.98 224861.57))");
}

} // namespace tut

0 comments on commit 6c1a7ce

Please sign in to comment.