This is an automated email from the ASF dual-hosted git repository. asf-gitbox-commits pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit 4e1d6226dabcc55e07aa5e35b7884179d7d7cb27 Author: jsorel <[email protected]> AuthorDate: Fri Oct 23 08:57:51 2026 +0200 fix(Delaunay): handle infinite loop when triangles become nearly flat --- .../sis/geometries/internal/shared/Planar.java | 71 +++++++++++++--- .../operation/triangulate/delaunay/Delaunay.java | 94 +++++++++++++++------- .../triangulate/delaunay/DelaunayTest.java | 49 ++++++++++- 3 files changed, 174 insertions(+), 40 deletions(-) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Planar.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Planar.java index 7ef7816082..c956bd5bda 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Planar.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Planar.java @@ -16,9 +16,14 @@ */ package org.apache.sis.geometries.internal.shared; +import java.util.ArrayDeque; import java.util.ArrayList; import java.util.Arrays; +import java.util.HashMap; +import java.util.HashSet; import java.util.List; +import java.util.Map; +import java.util.Set; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.Geometries; @@ -126,24 +131,72 @@ public final class Planar { } /** - * Keeps the triangles whose centroid is inside, or outside, the given rings. + * Keeps the triangles inside, or outside, a boundary made of edges of the triangulation. + * The triangle on the left of a boundary edge is inside and the one on its right is outside, + * and the sides are propagated to the neighbors through the other edges. This is exact even for + * very thin triangles along a boundary whose points are slightly off the lines of the rings. + * The triangles not connected to the boundary are classified by their centroid relative to the rings. * * @param xy the points coordinates. - * @param triangles three point indices per triangle. + * @param triangles counterclockwise triangles, three point indices per triangle. + * @param boundary directed boundary edges, two point indices per edge, with the interior on their left. * @param rings the rings, see {@link #inside(List, double, double)}. - * @param inside {@code true} for keeping triangles inside the rings, {@code false} for outside. + * @param inside {@code true} for keeping triangles inside, {@code false} for outside. * @return the kept triangles. */ - public static int[] keep(final double[] xy, final int[] triangles, final List<double[]> rings, final boolean inside) { + public static int[] keep(final double[] xy, final int[] triangles, final int[] boundary, + final List<double[]> rings, final boolean inside) + { + final Map<Long,Integer> byEdge = new HashMap<>(); + for (int i = 0; i < triangles.length; i++) { + byEdge.put(directed(triangles[i], triangles[i - i%3 + (i+1)%3]), i / 3); + } + final Set<Long> walls = new HashSet<>(); + final byte[] side = new byte[triangles.length / 3]; + final ArrayDeque<Integer> queue = new ArrayDeque<>(); + for (int i = 0; i < boundary.length; i += 2) { + final int u = boundary[i], v = boundary[i+1]; + walls.add(edgeKey(u, v)); + label(byEdge.get(directed(u, v)), (byte) 1, side, queue); + label(byEdge.get(directed(v, u)), (byte) -1, side, queue); + } + for (Integer t; (t = queue.poll()) != null;) { + for (int k = 0; k < 3; k++) { + final int a = triangles[t*3+k], b = triangles[t*3+(k+1)%3]; + if (!walls.contains(edgeKey(a, b))) { + label(byEdge.get(directed(b, a)), side[t], side, queue); + } + } + } int n = 0; - for (int i = 0; i < triangles.length; i += 3) { - final int a = triangles[i]*2, b = triangles[i+1]*2, c = triangles[i+2]*2; - if (inside(rings, (xy[a] + xy[b] + xy[c]) / 3, (xy[a+1] + xy[b+1] + xy[c+1]) / 3) == inside) { - System.arraycopy(triangles, i, triangles, n, 3); + final int[] kept = new int[triangles.length]; + for (int t = 0; t < side.length; t++) { + final int a = triangles[t*3]*2, b = triangles[t*3+1]*2, c = triangles[t*3+2]*2; + final boolean in = (side[t] != 0) ? side[t] > 0 + : inside(rings, (xy[a] + xy[b] + xy[c]) / 3, (xy[a+1] + xy[b+1] + xy[c+1]) / 3); + if (in == inside) { + System.arraycopy(triangles, t*3, kept, n, 3); n += 3; } } - return Arrays.copyOf(triangles, n); + return Arrays.copyOf(kept, n); + } + + /** + * Sets the side of a triangle if not already known, and adds it to the triangles to propagate. + */ + private static void label(final Integer triangle, final byte value, final byte[] side, final ArrayDeque<Integer> queue) { + if (triangle != null && side[triangle] == 0) { + side[triangle] = value; + queue.add(triangle); + } + } + + /** + * Returns a key for an edge which depends on its direction. + */ + private static long directed(final int a, final int b) { + return (((long) a) << 32) | (b & 0xFFFFFFFFL); } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/delaunay/Delaunay.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/delaunay/Delaunay.java index 2b2267b0a4..2ceb353160 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/delaunay/Delaunay.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/delaunay/Delaunay.java @@ -19,8 +19,10 @@ package org.apache.sis.geometries.operation.triangulate.delaunay; import java.util.ArrayDeque; import java.util.Arrays; import java.util.HashSet; +import java.util.Random; import java.util.Set; import java.util.concurrent.atomic.AtomicInteger; +import java.util.stream.IntStream; import java.util.logging.Level; import java.util.logging.Logger; import org.opengis.geometry.Envelope; @@ -32,6 +34,7 @@ import org.apache.sis.geometries.mesh.MeshPrimitive.Vertex; import org.apache.sis.geometries.operation.OperationException; import org.apache.sis.geometry.GeneralEnvelope; import org.apache.sis.maths.Array; +import org.apache.sis.maths.Cursor; import org.apache.sis.maths.Maths; import org.apache.sis.maths.NDArrays; import org.apache.sis.maths.Vector; @@ -126,11 +129,19 @@ public final class Delaunay { * Create delaunay triangulation. * This method should be called only once. * @param points not null - * @throws IllegalArgumentException if an index is incorrect + * @throws IllegalArgumentException if an index is incorrect or a position is not finite * @throws OperationException if an algorithm exception occurs */ public void build(MeshPrimitive.Points points) throws IllegalArgumentException, OperationException { ArgumentChecks.ensureNonNull("points", points); + //NaN or infinite values would make the geometric tests fail deep in the algorithm + final Cursor cursor = points.getPositions().cursor(); + while (cursor.next()) { + final Vector<?> position = cursor.samples(); + if (!Double.isFinite(position.get(0)) || !Double.isFinite(position.get(1))) { + throw new IllegalArgumentException("Position " + cursor.coordinate() + " is not finite: " + position); + } + } this.points = points; delaunayTriangulation(); } @@ -315,17 +326,21 @@ public final class Delaunay { // TRIANGULATE ///////////////////////////////////////////////////////// + /* + Insert the points in a random order, with a fixed seed for reproducible results. + Points in a spatially sorted order, for example along a ring or a grid row, would make + the k-d tree degenerate into a chain as deep as the number of points. + */ final Array index = points.getIndex(); - if (index == null) { - //use all points - for (long idx = 0, n = positions.getLength(); idx < n; idx++) { - insertPoint(idx); - } - } else { - //user defined points - for (int idx : index.toArrayInt()) { - insertPoint(idx); - } + final int[] order = (index != null) ? index.toArrayInt() // user defined points + : IntStream.range(0, Math.toIntExact(positions.getLength())).toArray(); + final Random random = new Random(order.length); + for (int i = order.length; --i > 0;) { + final int j = random.nextInt(i + 1); + final int t = order[i]; order[i] = order[j]; order[j] = t; + } + for (int idx : order) { + insertPoint(idx); } } @@ -357,7 +372,8 @@ public final class Delaunay { private void solve(long nearestVertex, Vertex vertex) { OrientedEdge edge = indexToEdge[Math.toIntExact(nearestVertex)]; - for (;;) { + for (int step = walkLimit(); ; step--) { + if (step < 0) throw new OperationException("No triangle found for a vertex, the triangulation is broken."); OrientedTriangle triangle = edge.getTriangle(); if (solve(triangle, vertex)) return; @@ -853,22 +869,17 @@ public final class Delaunay { final OrientedTriangle BDA = B_D.getTriangle(); final OrientedTriangle DBC = D_B.getTriangle(); final Vertex C = DBC.oppositePoint(D_B); - - final double BDAradius = BDA.getDelaunayCircleRadiusSquare(); - final double BDADistance = BDA.distanceToDelaunayCircleCenterSquare(C.getPosition().get(0), C.getPosition().get(1)); - if (BDADistance < BDAradius) { - /* - In the Delaunay algorithm we may encounter the case where two triangles - can be flipped infinitely. This happens when the four points of two adjacent - triangles are contained in the two circle formed by each triangle. - Ensure this is not an infinite loop where triangles are both valid. - */ - final Vertex A = BDA.oppositePoint(B_D); - final Vertex D = B_D.getEnd(); - final OrientedTriangle BCA = new OrientedTriangle(B_D.getStart(), C, A); - if (!BCA.isInDelaunayCircle(D.getPosition().get(0), D.getPosition().get(1))) { - swapTriangles(B_D, true, false); - } + final Vertex A = BDA.oppositePoint(B_D); + /* + Robust test, a circle computed in floating point is not accurate enough for small triangles + far from the origin, for example longitudes and latitudes, and could flip a non convex + quadrilateral, creating an inverted triangle. With an exact test, the quadrilateral is + convex when C is strictly inside the circle, and the flip cannot be undone by the next test. + The convexity is still verified to be safe with flat triangles. + */ + if (Maths.inCircle(BDA.a.getPosition(), BDA.b.getPosition(), BDA.c.getPosition(), C.getPosition()) + && isConvex(A, B_D.getStart(), C, B_D.getEnd())) { + swapTriangles(B_D, true, false); } } @@ -927,6 +938,10 @@ public final class Delaunay { final Vertex B = B_D.getStart(); final Vertex C = DBC.oppositePoint(D_B); final Vertex D = B_D.getEnd(); + if (!isConvex(A, B, C, D)) { + //the flip would create an inverted or flat triangle, breaking the triangulation + throw new OperationException("Cannot flip the diagonal of a non convex quadrilateral."); + } /* swap edge @@ -990,7 +1005,8 @@ public final class Delaunay { OrientedEdge edge = indexToEdge[Math.toIntExact(startIndex)]; - for (;;) { + for (int step = walkLimit(); ; step--) { + if (step < 0) throw new OperationException("No triangle found for a constraint, the triangulation is broken."); if (edge.getEnd().getIndex() == endIndex) { //edge already exist return edge; @@ -1070,7 +1086,25 @@ public final class Delaunay { * @return true if d is in the angle format by vectors ab and ac. */ private static boolean inScanArea(Vertex a, Vertex b, Vertex c, Vertex d) { - return Maths.lineSide(a.getPosition(), b.getPosition(), d.getPosition()) >= 0 + //strict, a point on line AB would create a flat triangle ABD when flipping BC + return Maths.lineSide(a.getPosition(), b.getPosition(), d.getPosition()) > 0 && Maths.lineSide(a.getPosition(), c.getPosition(), d.getPosition()) < 0; } + + /** + * Returns whether the diagonal BD of the counterclockwise quadrilateral ABCD can be replaced + * by the diagonal AC, which requires both new triangles ACD and ABC to be strictly counterclockwise. + */ + private static boolean isConvex(Vertex a, Vertex b, Vertex c, Vertex d) { + return Maths.lineSide(a.getPosition(), c.getPosition(), d.getPosition()) > 0 + && Maths.lineSide(a.getPosition(), b.getPosition(), c.getPosition()) > 0; + } + + /** + * Returns the maximal number of steps of a walk in the triangulation, a walk taking more steps + * is going around in circles because the triangulation is broken. + */ + private int walkLimit() { + return 3 * triangles.size() + 10; + } } diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/delaunay/DelaunayTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/delaunay/DelaunayTest.java index 57a045d658..3a9eae91a3 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/delaunay/DelaunayTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/delaunay/DelaunayTest.java @@ -21,9 +21,11 @@ import org.apache.sis.geometries.mesh.MeshPrimitive; import org.apache.sis.geometries.mesh.MeshPrimitiveComparator; import org.apache.sis.maths.NDArrays; import org.apache.sis.maths.Array; +import org.apache.sis.maths.Maths; import org.apache.sis.referencing.CommonCRS; // Test dependencies +import static org.junit.jupiter.api.Assertions.*; import org.junit.jupiter.api.Test; @@ -64,6 +66,8 @@ public class DelaunayTest { /** * Test triangulation of 4 points. + * The points are on a circle, so both diagonals give a valid Delaunay triangulation. + * Which one is chosen depends on the points insertion order. * 2 3 * +--+ * |\ | @@ -87,7 +91,11 @@ public class DelaunayTest { delaunay.build(points); final Array triangles = delaunay.getTrianglesIndex(); - compare(positions, new int[]{0,1,2, 1,3,2}, triangles); + try { + compare(positions, new int[]{0,1,2, 1,3,2}, triangles); + } catch (IllegalArgumentException e) { + compare(positions, new int[]{0,1,3, 0,3,2}, triangles); + } } /** @@ -333,6 +341,45 @@ public class DelaunayTest { }, triangles); } + /** + * Small triangles far from the origin, as longitudes and latitudes, with collinear points. + * Floating point circle tests are not accurate enough for those coordinates and could flip + * non convex quadrilaterals, creating inverted triangles. All triangles shall be counterclockwise + * and the constraints shall be edges of the triangulation. + */ + @Test + public void testSmallFarTriangles() { + final double x0 = 5.386590116164263, y0 = 43.31551467671114; + final double[] xy = { + x0, y0, // 0 + x0 + 3.082451586e-6, y0 + 3.082451585e-6, // 1, almost collinear with 0 and 4 + 5.386562347412109, 43.31554111681486, // 2 + 5.386562347412109, 43.315521240234375, // 3 + 5.386607056910576, 43.31553161745745, // 4 + x0 + 6.1649e-6, y0 + 6.1649e-6, // 5, on the same line + x0 + 1.2e-5, y0 - 4.1e-6, // 6 + x0 - 7.3e-6, y0 - 9.8e-6, // 7 + x0 + 2.2e-5, y0 + 3.0e-6, // 8 + x0 - 1.1e-5, y0 + 2.9e-5 // 9 + }; + final int[] constraints = {0,1, 1,5, 5,4, 3,0, 7,6, 6,8}; + final int[] triangles = Delaunay.triangulate(xy, constraints); + assertTrue(triangles.length >= 3*8, "Missing triangles."); + final java.util.Set<Long> edges = new java.util.HashSet<>(); + for (int i = 0; i < triangles.length; i += 3) { + final int a = triangles[i]*2, b = triangles[i+1]*2, c = triangles[i+2]*2; + assertTrue(Maths.lineSide(xy[a], xy[a+1], xy[b], xy[b+1], xy[c], xy[c+1]) > 0, "Inverted or flat triangle."); + for (int k = 0; k < 3; k++) { + final int u = triangles[i+k], v = triangles[i+(k+1)%3]; + edges.add(((long) Math.min(u, v) << 32) | Math.max(u, v)); + } + } + for (int i = 0; i < constraints.length; i += 2) { + final int u = constraints[i], v = constraints[i+1]; + assertTrue(edges.contains(((long) Math.min(u, v) << 32) | Math.max(u, v)), "Missing constraint " + u + "-" + v); + } + } + /** * Compare created index, ignoring order of elements. */
