Skip to content

Commit

Permalink
Port locationtech/jts#655, Fix buffer to use largest enclosed area fo…
Browse files Browse the repository at this point in the history 
…r invalid rings

Closes #732
  • Loading branch information
@pramsey
pramsey committed Dec 30, 2020
1 parent 5b722cf commit efe09d6
Show file tree
Hide file tree
Showing 7 changed files with 103 additions and 101 deletions.
26 changes: 26 additions & 0 deletions include/geos/algorithm/Orientation.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -89,6 +89,32 @@ class GEOS_DLL Orientation {
*/
static bool isCCW(const geom::CoordinateSequence* ring);

/**
* Tests if a ring defined by a CoordinateSequence is
* oriented counter-clockwise, using the signed area of the ring.
*
* * The list of points is assumed to have the first and last points equal.
* * This handles coordinate lists which contain repeated points.
* * This handles rings which contain collapsed segments
* (in particular, along the top of the ring).
* * This handles rings which are invalid due to self-intersection
*
* This algorithm is guaranteed to work with valid rings.
* For invalid rings (containing self-intersections),
* the algorithm determines the orientation of
* the largest enclosed area (including overlaps).
* This provides a more useful result in some situations, such as buffering.
*
* However, this approach may be less accurate in the case of
* rings with almost zero area.
* (Note that the orientation of rings with zero area is essentially
* undefined, and hence non-deterministic.)
*
* @param ring a CoordinateSequence forming a ring (with first and last point identical)
* @return true if the ring is oriented counter-clockwise.
*/
static bool isCCWArea(const geom::CoordinateSequence* ring);

};


Expand Down
85 changes: 4 additions & 81 deletions src/algorithm/Orientation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@
#include <cmath>
#include <vector>

#include <geos/algorithm/Area.h>
#include <geos/algorithm/Orientation.h>
#include <geos/algorithm/CGAlgorithmsDD.h>
#include <geos/geom/CoordinateSequence.h>
Expand Down Expand Up @@ -135,92 +136,14 @@ Orientation::isCCW(const geom::CoordinateSequence* ring)
}
}

#if 0
/* public static */
bool
Orientation::isCCW(const geom::CoordinateSequence* ring)
Orientation::isCCWArea(const geom::CoordinateSequence* ring)
{
// sanity check
if(ring->getSize() < 4) {
throw util::IllegalArgumentException("Ring has fewer than 4 points, so orientation cannot be determined");
}

// # of points without closing endpoint
const std::size_t nPts = ring->getSize() - 1;
assert(nPts >= 3); // This is here for scan-build

// find highest point
const geom::Coordinate* hiPt = &ring->getAt(0);
std::size_t hiIndex = 0;
for(std::size_t i = 1; i <= nPts; ++i) {
const geom::Coordinate* p = &ring->getAt(i);
if(p->y > hiPt->y) {
hiPt = p;
hiIndex = i;
}
}

// find distinct point before highest point
auto iPrev = hiIndex;
do {
if(iPrev == 0) {
iPrev = nPts;
}
iPrev = iPrev - 1;
}
while(ring->getAt(iPrev) == *hiPt && iPrev != hiIndex);

// find distinct point after highest point
auto iNext = hiIndex;
do {
iNext = (iNext + 1) % nPts;
}
while(ring->getAt(iNext) == *hiPt && iNext != hiIndex);

const geom::Coordinate* prev = &ring->getAt(iPrev);
const geom::Coordinate* next = &ring->getAt(iNext);

/*
* This check catches cases where the ring contains an A-B-A
* configuration of points.
* This can happen if the ring does not contain 3 distinct points
* (including the case where the input array has fewer than 4 elements),
* or it contains coincident line segments.
*/
if(prev->equals2D(*hiPt) || next->equals2D(*hiPt) ||
prev->equals2D(*next)) {
return false;
// MD - don't bother throwing exception,
// since this isn't a complete check for ring validity
//throw IllegalArgumentException("degenerate ring (does not contain 3 distinct points)");
}

int disc = Orientation::index(*prev, *hiPt, *next);

/*
* If disc is exactly 0, lines are collinear.
* There are two possible cases:
* (1) the lines lie along the x axis in opposite directions
* (2) the lines lie on top of one another
*
* (1) is handled by checking if next is left of prev ==> CCW
* (2) should never happen, so we're going to ignore it!
* (Might want to assert this)
*/
bool isCCW = false;
return algorithm::Area::ofRingSigned(ring) < 0;
}

if(disc == 0) {
// poly is CCW if prev x is right of next x
isCCW = (prev->x > next->x);
}
else {
// if area is positive, points are ordered CCW
isCCW = (disc > 0);
}

return isCCW;
}
#endif


} // namespace geos.algorithm
Expand Down
5 changes: 3 additions & 2 deletions src/operation/buffer/OffsetCurveSetBuilder.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -320,8 +320,9 @@ OffsetCurveSetBuilder::addRingSide(const CoordinateSequence* coord,
#if GEOS_DEBUG
std::cerr << "OffsetCurveSetBuilder::addPolygonRing: CCW: " << Orientation::isCCW(coord) << std::endl;
#endif
if(coord->size() >= LinearRing::MINIMUM_VALID_SIZE
&& Orientation::isCCW(coord)) {
if(coord->size() >= LinearRing::MINIMUM_VALID_SIZE &&
Orientation::isCCWArea(coord))
{
leftLoc = cwRightLoc;
rightLoc = cwLeftLoc;
#if GEOS_DEBUG
Expand Down
2 changes: 1 addition & 1 deletion tests/unit/Makefile.am
View file
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ geos_unit_SOURCES = \
algorithm/AngleTest.cpp \
algorithm/AreaTest.cpp \
algorithm/CGAlgorithms/computeOrientationTest.cpp \
algorithm/CGAlgorithms/isCCWTest.cpp \
algorithm/CGAlgorithms/OrientationIsCCWTest.cpp \
algorithm/CGAlgorithms/isPointInRingTest.cpp \
algorithm/CGAlgorithms/signedAreaTest.cpp \
algorithm/ConvexHullTest.cpp \
Expand Down
43 changes: 27 additions & 16 deletions ...unit/algorithm/CGAlgorithms/isCCWTest.cpp → ...thm/CGAlgorithms/OrientationIsCCWTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ struct test_isccw_data {
}

void
checkOrientationCCW(bool expectedCCW, const std::string& wkt)
checkCCW(bool expectedCCW, const std::string& wkt)
{
GeometryPtr geom(reader_.read(wkt));
geos::geom::Polygon* poly = dynamic_cast<geos::geom::Polygon*>(geom.get());
Expand All @@ -50,6 +50,17 @@ struct test_isccw_data {
ensure_equals("CoordinateSequence isCCW", expectedCCW, actualCCW);
}

void
checkCCWArea(bool expectedCCWArea, const std::string& wkt)
{
GeometryPtr geom(reader_.read(wkt));
geos::geom::Polygon* poly = dynamic_cast<geos::geom::Polygon*>(geom.get());
ensure("WKT must be POLYGON)", poly != nullptr);
const geos::geom::CoordinateSequence* cs = poly->getExteriorRing()->getCoordinatesRO();
bool actualCCWArea = Orientation::isCCWArea(cs);
ensure_equals("CoordinateSequence isCCWArea", expectedCCWArea, actualCCWArea);
}

void
checkHexOrientationCCW(bool expectedCCW, std::istringstream& wkt)
{
Expand All @@ -64,7 +75,7 @@ struct test_isccw_data {
typedef test_group<test_isccw_data> group;
typedef group::object object;

group test_isccw_group("geos::algorithm::CGAlgorithms::isCCW");
group test_isccw_group("geos::algorithm::CGAlgorithms::OrientationIsCCW");

//
// Test Cases
Expand All @@ -77,7 +88,7 @@ void object::test<1>
()
{
const std::string wkt("POLYGON ((60 180, 140 240, 140 240, 140 240, 200 180, 120 120, 60 180))");
checkOrientationCCW(false, wkt);
checkCCW(false, wkt);
}

// 2 - Test if coordinates of polygon are counter-clockwise oriented
Expand All @@ -87,7 +98,7 @@ void object::test<2>
()
{
const std::string wkt("POLYGON ((60 180, 140 120, 100 180, 140 240, 60 180))");
checkOrientationCCW(true, wkt);
checkCCW(true, wkt);
}

// 3 - Test the same polygon as in test No 2 but with duplicated top point
Expand All @@ -97,7 +108,7 @@ void object::test<3>
()
{
const std::string wkt("POLYGON ((60 180, 140 120, 100 180, 140 240, 140 240, 60 180))");
checkOrientationCCW(true, wkt);
checkCCW(true, wkt);
}

// 4 - Test orientation the narrow (almost collapsed) ring
Expand Down Expand Up @@ -133,7 +144,7 @@ void object::test<6>
()
{
const std::string wkt("POLYGON ((1 1, 9 1, 5 9, 1 1))");
checkOrientationCCW(true, wkt);
checkCCW(true, wkt);
}

// testFlatTopSegment
Expand All @@ -143,7 +154,7 @@ void object::test<7>
()
{
const std::string wkt("POLYGON ((100 200, 200 200, 200 100, 100 100, 100 200))");
checkOrientationCCW(false, wkt);
checkCCW(false, wkt);
}

// testFlatMultipleTopSegment
Expand All @@ -153,7 +164,7 @@ void object::test<8>
()
{
const std::string wkt("POLYGON ((100 200, 127 200, 151 200, 173 200, 200 200, 100 100, 100 200))");
checkOrientationCCW(false, wkt);
checkCCW(false, wkt);
}

// testDegenerateRingHorizontal
Expand All @@ -163,7 +174,7 @@ void object::test<9>
()
{
const std::string wkt("POLYGON ((100 200, 100 200, 200 200, 100 200))");
checkOrientationCCW(false, wkt);
checkCCW(false, wkt);
}

// testDegenerateRingAngled
Expand All @@ -173,7 +184,7 @@ void object::test<10>
()
{
const std::string wkt("POLYGON ((100 100, 100 100, 200 200, 100 100))");
checkOrientationCCW(false, wkt);
checkCCW(false, wkt);
}

// testDegenerateRingVertical
Expand All @@ -183,7 +194,7 @@ void object::test<11>
()
{
const std::string wkt("POLYGON ((200 100, 200 100, 200 200, 200 100))");
checkOrientationCCW(false, wkt);
checkCCW(false, wkt);
}

/**
Expand All @@ -196,7 +207,7 @@ void object::test<12>
()
{
const std::string wkt("POLYGON ((10 20, 61 20, 20 30, 50 60, 10 20))");
checkOrientationCCW(false, wkt);
checkCCW(false, wkt);
}

// testABATopFlatSegmentCollapse
Expand All @@ -206,7 +217,7 @@ void object::test<13>
()
{
const std::string wkt("POLYGON ((71 0, 40 40, 70 40, 40 40, 20 0, 71 0))");
checkOrientationCCW(true, wkt);
checkCCW(true, wkt);
}

// testABATopFlatSegmentCollapseMiddleStart
Expand All @@ -216,7 +227,7 @@ void object::test<14>
()
{
const std::string wkt("POLYGON ((90 90, 50 90, 10 10, 90 10, 50 90, 90 90))");
checkOrientationCCW(true, wkt);
checkCCW(true, wkt);
}

// testMultipleTopFlatSegmentCollapseSinglePoint
Expand All @@ -226,7 +237,7 @@ void object::test<15>
()
{
const std::string wkt("POLYGON ((100 100, 200 100, 150 200, 170 200, 200 200, 100 200, 150 200, 100 100))");
checkOrientationCCW(true, wkt);
checkCCW(true, wkt);
}

// testMultipleTopFlatSegmentCollapseFlatTop
Expand All @@ -236,7 +247,7 @@ void object::test<16>
()
{
const std::string wkt("POLYGON ((10 10, 90 10, 70 70, 90 70, 10 70, 30 70, 50 70, 10 10))");
checkOrientationCCW(true, wkt);
checkCCW(true, wkt);
}


Expand Down
18 changes: 18 additions & 0 deletions tests/unit/operation/buffer/BufferOpTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -3,6 +3,7 @@

// tut
#include <tut/tut.hpp>
#include <utility.h>
// geos
#include <geos/operation/buffer/BufferOp.h>
#include <geos/operation/buffer/BufferParameters.h>
Expand Down Expand Up @@ -447,5 +448,22 @@ void object::test<14>
ensure_equals(gBuffer->getGeometryTypeId(), geos::geom::GEOS_MULTIPOLYGON);
}


// This now works since buffer ring orientation is changed to use signed-area test.
// testBowtiePolygonLargestAreaRetained
template<>
template<>
void object::test<15>
()
{
std::string wkt0("POLYGON ((10 10, 50 10, 25 35, 35 35, 10 10))");
GeomPtr g0(wktreader.read(wkt0));
GeomPtr gresult = g0->buffer(0.0);
std::string wkt1("POLYGON ((10 10, 30 30, 50 10, 10 10))");
GeomPtr gexpected(wktreader.read(wkt1));
ensure_equals_geometry(gresult.get(), gexpected.get());
}


} // namespace tut

25 changes: 24 additions & 1 deletion tests/unit/simplify/DouglasPeuckerSimplifierTest.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,6 +2,7 @@
// Test Suite for geos::simplify::DouglasPeuckerSimplifierTest

#include <tut/tut.hpp>
#include <utility.h>
// geos
#include <geos/io/WKTReader.h>
#include <geos/io/WKTWriter.h>
Expand Down Expand Up @@ -353,7 +354,7 @@ void object::test<11>
// 13 - Polygon with inner ring whose extent is less than the simplify distance (#741)
template<>
template<>
void object::test<13>
void object::test<12>
()
{
std::string wkt_in("POLYGON ((0 0,0 1,1 1,0 0),(0.1 0.1,0.2 0.1,0.2 0.2,0.1 0.1))");
Expand All @@ -372,5 +373,27 @@ void object::test<13>
ensure(simplified->equalsExact(expected.get()));
}

/**
* Test that a polygon made invalid by simplification
* is fixed in a sensible way.
* Fixed by buffer(0) area-base orientation
* See https://github.com/locationtech/jts/issues/498
*/
template<>
template<>
void object::test<13>
()
{
std::string wkt_in("POLYGON ((21.32686 47.78723, 21.32386 47.79023, 21.32186 47.80223, 21.31486 47.81023, 21.32786 47.81123, 21.33986 47.80223, 21.33886 47.81123, 21.32686 47.82023, 21.32586 47.82723, 21.32786 47.82323, 21.33886 47.82623, 21.34186 47.82123, 21.36386 47.82223, 21.40686 47.81723, 21.32686 47.78723))");
std::string wkt_ex("POLYGON ((21.32686 47.78723, 21.31486 47.81023, 21.32786 47.81123, 21.33986 47.80223, 21.328068201892744 47.823286782334385, 21.33886 47.82623, 21.34186 47.82123, 21.40686 47.81723, 21.32686 47.78723))");
GeomPtr g(wktreader.read(wkt_in));
GeomPtr expected(wktreader.read(wkt_ex));
GeomPtr simplified = DouglasPeuckerSimplifier::simplify(g.get(), 0.0036);
ensure(simplified->isValid());
ensure_equals_geometry(simplified.get(), expected.get());
}



} // namespace tut

0 comments on commit efe09d6

Please sign in to comment.