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 7f9bdfd15da62d06d5a086dd2fd44225dcccd336 Author: jsorel <[email protected]> AuthorDate: Fri Oct 16 15:57:43 2026 +0200 feat(Geometry): simplify EarClipping --- .../internal/shared/CurveLinearizer.java | 12 +- .../geometries/internal/shared/DefaultPrism.java | 31 +- .../operation/triangulate/EarClipping.java | 320 ++++++------------ .../operation/triangulate/HoleToExterior.java | 357 +++++++-------------- .../operation/triangulate/HoleToExteriorTest.java | 29 ++ 5 files changed, 266 insertions(+), 483 deletions(-) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/CurveLinearizer.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/CurveLinearizer.java index 3ef811015a..5947bd65fa 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/CurveLinearizer.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/CurveLinearizer.java @@ -26,7 +26,6 @@ import org.opengis.geometry.Envelope; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.Curve; import org.apache.sis.geometries.DataPoints; -import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.Geometries; import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.Point; @@ -482,13 +481,8 @@ public final class CurveLinearizer { if (result.size() == n) { return points; } - final DataPointsType type = points.getType(); - final DataPoints copy = GeometryFactory.DEFAULT.createDataPoints(type, result.size()); - for (int i = 0; i < result.size(); i++) { - for (final String name : type.getAttributeNames()) { - copy.setAttribute(i, name, result.get(i).getAttribute(name)); - } - } - return copy; + //copy the points, so that the polyline does not depend on the given ones. + final DataPoints view = ConcatenatedDataPoints.of(result.toArray(Point[]::new)); + return view.subset(0, view.size()); } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultPrism.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultPrism.java index fdda173d48..faf8a6315f 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultPrism.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultPrism.java @@ -45,7 +45,6 @@ import org.apache.sis.maths.Array; import org.apache.sis.maths.NDArrays; import org.apache.sis.maths.SampleSystem; import org.apache.sis.maths.Vector; -import org.apache.sis.maths.Vectors; import org.apache.sis.measure.NumberRange; import org.apache.sis.referencing.CRS; import org.apache.sis.util.Utilities; @@ -403,30 +402,24 @@ public non-sealed class DefaultPrism extends AbstractGeometry implements Prism { * @param values extrusion value of each copied position. */ private DataPoints copy(final DataPoints source, final int[] indices, final double[] values) { - final DataPointsType type = source.getType(); - final SampleSystem system = SampleSystem.of(crs); + final DataPoints copy = source.subset(indices); final Map<String,Array> attributes = new LinkedHashMap<>(); - for (final String name : type.getAttributeNames()) { - final boolean isPosition = DataPointsType.ATT_POSITION.equals(name); - attributes.put(name, NDArrays.of(isPosition ? system : type.getAttributeSystem(name), - type.getAttributeType(name), indices.length)); + for (final String name : copy.getType().getAttributeNames()) { + attributes.put(name, copy.getAttributeArray(name)); } - final Array positions = attributes.get(DataPointsType.ATT_POSITION); - final Vector<?> position = Vectors.create(system, type.getAttributeType(DataPointsType.ATT_POSITION)); - final int dimension = source.getDimension(); + //the positions are in the reference system of the extrusion, which may have one more dimension. + final Array positions = NDArrays.of(SampleSystem.of(crs), + copy.getType().getAttributeType(DataPointsType.ATT_POSITION), indices.length); for (int i = 0; i < indices.length; i++) { - final Vector<?> p = source.getPosition(indices[i]); - for (int d = 0; d < dimension; d++) { - position.set(d, p.get(d)); + Vector<?> position = copy.getPosition(i); + if (axis < position.getDimension()) { + position.set(axis, values[i]); + } else { + position = position.extend(values[i]); } - position.set(axis, values[i]); positions.set(i, position); - for (final Map.Entry<String,Array> entry : attributes.entrySet()) { - if (entry.getValue() != positions) { - entry.getValue().set(i, source.getAttribute(indices[i], entry.getKey())); - } - } } + attributes.put(DataPointsType.ATT_POSITION, positions); return factory.createDataPoints(attributes); } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java index 5f85efcbbd..41c77e90a8 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java @@ -18,12 +18,12 @@ package org.apache.sis.geometries.operation.triangulate; import java.util.ArrayList; import java.util.Arrays; -import java.util.HashMap; +import java.util.LinkedHashMap; import java.util.List; import java.util.Map; +import java.util.stream.IntStream; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; -import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.mesh.MeshPrimitive; import org.apache.sis.geometries.surface.Polygon; import org.apache.sis.maths.Array; @@ -41,7 +41,7 @@ import org.apache.sis.util.ArraysExt; * * <p>The holes of the polygon are first joined to its exterior ring by {@link HoleToExterior}, * then the ring is triangulated. The triangulation works on the indices of the positions of - * the ring: the triangles reference the vertices instead of copying their positions, so that + * the polygon: the triangles reference the vertices instead of copying their positions, so that * all the attributes of the vertices are kept.</p> * * @author Johann Sorel @@ -57,19 +57,6 @@ public class EarClipping { public record Triangulation(DataPoints vertices, int[] triangles) { } - /** - * Values of all the attributes of a vertex, for finding the vertices which are identical. - */ - private record Key(double[] values) { - @Override public boolean equals(final Object other) { - return (other instanceof Key k) && Arrays.equals(values, k.values); - } - - @Override public int hashCode() { - return Arrays.hashCode(values); - } - } - /** * Type of the angle made by the ring at a vertex, computed when first needed. */ @@ -80,48 +67,38 @@ public class EarClipping { } /** - * Positions of the ring to triangulate, the last one at the same position as the first one. + * Positions of the polygon to triangulate, those of the exterior ring followed by those of the holes. */ private DataPoints source; /** - * Positions of the ring, with their coordinates in the order used by the triangulation. + * Positions of the polygon, with their coordinates in the order used by the triangulation. * The triangulation works on the two first coordinates, in the plane where the polygon has its largest area. */ private Vector<?>[] working; /** - * Indices of the three vertices of each triangle found so far. + * Indices in {@link #source} of the vertices of the ring joining the holes to the exterior ring, without + * repeated consecutive positions and without the closing position. The ring is clockwise in the plane + * of the triangulation. */ - private final List<int[]> triangles = new ArrayList<>(); + private int[] ring; - /** number of coordinates, may be inferior to coords length. */ - private int nbCoords; - /** indices of the vertices to triangulate. */ - private int[] coords; /** - * The type of angle made by the ring at each vertex to triangulate, - * stored to avoid computing it each time. + * For each vertex of the {@linkplain #ring}, the previous and next vertices which are not cut yet. + * Cutting an ear unlinks its vertex from the ring. */ - private Angle[] coordType; + private int[] previous, next; - //current segment analyzed - private int indexPrevious; - private int index; - private int indexNext; - private Vector<?> t1; //previous point - private Vector<?> t2; //point at index - private Vector<?> t3; //next point + /** + * The type of angle made by the ring at each vertex, stored to avoid computing it each time. + */ + private Angle[] angles; - private void reset(){ - //reset values - source = null; - working = null; - triangles.clear(); - nbCoords = 0; - coords = null; - coordType = null; - } + /** + * Indices in {@link #source} of the three vertices of each triangle found so far. + */ + private final List<int[]> triangles = new ArrayList<>(); /** * Triangulates the given polygon. @@ -130,14 +107,12 @@ public class EarClipping { * @return the vertices of the triangles and their indices. */ public Triangulation triangulateIndexed(Polygon geometry) { - reset(); - source = HoleToExterior.apply(geometry).getExteriorRing().getDataPoints(); - final int[] axes = workingAxes(source.getAttributeArray(DataPointsType.ATT_POSITION)); + source = geometry.getDataPoints(); + final int[] axes = workingAxes(geometry.getExteriorRing().getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION)); working = new Vector<?>[source.size()]; - for (int i = 0; i < working.length; i++) { - working[i] = source.getPosition(i).swizzle(axes); - } - run(); + Arrays.setAll(working, (i) -> source.getPosition(i).swizzle(axes)); + triangles.clear(); + run(HoleToExterior.ring(geometry)); return result(); } @@ -158,207 +133,116 @@ public class EarClipping { /** * Returns the order of the coordinates to use for triangulating a polygon with the given positions, * to be given to {@link Vector#swizzle(int...)}. The triangulation works on the two first coordinates - * of that order. The axis along which the normal of the polygon is the largest, computed with Newell's - * method, is moved after the two axes of the projection plane, so that a vertical polygon is not - * degenerate. If the polygon can be triangulated in its (x,y) plane, the order is unchanged. + * of that order. The axis along which the {@linkplain Maths#calculatePolygonNormal(Array) normal of the + * polygon} is the largest is moved after the two axes of the projection plane, so that a vertical polygon + * is not degenerate. If the polygon can be triangulated in its (x,y) plane, the order is unchanged. */ static int[] workingAxes(final Array positions) { final int dim = positions.getDimension(); - final int[] axes = new int[dim]; - for (int i = 0; i < dim; i++) { - axes[i] = i; - } - if (dim < 3) { - return axes; - } - final double[] c = positions.toArrayDouble(); - final int n = c.length / dim; - double nx = 0, ny = 0, nz = 0; - for (int i = 0; i < n; i++) { - final int a = i * dim, b = ((i + 1) % n) * dim; - nx += (c[a+1] - c[b+1]) * (c[a+2] + c[b+2]); - ny += (c[a+2] - c[b+2]) * (c[a] + c[b]); - nz += (c[a] - c[b]) * (c[a+1] + c[b+1]); + final int[] axes = ArraysExt.range(0, dim); + if (dim >= 3) { + final double[] normal = Maths.calculatePolygonNormal(positions); + final double nx = Math.abs(normal[0]), ny = Math.abs(normal[1]), nz = Math.abs(normal[2]); + final int drop = (nz >= nx && nz >= ny) ? 2 : (nx >= ny) ? 0 : 1; + //shift the following axes of the plane to the place of the dropped axis, which goes after them. + System.arraycopy(axes, drop + 1, axes, drop, 2 - drop); + axes[2] = drop; } - nx = Math.abs(nx); ny = Math.abs(ny); nz = Math.abs(nz); - if (nz >= nx && nz >= ny) { - return axes; - } - final int drop = (nx >= ny) ? 0 : 1; - int k = 0; - for (int i = 0; i < 3; i++) { - if (i != drop) axes[k++] = i; - } - axes[k] = drop; return axes; } /** * Builds the result from the triangles found. Only the vertices used by a triangle are kept, in the - * order of the ring, and the vertices having the same values for all attributes are merged. Those are - * the last position of the ring, which repeats the first one, and the positions repeated by the cuts - * joining the holes. + * order of their first use. The positions repeated by the cuts joining the holes have the same index, + * so they are a single vertex. */ private Triangulation result() { - final int count = working.length; - final boolean[] used = new boolean[count]; + final Map<Integer,Integer> vertices = new LinkedHashMap<>(); // index in the source → index in the result. + final int[] indices = new int[triangles.size() * 3]; + int k = 0; for (final int[] triangle : triangles) { for (final int vertex : triangle) { - used[vertex] = true; - } - } - final DataPointsType type = source.getType(); - final List<String> names = type.getAttributeNames(); - final Map<Key,Integer> distinct = new HashMap<>(); - final List<Integer> kept = new ArrayList<>(); - final int[] remap = new int[count]; - for (int i = 0; i < count; i++) { - if (used[i]) { - double[] values = new double[0]; - for (final String name : names) { - final double[] value = source.getAttribute(i, name).toArrayDouble(); - values = Arrays.copyOf(values, values.length + value.length); - System.arraycopy(value, 0, values, values.length - value.length, value.length); - } - final Integer previous = distinct.putIfAbsent(new Key(values), kept.size()); - if (previous != null) { - remap[i] = previous; - } else { - remap[i] = kept.size(); - kept.add(i); - } - } - } - final DataPoints vertices = GeometryFactory.DEFAULT.createDataPoints(type, kept.size()); - for (int i = 0; i < kept.size(); i++) { - for (final String name : names) { - vertices.setAttribute(i, name, source.getAttribute(kept.get(i), name)); + indices[k++] = vertices.computeIfAbsent(vertex, (key) -> vertices.size()); } } - final int[] indices = new int[triangles.size() * 3]; - for (int i = 0, k = 0; i < triangles.size(); i++) { - for (final int vertex : triangles.get(i)) { - indices[k++] = remap[vertex]; - } - } - return new Triangulation(vertices, indices); + return new Triangulation(source.subset(vertices.keySet().stream().mapToInt(Integer::intValue).toArray()), indices); } /** - * Triangulates the ring. + * Triangulates the given ring by cutting its ears one after the other. An ear is a convex vertex whose + * triangle with its two neighbors contains no other vertex. Cutting it removes the vertex from the + * ring, which stays a valid polygon, until only one triangle is left. + * + * @param joined indices in {@link #source} of the ring positions, the last one at the location of the first one. */ - private void run() { - if (working.length == 0) { - return; - } - final List<Integer> border = new ArrayList<>(working.length); - for (int i = 0; i < working.length; i++) { - border.add(i); - } - - //remove any neighor points overlaping - Vector<?> t = working[border.get(0)]; - for(int i=1,n=border.size();i<n;i++){ - final Vector<?> candidate = working[border.get(i)]; - if(candidate.equals(t)){ - border.remove(i); - i--; - n--; - }else{ - t = candidate; + private void run(final int[] joined) { + //remove the repeated consecutive positions, and the closing position which repeats the first one. + ring = IntStream.range(0, joined.length - 1).filter((i) -> !working[joined[i]].equals(working[joined[i + 1]])) + .map((i) -> joined[i]).toArray(); + if (!Maths.isClockWise(Arrays.stream(joined).mapToObj((i) -> working[i]).toList())) { + ArraysExt.reverse(ring); + } + final int n = ring.length; + previous = new int[n]; + next = new int[n]; + for (int i = 0; i < n; i++) { + previous[i] = (i + n - 1) % n; + next[i] = (i + 1) % n; + } + angles = new Angle[n]; + Arrays.fill(angles, Angle.UNKNOWN); + + //walk around the ring, until a whole turn finds no ear, which should not happen if the geometry is correct. + int vertex = 0; + for (int remaining = n, attempts = 0; remaining >= 3 && attempts < remaining; vertex = next[vertex]) { + if (isConvex(vertex) && isEar(vertex)) { + triangles.add(new int[] {ring[previous[vertex]], ring[vertex], ring[next[vertex]]}); + //unlink the vertex, the angles at its neighbors may change. + next[previous[vertex]] = next[vertex]; + previous[next[vertex]] = previous[vertex]; + angles[previous[vertex]] = Angle.UNKNOWN; + angles[next[vertex]] = Angle.UNKNOWN; + remaining--; + attempts = 0; + } else { + attempts++; } } + } - nbCoords = border.size(); - coordType = new Angle[nbCoords]; - Arrays.fill(coordType, Angle.UNKNOWN); - coords = border.stream().mapToInt(Integer::intValue).toArray(); - - //flip coordinates if not clockwise - if(!Maths.isClockWise(Arrays.stream(coords).mapToObj((i) -> working[i]).toList())){ - ArraysExt.reverse(coords); - } - - //now cut ears progressively - int nbcut = 0; - while(nbCoords > 3){ - nbcut = 0; - for(index=0;index<nbCoords-2 && nbCoords>3;){ - indexPrevious = (index==0) ? nbCoords-2 : index-1 ; - indexNext = index+1; - - t1 = working[coords[indexPrevious]]; - t2 = working[coords[index]]; - t3 = working[coords[indexNext]]; - - if(isConvex(index) && isEar()){ - //reset angle type, we remove a point so it might change - coordType[indexPrevious]=Angle.UNKNOWN; - coordType[indexNext]=Angle.UNKNOWN; - - triangles.add(new int[]{coords[indexPrevious], coords[index], coords[indexNext]}); - Maths.removeWithin(coords, index); - Maths.removeWithin(coordType, index); - nbCoords--; - nbcut++; - //we do not increment since we have remove the point - }else{ - index++; - } - } - - if(nbcut == 0){ - //this should not happen if the geometry is correct - //System.out.println("Triangulation failed. no ear to cut."); - return; - } - } + /** + * Returns the working position of the given vertex of the ring. + */ + private Vector<?> position(final int vertex) { + return working[ring[vertex]]; } /** - * @return true if the ring is convex at the vertex of given index + * @return true if the ring is convex at the given vertex */ - private boolean isConvex(int idx){ - if(coordType[idx]==Angle.UNKNOWN){ - //calculate angle type - final Vector<?> s1 = working[coords[(idx==0) ? (nbCoords-2) : (idx-1)]]; - final Vector<?> s2 = working[coords[idx]]; - final Vector<?> s3 = working[coords[idx+1]]; - final double side = Maths.lineSide(s1, s3, s2); - coordType[idx] = side>0 ? Angle.CONVEX : Angle.CONCAVE; + private boolean isConvex(final int vertex) { + if (angles[vertex] == Angle.UNKNOWN) { + final boolean convex = Maths.isConvex(position(previous[vertex]), position(vertex), position(next[vertex]), true); + angles[vertex] = convex ? Angle.CONVEX : Angle.CONCAVE; } - return coordType[idx] == Angle.CONVEX; + return angles[vertex] == Angle.CONVEX; } /** - * Check this segment triangle is an ear. - * Does not contain any other point. - * @return true if no other vertex is inside the triangle + * Returns whether the triangle made by the given vertex and its two neighbors contains no other vertex + * of the ring. Only the concave vertices need to be tested. A cut joining a hole repeats its end points, + * so a copy located on a corner of the triangle does not prevent the ear. */ - private boolean isEar() { - final double[] a = new double[]{t1.get(0), t1.get(1)}; - final double[] b = new double[]{t2.get(0), t2.get(1)}; - final double[] c = new double[]{t3.get(0), t3.get(1)}; - - for(int i=0; i<nbCoords-1; i++){ - //test only concave points - if (!isConvex(i)){ - //check it's not one of the current segment - if(i!=index && i!=indexNext && i!=indexPrevious){ - final Vector<?> v = working[coords[i]]; - final double[] p = new double[]{v.get(0),v.get(1)}; - //a hole bridge duplicates its end points, a copy located on a triangle vertex does not block the ear - if (Arrays.equals(p, a) || Arrays.equals(p, b) || Arrays.equals(p, c)) { - continue; - } - //check if it's in the segment triangle - if(Maths.inTriangle(a, b, c, p)){ - return false; - } - } + private boolean isEar(final int vertex) { + final Vector<?> t1 = position(previous[vertex]); + final Vector<?> t2 = position(vertex); + final Vector<?> t3 = position(next[vertex]); + for (int i = next[next[vertex]]; i != previous[vertex]; i = next[i]) { + final Vector<?> p = position(i); + if (!isConvex(i) && !p.equals(t1) && !p.equals(t2) && !p.equals(t3) && Maths.inTriangle(t1, t2, t3, p)) { + return false; } } - return true; } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/HoleToExterior.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/HoleToExterior.java index 4ee173c257..91b2b915bd 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/HoleToExterior.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/HoleToExterior.java @@ -21,6 +21,7 @@ import java.util.Arrays; import java.util.Collections; import java.util.Comparator; import java.util.List; +import java.util.stream.IntStream; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.GeometryFactory; @@ -38,11 +39,9 @@ import org.apache.sis.maths.Vector; * rightmost one to the leftmost one. Each hole is joined by a cut going from a vertex of the ring to the * rightmost position of the hole, around the hole, and back to the same vertex of the ring. The vertex * of the ring is one which can be reached without crossing any edge, found as described by David Eberly - * in <cite>Triangulation by Ear Clipping</cite>. The two sides of a cut are made of the same positions, - * so the resulting ring touches itself along the cut, and no position is added.</p> - * - * <p>If no such vertex is found, which may happen with an invalid polygon, the cut starts at a new - * position on the nearest segment of the ring, whose attributes are interpolated along that segment.</p> + * in <cite>Triangulation by Ear Clipping</cite>. If there is none, which may happen with an invalid polygon, + * the cut starts from the nearest end of the segment of the ring nearest to the hole. The two sides of a cut + * are made of the same positions, so the resulting ring touches itself along the cut, and no position is added.</p> * * <p>The positions are compared in the coordinate plane where the polygon has its largest area, * so that a vertical polygon is not degenerate. This is the first stage of {@link EarClipping}.</p> @@ -52,144 +51,123 @@ import org.apache.sis.maths.Vector; public final class HoleToExterior { /** - * Positions of the polygon, of the exterior ring followed by those of the holes. - * The vertex of index <var>i</var> is the position <var>i</var> of this {@code DataPoints} - * if <var>i</var> is smaller than its size, or a vertex added for a cut otherwise. + * The polygon whose holes are joined to the exterior ring. */ - private final DataPoints source; + private final Polygon polygon; /** - * Positions of all the vertices, with their coordinates in the order used for comparing them. + * Positions of the polygon, the exterior ring first, with their coordinates in the order used for comparing them. */ - private final List<Vector<?>> working = new ArrayList<>(); + private final Vector<?>[] working; /** - * For each vertex added to start a cut, the indices of the positions of the polygon it is interpolated from. - * The vertex at index <var>i</var> in this list has the index {@code source.size() + i}. + * Whether the exterior ring is clockwise. The joined ring has the same orientation. */ - private final List<int[]> addedIndices = new ArrayList<>(); + private final boolean clockwise; /** - * For each vertex added to start a cut, the weights of the positions of the polygon it is interpolated from. + * Prepares the joining of the holes of the given polygon. */ - private final List<double[]> addedWeights = new ArrayList<>(); + private HoleToExterior(final Polygon polygon) { + this.polygon = polygon; + final DataPoints exterior = polygon.getExteriorRing().getDataPoints(); + final DataPoints positions = polygon.getDataPoints(); + final int[] axes = EarClipping.workingAxes(exterior.getAttributeArray(DataPointsType.ATT_POSITION)); + working = new Vector<?>[positions.size()]; + Arrays.setAll(working, (i) -> positions.getPosition(i).swizzle(axes)); + clockwise = Maths.isClockWise(Arrays.asList(working).subList(0, exterior.size())); + } /** - * Prepares the joining of the holes of the given polygon. + * Returns the exterior ring of the given polygon with its holes joined to it, as the indices of its positions + * in {@code polygon.getDataPoints()}, where the positions of the exterior ring are followed by those of the holes. + * The positions at both ends of the cuts are repeated, but no position is added. The ring is closed: its last + * index is the one of the last position of the exterior ring, which is at the same location as the first one. + * + * @param polygon the polygon whose holes are to be joined to the exterior ring. + * @return indices of the positions of the joined ring, or an empty array if the exterior ring is empty. */ - private HoleToExterior(final Polygon polygon) { - source = polygon.getDataPoints(); - final int[] axes = EarClipping.workingAxes( - polygon.getExteriorRing().getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION)); - for (int i = 0, n = source.size(); i < n; i++) { - working.add(source.getPosition(i).swizzle(axes)); + public static int[] ring(final Polygon polygon) { + if (polygon.getExteriorRing().getDataPoints().isEmpty()) { + return new int[0]; } + return new HoleToExterior(polygon).join(); } /** - * Returns a polygon without hole covering the same area as the given polygon. - * The exterior ring of the returned polygon is the exterior ring of the given polygon with the holes - * joined to it. Its positions are those of the given polygon, with all their attributes, some of them - * repeated at both ends of the cuts. Positions are added only for the cuts which could not start - * from a vertex. + * Returns a polygon without hole covering the same area as the given polygon. The exterior ring of the + * returned polygon is the {@linkplain #ring(Polygon) joined ring}, made of the positions of the given + * polygon with all their attributes. * * @param polygon the polygon whose holes are to be joined to the exterior ring. * @return the given polygon if it has no hole, or a new polygon without hole otherwise. */ public static Polygon apply(final Polygon polygon) { - if (polygon.getExteriorRing().getDataPoints().isEmpty()) { + if (polygon.getExteriorRing().getDataPoints().isEmpty() + || polygon.getInteriorRings().stream().allMatch((hole) -> hole.getDataPoints().isEmpty())) { return polygon; } - for (final LinearRing hole : polygon.getInteriorRings()) { - if (!hole.getDataPoints().isEmpty()) { - return new HoleToExterior(polygon).join(polygon); - } - } - return polygon; + final GeometryFactory factory = GeometryFactory.DEFAULT; + return factory.createPolygon(factory.createLinearRing(polygon.getDataPoints().subset(ring(polygon))), null); } /** - * Joins the holes of the given polygon, which is the polygon given to the constructor. + * Joins the holes of the polygon to its exterior ring. + * + * @return indices of the positions of the joined ring. */ - private Polygon join(final Polygon polygon) { + private int[] join() { //the positions of the polygon are those of the exterior ring followed by those of the holes. int start = polygon.getExteriorRing().getDataPoints().size(); final List<Integer> border = range(0, start); final List<List<Integer>> holes = new ArrayList<>(); for (final LinearRing hole : polygon.getInteriorRings()) { final int size = hole.getDataPoints().size(); - if (size != 0) { + if (size > 1) { holes.add(range(start, size)); } start += size; } - //the holes shall turn in the opposite direction of the exterior ring, so that the joined ring has a single orientation. - final double orientation = Math.signum(area(border)); + //the holes shall turn in the opposite direction of the exterior ring, so that the joined ring has a single + //orientation. Each hole is reordered to start and end at its vertex of maximum x, where its cut starts. for (final List<Integer> loop : holes) { - if (Math.signum(area(loop)) == orientation) { + if (Maths.isClockWise(positions(loop)) == clockwise) { Collections.reverse(loop); } + loop.remove(loop.size() - 1); + Collections.rotate(loop, -loop.indexOf(Collections.max(loop, Comparator.comparingDouble(this::x)))); + loop.add(loop.get(0)); } - //sort the holes by decreasing maximum x value: the region on the right of a hole is then + //join the holes by decreasing maximum x: the region on the right of a hole is then //bounded by the exterior ring and by the holes already joined to it. - holes.sort(Comparator.comparingDouble((List<Integer> h) -> - h.stream().mapToDouble(this::x).max().orElse(Double.NaN)).reversed()); - - //attach holes to the main geometry + holes.sort(Comparator.comparingDouble((List<Integer> loop) -> x(loop.get(0))).reversed()); for (final List<Integer> loop : holes) { - //remove the duplicated end point of the hole, the cut comes back to the start point instead. - loop.remove(loop.size()-1); - if (loop.isEmpty()) { - continue; - } - //the cut starts from the hole position of maximum x. - int index = 0; - for (int k = 1; k < loop.size(); k++) { - if (x(loop.get(k)) > x(loop.get(index))) { - index = k; - } - } - final int rightmost = loop.get(index); - final List<Integer> cut = new ArrayList<>(loop.size() + 2); - cut.addAll(loop.subList(index, loop.size())); - cut.addAll(loop.subList(0, index + 1)); - - final int visible = visibleVertex(border, rightmost, orientation); - if (visible >= 0) { - //the cut goes from a position of the ring to the hole and back, no position is added. - cut.add(border.get(visible)); - border.addAll(visible + 1, cut); - } else { - //no visible position, which may happen with an invalid polygon. The cut starts on the nearest - //segment of the ring, at a new position whose coordinates and attributes are interpolated. - final Vector<?> m = working.get(rightmost); - final double[] buffer1 = new double[2]; - final double[] buffer2 = new double[2]; - final double[] ratio = new double[2]; - final int[] offset = new int[2]; - Distance.nearest(positions(border), buffer1, List.of(m, m), buffer2, ratio, offset, 0.000000001); - final int bridge = interpolate(border.get(offset[0]), border.get(offset[0] + 1), ratio[0]); - cut.add(0, bridge); - cut.add(bridge); - border.addAll(offset[0] + 1, cut); - } + //the cut goes from a vertex of the ring around the hole, and back to that vertex. + final int vertex = cutStart(border, loop.get(0)); + loop.add(border.get(vertex)); + border.addAll(vertex + 1, loop); } + return border.stream().mapToInt(Integer::intValue).toArray(); + } - //build the ring, with the attributes of the positions of the polygon or interpolated for the added ones. - final DataPointsType type = source.getType(); - final DataPoints ring = GeometryFactory.DEFAULT.createDataPoints(type, border.size()); - for (int i = 0; i < border.size(); i++) { - final int vertex = border.get(i); - final int added = vertex - source.size(); - for (final String name : type.getAttributeNames()) { - ring.setAttribute(i, name, (added < 0) ? source.getAttribute(vertex, name) - : source.interpolateAttribute(name, addedIndices.get(added), addedWeights.get(added))); - } + /** + * Returns the position in the ring of the vertex where the cut toward the given hole vertex starts. + * It is a vertex visible from the hole vertex if there is one or, which may happen with an invalid + * polygon, the nearest end of the segment of the ring nearest to the hole vertex. + * + * @param border the ring, the last position being the same vertex as the first one. + * @param hole the hole vertex of maximum x. + * @return position in {@code border} of the vertex where the cut starts. + */ + private int cutStart(final List<Integer> border, final int hole) { + final int visible = visibleVertex(border, hole); + if (visible >= 0) { + return visible; } - final GeometryFactory factory = GeometryFactory.DEFAULT; - return factory.createPolygon(factory.createLinearRing(ring), null); + final Distance.SegmentPoint nearest = Distance.nearestSegment(positions(border), working[hole]); + return occurrence(border, nearest.segment() + (nearest.ratio() > 0.5 ? 1 : 0), hole); } /** @@ -201,112 +179,87 @@ public final class HoleToExterior { * of the hole vertex, the point reached and that end. Then the reflex vertex closest in angle to the ray * is visible. * - * @param border the ring, the last position being the same vertex as the first one. - * @param hole the hole vertex of maximum x. - * @param orientation 1 if the ring is counterclockwise, -1 if it is clockwise. + * @param border the ring, the last position being the same vertex as the first one. + * @param hole the hole vertex of maximum x. * @return position in {@code border} of a visible vertex, or -1. */ - private int visibleVertex(final List<Integer> border, final int hole, final double orientation) { + private int visibleVertex(final List<Integer> border, final int hole) { final double mx = x(hole), my = y(hole); double ix = Double.POSITIVE_INFINITY; int hit = -1; // position of the vertex reached by the ray, if any. int edge = -1; // position of the start of the segment reached by the ray, if any. for (int i = 1; i < border.size(); i++) { final int a = border.get(i-1), b = border.get(i); - final double ay = y(a), by = y(b); - for (final int k : new int[] {i-1, i}) { - final int v = border.get(k); - if (y(v) == my && x(v) >= mx && x(v) < ix) { - ix = x(v); hit = k; edge = -1; - } - } - if ((ay < my && by > my) || (ay > my && by < my)) { - final double x = x(a) + (my - ay) * (x(b) - x(a)) / (by - ay); + if (y(b) == my && x(b) >= mx && x(b) < ix) { + ix = x(b); hit = i; edge = -1; + } else if ((y(a) < my && y(b) > my) || (y(a) > my && y(b) < my)) { + final double x = x(a) + (my - y(a)) * (x(b) - x(a)) / (y(b) - y(a)); if (x >= mx && x < ix) { ix = x; hit = -1; edge = i-1; } } } if (hit >= 0) { - return occurrence(border, hit, hole, orientation); + return occurrence(border, hit, hole); } if (edge < 0) { return -1; } - int best = (x(border.get(edge + 1)) > x(border.get(edge))) ? edge + 1 : edge; - final int end = border.get(best); - final double[] m = {mx, my}, reached = {ix, my}, p = {x(end), y(end)}; - double bestAngle = angle(hole, end); - double bestDistance = distance(hole, end); - for (int k = 0; k < border.size() - 1; k++) { - final int v = border.get(k); - if ((x(v) == p[0] && y(v) == p[1]) || !isReflex(border, k, orientation) - || !Maths.inTriangle(m, reached, p, new double[] {x(v), y(v)})) { - continue; - } - final double angle = angle(hole, v); - final double distance = distance(hole, v); - if (angle < bestAngle || (angle == bestAngle && distance < bestDistance)) { - best = k; - bestAngle = angle; - bestDistance = distance; - } - } - return occurrence(border, best, hole, orientation); + final int end = (x(border.get(edge + 1)) > x(border.get(edge))) ? edge + 1 : edge; + final Vector<?> m = working[hole], p = working[border.get(end)]; + final Vector<?> reached = m.copy(); + reached.set(0, ix); + final Comparator<Integer> closestToRay = Comparator.comparingDouble((Integer k) -> angle(hole, border.get(k))) + .thenComparingDouble((k) -> Distance.distanceSquare(m, working[border.get(k)])); + final int visible = IntStream.range(0, border.size() - 1).boxed().filter((k) -> { + final Vector<?> v = working[border.get(k)]; + return !v.equals(p) && isReflex(border, k) && Maths.inTriangle(m, reached, p, v); + }).min(closestToRay).orElse(end); + return occurrence(border, visible, hole); } /** * Returns the position in the ring of the vertex at the given position, or of another occurrence of it * such as those added by previous cuts, from which the segment toward the hole vertex enters the polygon. */ - private int occurrence(final List<Integer> border, int position, final int hole, final double orientation) { + private int occurrence(final List<Integer> border, int position, final int hole) { if (position == border.size() - 1) { position = 0; // the last position is the same vertex as the first one. } - final int vertex = border.get(position); - if (inCone(border, position, hole, orientation)) { - return position; - } - for (int k = 0; k < border.size() - 1; k++) { - final int v = border.get(k); - if (x(v) == x(vertex) && y(v) == y(vertex) && inCone(border, k, hole, orientation)) { - return k; + if (!inCone(border, position, hole)) { + final Vector<?> vertex = working[border.get(position)]; + for (int k = 0; k < border.size() - 1; k++) { + if (working[border.get(k)].equals(vertex) && inCone(border, k, hole)) { + return k; + } } } return position; } /** - * Returns whether the segment from the vertex at the given position of the ring toward the given point - * starts inside the polygon, which is between the two segments of the ring at that vertex. - * This is the {@code InCone} test of Joseph O'Rourke in <cite>Computational Geometry in C</cite>. + * Returns whether the segment from the vertex at the given position of the ring toward the given vertex + * starts inside the polygon. */ - private boolean inCone(final List<Integer> border, final int position, final int point, final double orientation) { - final int previous = border.get(position == 0 ? border.size() - 2 : position - 1); - final int vertex = border.get(position); - final int next = border.get(position + 1); - if (orientation * cross(vertex, next, previous) >= 0) { - return orientation * cross(vertex, point, previous) > 0 - && orientation * cross(point, vertex, next) > 0; - } - return !(orientation * cross(vertex, point, next) >= 0 - && orientation * cross(point, vertex, previous) >= 0); + private boolean inCone(final List<Integer> border, final int position, final int point) { + return Maths.inCone(working[previous(border, position)], working[border.get(position)], + working[border.get(position + 1)], working[point], clockwise); } /** - * Returns whether the ring turns in the opposite direction of its orientation at the given position. + * Returns whether the ring turns toward its outside at the given position, + * which is a convex vertex for the ring traversed in the opposite direction. */ - private boolean isReflex(final List<Integer> border, final int position, final double orientation) { - final int previous = border.get(position == 0 ? border.size() - 2 : position - 1); - return orientation * cross(previous, border.get(position), border.get(position + 1)) < 0; + private boolean isReflex(final List<Integer> border, final int position) { + return Maths.isConvex(working[previous(border, position)], working[border.get(position)], + working[border.get(position + 1)], !clockwise); } /** - * Returns the cross product of the vectors from <var>p</var> to <var>q</var> and from <var>p</var> to <var>r</var>, - * positive if <var>r</var> is on the left of the line from <var>p</var> to <var>q</var>. + * Returns the vertex before the given position of the ring, the last position being the same vertex as the first one. */ - private double cross(final int p, final int q, final int r) { - return (x(q) - x(p)) * (y(r) - y(p)) - (y(q) - y(p)) * (x(r) - x(p)); + private static int previous(final List<Integer> border, final int position) { + return border.get((position == 0) ? border.size() - 2 : position - 1); } /** @@ -316,105 +269,35 @@ public final class HoleToExterior { return Math.abs(Math.atan2(y(to) - y(from), x(to) - x(from))); } - /** - * Returns the squared distance between the two given vertices. - */ - private double distance(final int from, final int to) { - final double dx = x(to) - x(from), dy = y(to) - y(from); - return dx*dx + dy*dy; - } - - /** - * Returns twice the signed area of the given ring, positive if it is counterclockwise. - */ - private double area(final List<Integer> ring) { - double area = 0; - for (int i = 0, n = ring.size(); i < n; i++) { - final int a = ring.get(i), b = ring.get((i + 1) % n); - area += x(a) * y(b) - x(b) * y(a); - } - return area; - } - /** * Returns the first working coordinate of the given vertex. */ private double x(final int vertex) { - return working.get(vertex).get(0); + return working[vertex].get(0); } /** * Returns the second working coordinate of the given vertex. */ private double y(final int vertex) { - return working.get(vertex).get(1); - } - - /** - * Returns the indices from {@code start} inclusive to {@code start + count} exclusive, in a mutable list. - */ - private static List<Integer> range(final int start, final int count) { - final List<Integer> indices = new ArrayList<>(count); - for (int i = 0; i < count; i++) { - indices.add(start + i); - } - return indices; + return working[vertex].get(1); } /** * Returns the working positions of the given vertices. */ private List<Vector<?>> positions(final List<Integer> vertices) { - final List<Vector<?>> positions = new ArrayList<>(vertices.size()); - for (final int vertex : vertices) { - positions.add(working.get(vertex)); - } - return positions; + return vertices.stream().<Vector<?>>map((vertex) -> working[vertex]).toList(); } /** - * Adds a vertex on the segment joining the two given vertices, and returns its index. - * Its working position is interpolated linearly, and its attributes will be interpolated - * from the positions of the polygon the two vertices are made of. - * - * @param start index of the vertex at ratio 0. - * @param end index of the vertex at ratio 1. - * @param ratio position of the new vertex between the two vertices, from 0 to 1. - * @return index of the new vertex. + * Returns the indices from {@code start} inclusive to {@code start + count} exclusive, in a mutable list. */ - private int interpolate(final int start, final int end, final double ratio) { - final Vector<?> p0 = working.get(start); - final Vector<?> p1 = working.get(end); - final Vector<?> position = p0.copy(); - for (int d = 0; d < position.getDimension(); d++) { - position.set(d, p0.get(d) + ratio * (p1.get(d) - p0.get(d))); + private static List<Integer> range(final int start, final int count) { + final List<Integer> indices = new ArrayList<>(count); + for (int i = 0; i < count; i++) { + indices.add(start + i); } - final int[] i0 = sourceIndices(start), i1 = sourceIndices(end); - final double[] w0 = sourceWeights(start), w1 = sourceWeights(end); - final int[] indices = Arrays.copyOf(i0, i0.length + i1.length); - System.arraycopy(i1, 0, indices, i0.length, i1.length); - final double[] weights = new double[indices.length]; - for (int i = 0; i < w0.length; i++) weights[i] = (1 - ratio) * w0[i]; - for (int i = 0; i < w1.length; i++) weights[i0.length + i] = ratio * w1[i]; - working.add(position); - addedIndices.add(indices); - addedWeights.add(weights); - return working.size() - 1; - } - - /** - * Returns the indices of the positions of the polygon the given vertex is interpolated from. - */ - private int[] sourceIndices(final int vertex) { - final int added = vertex - source.size(); - return (added < 0) ? new int[] {vertex} : addedIndices.get(added); - } - - /** - * Returns the weights of the positions of the polygon the given vertex is interpolated from. - */ - private double[] sourceWeights(final int vertex) { - final int added = vertex - source.size(); - return (added < 0) ? new double[] {1} : addedWeights.get(added); + return indices; } } diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/HoleToExteriorTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/HoleToExteriorTest.java index 1b8d9802e5..0240a699d2 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/HoleToExteriorTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/HoleToExteriorTest.java @@ -133,6 +133,35 @@ public class HoleToExteriorTest { assertAttribute(ring); } + /** + * The joined ring is given as indices of the positions of the polygon, those of the exterior ring first. + * It ends with the closing position of the exterior ring. + */ + @Test + public void testRingIndices() { + final CoordinateReferenceSystem crs = CommonCRS.WGS84.normalizedGeographic(); + final Polygon polygon = GeometryFactory.DEFAULT.createPolygon(ring(crs, 0,0, 10,0, 10,10, 0,10, 0,0), + List.of(ring(crs, 2,3, 2,5, 4,5, 4,3, 2,3))); + //the cut goes from (10,0) to the hole position (4,5) of index 7, around the hole, and back. + assertArrayEquals(new int[] {0, 1, 7, 8, 5, 6, 7, 1, 2, 3, 4}, HoleToExterior.ring(polygon)); + assertArrayEquals(new int[] {0, 1, 2, 3}, HoleToExterior.ring( + GeometryFactory.DEFAULT.createPolygon(ring(crs, 0,0, 1,0, 1,1, 0,0), null)), "Without hole"); + } + + /** + * A hole outside of the exterior ring makes an invalid polygon, where no vertex is visible from the hole. + * The cut shall start from the nearest end of the nearest segment of the exterior ring, without adding position. + */ + @Test + public void testHoleOutside() { + final CoordinateReferenceSystem crs = CommonCRS.WGS84.normalizedGeographic(); + final Polygon polygon = GeometryFactory.DEFAULT.createPolygon(ring(crs, 0,0, 4,0, 4,4, 0,4, 0,0), + List.of(ring(crs, 6,1, 6,3, 8,3, 8,1, 6,1))); + //the hole position (8,3) of index 7 is nearest to the segment from (4,0) to (4,4), at (4,3). + assertArrayEquals(new int[] {0, 1, 2, 7, 8, 5, 6, 7, 2, 3, 4}, HoleToExterior.ring(polygon)); + assertNoAddedPosition(polygon, HoleToExterior.apply(polygon).getExteriorRing().getDataPoints()); + } + /** * A reflex vertex of the exterior ring hides the end of the segment reached by the ray from the hole, * so the cut shall go to that reflex vertex.
