diff --git a/include/geos/geom/LineSegment.h b/include/geos/geom/LineSegment.h index 714c123d9c..ea1e0a7881 100644 --- a/include/geos/geom/LineSegment.h +++ b/include/geos/geom/LineSegment.h @@ -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 { diff --git a/include/geos/operation/buffer/BufferCurveSetBuilder.h b/include/geos/operation/buffer/BufferCurveSetBuilder.h index 2c0ddc0616..e731840bfe 100644 --- a/include/geos/operation/buffer/BufferCurveSetBuilder.h +++ b/include/geos/operation/buffer/BufferCurveSetBuilder.h @@ -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. diff --git a/src/operation/buffer/BufferCurveSetBuilder.cpp b/src/operation/buffer/BufferCurveSetBuilder.cpp index 239b28a8d9..61a5b81329 100644 --- a/src/operation/buffer/BufferCurveSetBuilder.cpp +++ b/src/operation/buffer/BufferCurveSetBuilder.cpp @@ -38,10 +38,12 @@ #include #include #include +#include #include // for min #include #include +#include #include #include #include @@ -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*/ diff --git a/tests/unit/geom/LineSegmentTest.cpp b/tests/unit/geom/LineSegmentTest.cpp index 31bb2dc7c6..71cb7892d4 100644 --- a/tests/unit/geom/LineSegmentTest.cpp +++ b/tests/unit/geom/LineSegmentTest.cpp @@ -11,6 +11,10 @@ // std #include +using geos::geom::Coordinate; +using geos::geom::CoordinateXY; +using geos::geom::LineSegment; + namespace tut { // // Test Group @@ -18,9 +22,6 @@ namespace tut { 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; @@ -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) {} @@ -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 diff --git a/tests/unit/operation/buffer/BufferOpTest.cpp b/tests/unit/operation/buffer/BufferOpTest.cpp index a324bfb2ad..478704424f 100644 --- a/tests/unit/operation/buffer/BufferOpTest.cpp +++ b/tests/unit/operation/buffer/BufferOpTest.cpp @@ -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 geom = wktreader.read(wkt); + std::unique_ptr 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 geom = wktreader.read(wkt); + std::unique_ptr actual = geom->buffer(dist); + std::unique_ptr expected = wktreader.read(wktExpected); + ensure_equals_geometry(expected.get(), actual.get(), tolerance); + } + private: // noncopyable test_bufferop_data(test_bufferop_data const& other) = delete; @@ -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> @@ -534,7 +546,7 @@ void object::test<18> ensure( 0 == dynamic_cast(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> @@ -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)); @@ -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