diff --git a/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java b/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java index 19e25a0b..f656251c 100644 --- a/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java +++ b/planetiler-core/src/main/java/com/onthegomap/planetiler/geo/GeoUtils.java @@ -408,8 +408,11 @@ public class GeoUtils { * {@code ring}, ignoring repeated and collinear points. */ public static boolean isConvex(LinearRing ring) { + double threshold = 1e-3; + double minPointsToCheck = 10; CoordinateSequence seq = ring.getCoordinateSequence(); - if (seq.size() <= 3) { + int size = seq.size(); + if (size <= 3) { return false; } @@ -418,7 +421,7 @@ public class GeoUtils { double c0y = seq.getY(0); double c1x = Double.NaN, c1y = Double.NaN; int i; - for (i = 1; i < seq.size(); i++) { + for (i = 1; i < size; i++) { c1x = seq.getX(i); c1y = seq.getY(i); if (c1x != c0x || c1y != c0y) { @@ -429,29 +432,39 @@ public class GeoUtils { double dx1 = c1x - c0x; double dy1 = c1y - c0y; - int sign = 0; + double negZ = 1e-20, posZ = 1e-20; - for (; i < seq.size(); i++) { - double c2x = seq.getX(i); - double c2y = seq.getY(i); + // need to wrap around to make sure the triangle formed by last and first points does not change sign + for (; i <= size + 1; i++) { + // first and last point should be the same, so skip index 0 + int idx = i < size ? i : (i + 1 - size); + double c2x = seq.getX(idx); + double c2y = seq.getY(idx); double dx2 = c2x - c1x; double dy2 = c2y - c1y; double z = dx1 * dy2 - dy1 * dx2; - // if z == 0 (with small delta to account for rounding errors) then keep skipping - // points to ignore repeated or collinear points - if (Math.abs(z) < 1e-10) { - continue; + double absZ = Math.abs(z); + + // look for sign changes in the triangles formed by sequential points + // but, we want to allow for rounding errors and small concavities relative to the overall shape + // so track the largest positive and negative threshold for triangle area and compare them once we + // have enough points + boolean extendedBounds = false; + if (z < 0 && absZ > negZ) { + negZ = absZ; + extendedBounds = true; + } else if (z > 0 && absZ > posZ) { + posZ = absZ; + extendedBounds = true; } - int s = z >= 0d ? 1 : -1; - if (sign == 0) { - // on the first non-repeated, non-collinear points, store sign of the area for comparison - sign = s; - } else if (sign != s) { - // the sign of this triangle has changed, not convex - return false; + if (i == minPointsToCheck || (i > minPointsToCheck && extendedBounds)) { + double ratio = negZ < posZ ? negZ / posZ : posZ / negZ; + if (ratio > threshold) { + return false; + } } c1x = c2x; @@ -459,7 +472,7 @@ public class GeoUtils { dx1 = dx2; dy1 = dy2; } - return true; + return (negZ < posZ ? negZ / posZ : posZ / negZ) < threshold; } /** diff --git a/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java b/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java index 84e46d72..206710f5 100644 --- a/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java +++ b/planetiler-core/src/test/java/com/onthegomap/planetiler/geo/GeoUtilsTest.java @@ -206,7 +206,7 @@ class GeoUtilsTest { } @Test - void testBarelyConcaveRectangle() { + void testBarelyConcaveTriangle() { assertConvex(false, newLinearRing( 0, 0, 1, 0, @@ -216,6 +216,24 @@ class GeoUtilsTest { )); } + @Test + void testAllowVerySmallConcavity() { + assertConvex(true, newLinearRing( + 0, 0, + 1, 0, + 1, 1, + 0.5001, 0.5, + 0, 0 + )); + assertConvex(true, newLinearRing( + 0, 0, + 1, 0, + 1, 1, + 0.5, 0.4999, + 0, 0 + )); + } + @Test void test5PointsConcave() { assertConvex(false, newLinearRing( @@ -295,7 +313,11 @@ class GeoUtilsTest { LinearRing flipped = flip ? (LinearRing) AffineTransformation.scaleInstance(-1, 1).transform(rotated) : rotated; for (boolean reverse : new boolean[]{false, true}) { LinearRing reversed = reverse ? flipped.reverse() : flipped; - assertEquals(isConvex, isConvex(reversed), "rotation=" + rotation + " flip=" + flip + " reverse=" + reverse); + for (double scale : new double[]{1, 1e-2, 1 / Math.pow(2, 14) / 4096}) { + LinearRing scaled = (LinearRing) AffineTransformation.scaleInstance(scale, scale).transform(reversed); + assertEquals(isConvex, isConvex(scaled), + "rotation=" + rotation + " flip=" + flip + " reverse=" + reverse + " scale=" + scale); + } } } }