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 0f3587edc5376a7d23f338b883a84de00eab4a2c Author: jsorel <[email protected]> AuthorDate: Fri Oct 16 14:56:04 2026 +0200 feat(Geometry): separate earclipping and holeremover as two separate classes --- .../sis/geometries/operation/ToPrimitive.java | 63 +- .../operation/triangulate/EarClipping.java | 634 ++++++--------------- .../operation/triangulate/HoleToExterior.java | 420 ++++++++++++++ .../sis/geometries/operation/ToPrimitiveTest.java | 94 ++- .../operation/triangulate/HoleToExteriorTest.java | 178 ++++++ 5 files changed, 899 insertions(+), 490 deletions(-) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/ToPrimitive.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/ToPrimitive.java index 39f7bb5d6a..b1579b30d6 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/ToPrimitive.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/ToPrimitive.java @@ -19,8 +19,10 @@ package org.apache.sis.geometries.operation; import java.util.ArrayList; import java.util.Collections; import java.util.HashMap; +import java.util.LinkedHashSet; import java.util.List; import java.util.Map; +import java.util.Set; import javax.measure.Quantity; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.DataPoints; @@ -38,7 +40,6 @@ import org.apache.sis.geometries.surface.Polygon; import org.apache.sis.geometries.surface.PolyhedralSurface; import org.apache.sis.maths.Array; import org.apache.sis.maths.NDArrays; -import org.apache.sis.maths.Vector; /** @@ -122,30 +123,56 @@ public final class ToPrimitive { } /** - * Triangulates the given polygons into a single mesh of triangles, without index. + * Triangulates the given polygons into a single indexed mesh of triangles. + * The mesh holds the attributes which are common to all the polygons. */ private static MeshPrimitive.Triangles triangles(final CoordinateReferenceSystem crs, final List<Polygon> faces) { - final int dimension = crs.getCoordinateSystem().getDimension(); final EarClipping triangulator = new EarClipping(); - final List<double[]> triangles = new ArrayList<>(); + final List<EarClipping.Triangulation> triangulations = new ArrayList<>(faces.size()); + final Set<String> names = new LinkedHashSet<>(); + names.add(DataPointsType.ATT_POSITION); + int vertexCount = 0, indexCount = 0; for (final Polygon face : faces) { - //the returned list is reused by the next call, so the coordinates are copied now. - for (final Vector<?>[] triangle : triangulator.triangulate(face)) { - final double[] coordinates = new double[3 * dimension]; - for (int v = 0; v < 3; v++) { - for (int d = 0; d < dimension; d++) { - coordinates[v * dimension + d] = triangle[v].get(d); - } - } - triangles.add(coordinates); + final EarClipping.Triangulation triangulation = triangulator.triangulateIndexed(face); + final List<String> faceNames = triangulation.vertices().getType().getAttributeNames(); + if (triangulations.isEmpty()) { + names.addAll(faceNames); + } else { + names.retainAll(faceNames); } + triangulations.add(triangulation); + vertexCount += triangulation.vertices().size(); + indexCount += triangulation.triangles().length; } - final double[] positions = new double[triangles.size() * 3 * dimension]; - for (int i = 0; i < triangles.size(); i++) { - System.arraycopy(triangles.get(i), 0, positions, i * 3 * dimension, 3 * dimension); - } + /* + * The vertices of all faces are appended in a single array for each attribute, + * and the indices of each face are shifted by the number of vertices before it. + */ final MeshPrimitive.Triangles primitive = new MeshPrimitive.Triangles(); - primitive.setPositions(NDArrays.of(crs, positions)); + if (triangulations.isEmpty()) { + primitive.setPositions(NDArrays.of(crs, new double[0])); + return primitive; + } + final DataPointsType type = triangulations.get(0).vertices().getType(); + final int[] index = new int[indexCount]; + for (final String name : names) { + final Array array = NDArrays.of(type.getAttributeSystem(name), type.getAttributeType(name), vertexCount); + int offset = 0; + for (final EarClipping.Triangulation triangulation : triangulations) { + final DataPoints vertices = triangulation.vertices(); + array.set(offset, vertices.getAttributeArray(name), 0, vertices.size()); + offset += vertices.size(); + } + primitive.setAttribute(name, array); + } + int vertexOffset = 0, indexOffset = 0; + for (final EarClipping.Triangulation triangulation : triangulations) { + for (final int i : triangulation.triangles()) { + index[indexOffset++] = vertexOffset + i; + } + vertexOffset += triangulation.vertices().size(); + } + primitive.setIndex(NDArrays.packIntegerDataType(NDArrays.of(1, index))); return primitive; } 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 70023f1917..5f85efcbbd 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,19 +18,18 @@ package org.apache.sis.geometries.operation.triangulate; import java.util.ArrayList; import java.util.Arrays; -import java.util.Comparator; +import java.util.HashMap; import java.util.List; -import java.util.stream.Collectors; -import org.opengis.referencing.crs.CoordinateReferenceSystem; +import java.util.Map; +import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; -import org.apache.sis.geometries.curve.LineString; +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; import org.apache.sis.maths.Maths; import org.apache.sis.maths.NDArrays; import org.apache.sis.maths.Vector; -import static org.apache.sis.maths.Vectors.*; import org.apache.sis.util.ArraysExt; @@ -40,109 +39,137 @@ import org.apache.sis.util.ArraysExt; * Triangulation using ear clipping algorithm. * http://en.wikipedia.org/wiki/Polygon_triangulation * + * <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 + * all the attributes of the vertices are kept.</p> + * * @author Johann Sorel */ public class EarClipping { - private static final class SimplePolygon{ - private LineString outter; - private final List<LineString> inners = new ArrayList(); - /** - * direction of the outter loop - */ - private boolean outterIsClockwise; + /** + * Result of a triangulation: the vertices of the triangles, and the triangles as indices of those vertices. + * + * @param vertices the vertices used by the triangles, with all the attributes of the polygon. + * @param triangles the indices in {@code vertices} of the three corners of each triangle, one triangle after the other. + */ + 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); + } } - private final List<Vector[]> triangles = new ArrayList<>(); + /** + * Type of the angle made by the ring at a vertex, computed when first needed. + */ + private enum Angle { + UNKNOWN, + CONVEX, + CONCAVE + } /** - * A path can be composed of multiple parts. - * Each part is a closed polyline with inner holes. + * Positions of the ring to triangulate, the last one at the same position as the first one. */ - private final List parts = new ArrayList(); + private DataPoints source; - /** number of coordinates, may be inferior to coords length. */ - private int nbCoords; - /** current coordinates to triangulate. */ - private Vector[] coords; /** - * store the type of angle made by the coordinate - * avoid recalculate it each time. - * 0 : not calculated - * 1 : convex - * 2 : concave + * Positions of the ring, 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 int[] coordType; + private Vector<?>[] working; /** - * Order of the coordinates used by the triangulation, as indices of the coordinates of the - * polygon: the triangulation works on the two first ones. The polygon is projected on the - * coordinate plane where it has its largest area, so that a vertical polygon is not degenerate. - * {@code null} if the polygon is triangulated in its (x,y) plane, without reordering. + * Indices of the three vertices of each triangle found so far. */ - private int[] axes; + private final List<int[]> triangles = new ArrayList<>(); + + /** 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. + */ + private Angle[] coordType; //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 + private Vector<?> t1; //previous point + private Vector<?> t2; //point at index + private Vector<?> t3; //next point private void reset(){ //reset values + source = null; + working = null; triangles.clear(); nbCoords = 0; coords = null; coordType = null; - parts.clear(); } /** + * Triangulates the given polygon. * * @param geometry geometry to triangulate - * @return List of Coordinate[] for each triangle. + * @return the vertices of the triangles and their indices. */ - public List<Vector[]> triangulate(Polygon geometry) { - + public Triangulation triangulateIndexed(Polygon geometry) { reset(); - - final SimplePolygon part = new SimplePolygon(); - //copy collection to avoid modifications - part.outter = geometry.getExteriorRing(); - - final int nbHole = geometry.getNumInteriorRing(); - - for(int i=0;i<nbHole;i++){ - final LineString inner = geometry.getInteriorRingN(i); - part.inners.add(inner); + source = HoleToExterior.apply(geometry).getExteriorRing().getDataPoints(); + final int[] axes = workingAxes(source.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(); + return result(); + } - axes = workingAxes(geometry.getExteriorRing().getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION)); - run(part); - - if (axes != null) { - //triangles share their vertices, so each one is restored in a new vector. - for (final Vector[] triangle : triangles) { - for (int i = 0; i < triangle.length; i++) { - triangle[i] = fromWorking(triangle[i]); - } - } + /** + * Triangulates the given polygon into an indexed mesh holding all the attributes of the polygon. + */ + public MeshPrimitive.Triangles toMesh(Polygon polygon) { + final Triangulation triangulation = triangulateIndexed(polygon); + final DataPoints vertices = triangulation.vertices(); + final MeshPrimitive.Triangles mesh = new MeshPrimitive.Triangles(); + for (final String name : vertices.getType().getAttributeNames()) { + mesh.setAttribute(name, vertices.getAttributeArray(name)); } - return triangles; + mesh.setIndex(NDArrays.packIntegerDataType(NDArrays.of(1, triangulation.triangles()))); + return mesh; } /** * Returns the order of the coordinates to use for triangulating a polygon with the given positions, - * or {@code null} if the polygon can be triangulated in its (x,y) plane. 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. + * 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. */ - private static int[] workingAxes(final Array positions) { + 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 null; + return axes; } final double[] c = positions.toArrayDouble(); final int n = c.length / dim; @@ -155,140 +182,86 @@ public class EarClipping { } nx = Math.abs(nx); ny = Math.abs(ny); nz = Math.abs(nz); if (nz >= nx && nz >= ny) { - return null; + return axes; } final int drop = (nx >= ny) ? 0 : 1; - final int[] axes = new int[dim]; int k = 0; for (int i = 0; i < 3; i++) { if (i != drop) axes[k++] = i; } - axes[k++] = drop; - for (int i = 3; i < dim; i++) { - axes[k++] = i; - } + axes[k] = drop; return axes; } /** - * Returns a copy of the given position with its coordinates in the order used by the triangulation. + * 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. */ - private Vector<?> toWorking(final Vector<?> position) { - final Vector<?> copy = position.copy(); - if (axes != null) { - for (int i = 0; i < axes.length; i++) { - copy.set(i, position.get(axes[i])); + private Triangulation result() { + final int count = working.length; + final boolean[] used = new boolean[count]; + for (final int[] triangle : triangles) { + for (final int vertex : triangle) { + used[vertex] = true; } } - return copy; - } - - /** - * Returns a copy of the given position with its coordinates restored in the order of the polygon. - */ - private Vector<?> fromWorking(final Vector<?> position) { - final Vector<?> copy = position.copy(); - for (int i = 0; i < axes.length; i++) { - copy.set(axes[i], position.get(i)); + 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); + } + } } - return copy; + 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)); + } + } + 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); } /** - * TODO not efficient, improve performances. + * Triangulates the ring. */ - public MeshPrimitive.Triangles toMesh(Polygon polygon) { - final List<Vector[]> list = triangulate(polygon); - - final CoordinateReferenceSystem crs = polygon.getCoordinateReferenceSystem(); - final int dimension = crs.getCoordinateSystem().getDimension(); - final Array positions = NDArrays.of(crs, new double[list.size()*3*dimension]); - for (int i = 0, k = 0, n = list.size(); i < n; i++, k+=3) { - final Vector[] t = list.get(i); - positions.set(k+0, t[0]); - positions.set(k+1, t[1]); - positions.set(k+2, t[2]); + private void run() { + if (working.length == 0) { + return; } - - final MeshPrimitive.Triangles triangles = new MeshPrimitive.Triangles(); - triangles.setPositions(positions); - - return triangles; - } - - private void run(SimplePolygon part){ - - //build a single geometry linking inner holes. - //the array stream reuses a single vector for all positions, so each one is copied. - final List<Vector<?>> borderCoords = new ArrayList<>(); - part.outter.getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION).stream(false) - .map(this::toWorking).forEach(borderCoords::add); - - //copy the holes in the same way, in mutable lists since their last point is removed below, - //and sort them by decreasing minimum x value. - final List<List<Vector<?>>> holes = new ArrayList<>(); - for (final LineString inner : part.inners) { - holes.add(inner.getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION).stream(false) - .map(this::toWorking).collect(Collectors.toCollection(ArrayList::new))); - } - holes.sort(Comparator.comparingDouble((List<Vector<?>> h) -> - h.stream().mapToDouble((v) -> v.get(0)).min().orElse(Double.NaN)).reversed()); - - //attach holes to the main geometry - for (final List<Vector<?>> loop : holes) { - //we must find the minimum x coordinate in the inner loop - int index = 0; - Vector min = loop.get(index); - for(int k=1,p=loop.size();k<p;k++){ - Vector candidate = (Vector) loop.get(k); - if (candidate.get(0) < min.get(0)) { - min = candidate; - index = k; - } - } - - //now find the closest point on the outter loop - final List line2 = new ArrayList(); - line2.add(min);line2.add(min); - final double[] buffer1 = new double[2]; - final double[] buffer2 = new double[2]; - final double[] ratio = new double[2]; - final int[] offset = new int[2]; - nearest(borderCoords, buffer1, line2, buffer2, ratio, offset, 0.000000001); - - //the cut starts on the outter loop segment nearest to the hole. The point found above has only - //the two working coordinates, all coordinates are interpolated along that segment. - final Vector<?> start = borderCoords.get(offset[0]); - final Vector<?> end = borderCoords.get(offset[0] + 1); - final double length = Math.hypot(end.get(0) - start.get(0), end.get(1) - start.get(1)); - final double ratioAlong = (length > 0) ? Math.hypot(buffer1[0] - start.get(0), buffer1[1] - start.get(1)) / length : 0; - final Vector<?> bridge = start.copy(); - for (int d = 0; d < bridge.getDimension(); d++) { - bridge.set(d, start.get(d) + ratioAlong * (end.get(d) - start.get(d))); - } - - //make a cut from outter loop to inner loop - int insertIndex = offset[0]+1; - borderCoords.add(insertIndex, bridge); - insertIndex++; - //remove the duplicateion inner loop end point - loop.remove(loop.size()-1); - for(int k=index,p=loop.size();k<p;k++,insertIndex++){ - borderCoords.add(insertIndex,loop.get(k)); - } - for(int k=0;k<=index;k++,insertIndex++){ - borderCoords.add(insertIndex,loop.get(k)); - } - borderCoords.add(insertIndex, bridge.copy()); - + 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 = borderCoords.get(0); - for(int i=1,n=borderCoords.size();i<n;i++){ - Vector candidate = borderCoords.get(i); + 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)){ - borderCoords.remove(i); + border.remove(i); i--; n--; }else{ @@ -296,16 +269,13 @@ public class EarClipping { } } - - - final boolean clockwise = Maths.isClockWise(borderCoords); - - nbCoords = borderCoords.size(); - coordType = new int[nbCoords]; - coords = borderCoords.toArray(Vector[]::new); + 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(!clockwise){ + if(!Maths.isClockWise(Arrays.stream(coords).mapToObj((i) -> working[i]).toList())){ ArraysExt.reverse(coords); } @@ -317,18 +287,18 @@ public class EarClipping { indexPrevious = (index==0) ? nbCoords-2 : index-1 ; indexNext = index+1; - t1 = coords[indexPrevious]; - t2 = coords[index]; - t3 = coords[indexNext]; + t1 = working[coords[indexPrevious]]; + t2 = working[coords[index]]; + t3 = working[coords[indexNext]]; - if(isConvex() && isEar()){ + if(isConvex(index) && isEar()){ //reset angle type, we remove a point so it might change - coordType[indexPrevious]=0; - coordType[indexNext]=0; + coordType[indexPrevious]=Angle.UNKNOWN; + coordType[indexNext]=Angle.UNKNOWN; - triangles.add(new Vector[]{t1,t2,t3}); - removeWithin(coords, index); - removeWithin(coordType, index); + 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 @@ -346,39 +316,26 @@ public class EarClipping { } /** - * @return true if currentsegment is convex - */ - private boolean isConvex(){ - if(coordType[index]==0){ - //calculate angle type - final double side = Maths.lineSide(t1, t3, t2); - coordType[index] = side>0 ? 1 : 2; - } - return coordType[index] == 1; - } - - /** - * @return true if segment is convex + * @return true if the ring is convex at the vertex of given index */ private boolean isConvex(int idx){ - if(coordType[idx]==0){ + if(coordType[idx]==Angle.UNKNOWN){ //calculate angle type - final Vector s1 = coords[(idx==0) ? (nbCoords-2) : (idx-1)]; - final Vector s2 = coords[idx]; - final Vector s3 = coords[idx+1]; + 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 ? 1 : 2; + coordType[idx] = side>0 ? Angle.CONVEX : Angle.CONCAVE; } - return coordType[idx] == 1; + return coordType[idx] == Angle.CONVEX; } /** * Check this segment triangle is an ear. * Does not contain any other point. - * @return - * @throws OperationException + * @return true if no other vertex is inside the triangle */ - private boolean isEar() throws IllegalArgumentException { + 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)}; @@ -388,7 +345,8 @@ public class EarClipping { if (!isConvex(i)){ //check it's not one of the current segment if(i!=index && i!=indexNext && i!=indexPrevious){ - final double[] p = new double[]{coords[i].get(0),coords[i].get(1)}; + 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; @@ -403,250 +361,4 @@ public class EarClipping { return true; } - - - /** - * Remove element at index. - * Move elements after index by -1. - * Last element is set to null - * - * @param values original array - * @param index removed element index - * @return removed object - */ - public static Object removeWithin(Object[] values, int index) { - final Object removedValue = values[index]; - if(index+1 < values.length){ - copy(values, index+1, values.length-index-1, values, index); - } - values[values.length-1] = null; //remove reference - return removedValue; - } - - /** - * Remove element at index. - * Move elements after index by -1. - * Last element is set to zero - * - * @param values original array - * @param index removed element index - * @return removed value - */ - public static int removeWithin(int[] values, int index) { - final int removedValue = values[index]; - if(index+1 < values.length){ - copy(values, index+1, values.length-index-1, values, index); - } - values[values.length-1] = 0; - return removedValue; - } - - /** - * Copy a portion of given source array of length elements from soffset index, to - * the target array beginning to toffset index. - * If needed, the returned array is padded with zeros to obtain the required - * length. - * - * @param source not null - * @param soffset offset to start copying - * @param length number of elemnts to copy - * @param target array where to insert new elements - * @param toffset start insert index - * @return the target array with copied values. - */ - public static Object[] copy(Object[] source, int soffset, int length, Object[] target, int toffset){ - for(int i=soffset;i<soffset+length;i++,toffset++){ - target[toffset] = source[i]; - } - return target; - } - - /** - * Copy a portion of given source array of length elements from soffset index, to - * the target array beginning to toffset index. - * If needed, the returned array is padded with zeros to obtain the required - * length. - * - * @param source not null - * @param soffset offset to start copying - * @param length number of elemnts to copy - * @param target array where to insert new elements - * @param toffset start insert index - * @return the target array with copied values. - */ - public static int[] copy(int[] source, int soffset, int length, int[] target, int toffset){ - for(int i=soffset;i<soffset+length;i++,toffset++){ - target[toffset] = source[i]; - } - return target; - } - - /** - * Copy a portion of given source array of length elements from soffset index, to - * the target array beginning to toffset index. - * If needed, the returned array is padded with zeros to obtain the required - * length. - * - * @param source not null - * @param soffset offset to start copying - * @param length number of elemnts to copy - * @param target array where to insert new elements - * @param toffset start insert index - * @return the target array with copied values. - */ - public static double[] copy(double[] source, int soffset, int length, double[] target, int toffset){ - for(int i=soffset;i<soffset+length;i++,toffset++){ - target[toffset] = source[i]; - } - return target; - } - - /** - * Find nearest points between two lines. - * <br> - * See Distance. - * - * @param line1Coords line 1 coordinates - * @param buffer1 buffer to store nearest point on line 1 - * @param line2Coords line 2 coordinates - * @param buffer2 buffer to store nearest point on line 2 - * @param ratio neareast points ratio between line start and end - * @param offset : will store the segment offset of the nearest points - * @param epsilon tolerance factor - */ - public static void nearest(List<Vector<?>> line1Coords, double[] buffer1, - List<Vector<?>> line2Coords, double[] buffer2, - double[] ratio, int[] offset, double epsilon){ - - double distance = Double.MAX_VALUE; - - final double[] tempRatio = new double[2]; - final double[] tempC1 = new double[buffer1.length]; - final double[] tempC2 = new double[buffer1.length]; - - final int nb1 = line1Coords.size()-1; - final int nb2 = line2Coords.size()-1; - - for(int i=0;i<nb1;i++){ - final Vector<?> s1 = line1Coords.get(i); - final Vector<?> e1 = line1Coords.get(i+1); - for(int k=0;k<nb2;k++){ - final Vector<?> s2 = line2Coords.get(k); - final Vector<?> e2 = line2Coords.get(k+1); - - final double dist = Math.sqrt(distanceSquare( - new double[]{s1.get(0),s1.get(1)}, new double[]{e1.get(0),e1.get(1)}, tempC1, - new double[]{s2.get(0),s2.get(1)}, new double[]{e2.get(0),e2.get(1)}, tempC2, - tempRatio,epsilon)); - if(dist<distance){ - //keep informations - distance = dist; - offset[0] = i; - offset[1] = k; - copy(tempC1, 0, tempC1.length, buffer1, 0); - copy(tempC2, 0, tempC2.length, buffer2, 0); - } - } - } - - } - - /** - * Calculate square distance between two line segments. - * - * Adapted from book : Real-TimeCollision Detection by Christer Ericson - * (ClosestPtSegmentSegment p.149) - * - * @param line1Start line 1 start point - * @param line1End line 1 start point - * @param buffer1 closest point on line 1 - * @param line2Start line 1 start point - * @param line2End line 1 start point - * @param buffer2 closest point on line 2 - * @param ratio size 2 , for each line, - * ratio [0..1] of the closest point position between start and end points. - * @param epsilon tolerance - * @return distance - */ - public static double distanceSquare(double[] line1Start, double[] line1End, double[] buffer1, - double[] line2Start, double[] line2End, double[] buffer2, - double[] ratio, double epsilon) { - final double[] d1 = subtract(line1End,line1Start); // Direction vector of segment S1 - final double[] d2 = subtract(line2End,line2Start); // Direction vector of segment S2 - final double[] r = subtract(line1Start, line2Start); - final double a = dot(d1, d1); // Squared length of segment S1, always nonnegative - final double e = dot(d2, d2); // Squared length of segment S2, always nonnegative - final double f = dot(d2, r); - // Check if either or both segments degenerate into points - if (a <= epsilon && e <= epsilon) { - // Both segments degenerate into points - ratio[0] = 0; - ratio[1] = 0; - copy(line1Start,0,line1Start.length,buffer1,0); - copy(line2Start,0,line2Start.length,buffer2,0); - final double[] t = subtract(buffer1,buffer2); - return dot(t,t); - } - if (a <= epsilon) { - // First segment degenerates into a point - ratio[0] = 0; - ratio[1] = clip( f/e, 0, 1); - } else { - final double c = dot(d1, r); - if (e <= epsilon) { - // Second segment degenerates into a point - ratio[0] = clip( -c/a, 0, 1); - ratio[1] = 0; - } else { - // The general nondegenerate case starts here - final double b = dot(d1, d2); - final double denom = a*e-b*b; // Always nonnegative - // If segments not parallel, compute closest point on L1 to L2 and - // clamp to segment S1. Else pick arbitrary s (here 0) - if (denom != 0d) { - ratio[0] = clip((b*f - c*e) / denom, 0.0f, 1.0f); - } else { - ratio[0] = 0; - } - // Compute point on L2 closest to S1(s) using - // t = Dot((P1 + D1*s) - P2,D2) / Dot(D2,D2) = (b*s + f) / e - ratio[1] = (b*ratio[0] + f) / e; - - //If t in [0,1] done. Else clamp t, recompute s for the new value - //of t using s = Dot((P2 + D2*t) - P1,D1) / Dot(D1,D1)= (t*b - c) / a - //and clamp s to [0, 1] - if(ratio[1] < 0){ - ratio[1] = 0; - ratio[0] = clip(-c / a, 0.0f, 1.0f); - }else if(ratio[1] > 1){ - ratio[1] = 1; - ratio[0] = clip((b - c) / a, 0.0f, 1.0f); - } - - } - } - - add(line1Start, scale(d1, ratio[0]), buffer1); - add(line2Start, scale(d2, ratio[1]), buffer2); - final double[] t = subtract(buffer1,buffer2); - return dot(t,t); - } - - /** - * Clip returns: - * <ul> - * <li>min if the value is lower than min</li> - * <li>value if the value is between min and max</li> - * <li>max if the value is higher than max</li> - * </ul> - * @param val the value - * @param min minimum value - * @param max maximum value - * @return clipped value - */ - public static double clip(double val, double min, double max){ - if(val<min) return min; - if(val>max) return max; - return val; - } } 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 new file mode 100644 index 0000000000..4ee173c257 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/HoleToExterior.java @@ -0,0 +1,420 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.sis.geometries.operation.triangulate; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Collections; +import java.util.Comparator; +import java.util.List; +import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.DataPointsType; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.curve.LinearRing; +import org.apache.sis.geometries.operation.Distance; +import org.apache.sis.geometries.surface.Polygon; +import org.apache.sis.maths.Maths; +import org.apache.sis.maths.Vector; + + +/** + * Joins the holes of a polygon to its exterior ring, which gives a polygon without hole covering the same area. + * + * <p>The holes are first oriented in the opposite direction of the exterior ring, then joined from the + * 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> + * + * <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> + * + * @author Johann Sorel (Geomatys) + */ +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. + */ + private final DataPoints source; + + /** + * Positions of all the vertices, with their coordinates in the order used for comparing them. + */ + private final List<Vector<?>> working = new ArrayList<>(); + + /** + * 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}. + */ + private final List<int[]> addedIndices = new ArrayList<>(); + + /** + * For each vertex added to start a cut, the weights of the positions of the polygon it is interpolated from. + */ + private final List<double[]> addedWeights = new ArrayList<>(); + + /** + * Prepares the joining of the holes of the given polygon. + */ + 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)); + } + } + + /** + * 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. + * + * @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()) { + return polygon; + } + for (final LinearRing hole : polygon.getInteriorRings()) { + if (!hole.getDataPoints().isEmpty()) { + return new HoleToExterior(polygon).join(polygon); + } + } + return polygon; + } + + /** + * Joins the holes of the given polygon, which is the polygon given to the constructor. + */ + private Polygon join(final Polygon polygon) { + //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) { + 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)); + for (final List<Integer> loop : holes) { + if (Math.signum(area(loop)) == orientation) { + Collections.reverse(loop); + } + } + + //sort the holes by decreasing maximum x value: 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 + 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); + } + } + + //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))); + } + } + final GeometryFactory factory = GeometryFactory.DEFAULT; + return factory.createPolygon(factory.createLinearRing(ring), null); + } + + /** + * Returns the position in the ring of a vertex which can be joined to the given hole vertex by a segment + * crossing no edge, or -1 if none is found. This is the search of mutually visible vertices described by + * David Eberly in <cite>Triangulation by Ear Clipping</cite>: a ray is cast from the hole vertex toward + * increasing x, which first reaches the ring at a vertex or inside a segment. In the later case, the end + * of the segment of largest x is visible unless a reflex vertex of the ring is inside the triangle made + * 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. + * @return position in {@code border} of a visible vertex, or -1. + */ + private int visibleVertex(final List<Integer> border, final int hole, final double orientation) { + 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 (x >= mx && x < ix) { + ix = x; hit = -1; edge = i-1; + } + } + } + if (hit >= 0) { + return occurrence(border, hit, hole, orientation); + } + 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); + } + + /** + * 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) { + 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; + } + } + 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>. + */ + 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); + } + + /** + * Returns whether the ring turns in the opposite direction of its orientation at the given position. + */ + 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; + } + + /** + * 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>. + */ + 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)); + } + + /** + * Returns the angle between the increasing x direction and the segment between the two given vertices. + */ + private double angle(final int from, final int to) { + 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); + } + + /** + * 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; + } + + /** + * 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; + } + + /** + * 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. + */ + 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))); + } + 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); + } +} diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/ToPrimitiveTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/ToPrimitiveTest.java index f239dcd651..1bd3e9c2c6 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/ToPrimitiveTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/ToPrimitiveTest.java @@ -16,7 +16,9 @@ */ package org.apache.sis.geometries.operation; +import java.util.LinkedHashMap; import java.util.List; +import java.util.Map; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.Geometry; import org.apache.sis.geometries.GeometryFactory; @@ -136,10 +138,62 @@ public class ToPrimitiveTest { assertTrue(result instanceof MeshPrimitive); final MeshPrimitive primitive = (MeshPrimitive) result; assertEquals(MeshPrimitive.Type.TRIANGLES, primitive.getType()); - assertEquals(6, primitive.getPositions().getLength(), "Expected 2 triangles."); + assertEquals(6, primitive.getIndex().getLength(), "Expected 2 triangles."); assertEquals(4, trianglesArea(primitive), 1E-12); } + /** + * The triangles of a polygon reference its vertices, which keep all their attributes. + * The vertices added to join a hole have their attributes interpolated, so a linear + * function of the position keeps its value at every vertex. + */ + @Test + public void testPolygonAttributes() { + final GeometryFactory f = GeometryFactory.DEFAULT; + final Polygon square = f.createPolygon(attributedRing(0,0, 2,0, 2,2, 0,2, 0,0), null); + final MeshPrimitive mesh = (MeshPrimitive) new GeometryProcessor().toPrimitive(square, null, null); + assertEquals(4, mesh.getPositions().getLength(), "The closing position shall not be a separate vertex."); + assertEquals(6, mesh.getIndex().getLength(), "Expected 2 triangles."); + assertAttribute(mesh); + + final Polygon holed = f.createPolygon(attributedRing(0,0, 10,0, 10,10, 0,10, 0,0), + List.of(attributedRing(2,3, 2,5, 4,5, 4,3, 2,3))); + final MeshPrimitive holedMesh = (MeshPrimitive) new GeometryProcessor().toPrimitive(holed, null, null); + assertEquals(96, trianglesArea(holedMesh), 1E-9); + assertAttribute(holedMesh); + } + + /** + * Name of the attribute holding {@code x + 10*y} at each position, for checking that it is kept. + */ + private static final String ATT_TEST = "TEST"; + + /** + * Creates a ring in the two-dimensional test CRS, with a {@value #ATT_TEST} attribute. + */ + private static LinearRing attributedRing(final double... coordinates) { + final double[] values = new double[coordinates.length / 2]; + for (int i = 0; i < values.length; i++) { + values[i] = coordinates[2*i] + 10 * coordinates[2*i + 1]; + } + final Map<String,Array> attributes = new LinkedHashMap<>(); + attributes.put(DataPointsType.ATT_POSITION, NDArrays.of(CRS2D, coordinates)); + attributes.put(ATT_TEST, NDArrays.of(1, values)); + return GeometryFactory.DEFAULT.createLinearRing(GeometryFactory.DEFAULT.createDataPoints(attributes)); + } + + /** + * Verifies that the {@value #ATT_TEST} attribute of each vertex of the given mesh is {@code x + 10*y}. + */ + private static void assertAttribute(final MeshPrimitive mesh) { + final double[] c = mesh.getPositions().toArrayDouble(); + final double[] values = mesh.getAttribute(ATT_TEST).toArrayDouble(); + assertEquals(c.length / 2, values.length); + for (int i = 0; i < values.length; i++) { + assertEquals(c[2*i] + 10 * c[2*i + 1], values[i], 1E-9); + } + } + /** * Tests the triangulation of a polygon with a hole, which requires to cut the exterior * ring toward the hole. The triangles shall cover the polygon area without the hole. @@ -178,7 +232,7 @@ public class ToPrimitiveTest { final MeshPrimitive primitive = (MeshPrimitive) new GeometryProcessor().toPrimitive(polygon, null, null); assertEquals(MeshPrimitive.Type.TRIANGLES, primitive.getType()); assertEquals(3, primitive.getPositions().getDimension()); - assertEquals(6, primitive.getPositions().getLength(), "Expected 2 triangles."); + assertEquals(6, primitive.getIndex().getLength(), "Expected 2 triangles."); final double[] c = primitive.getPositions().toArrayDouble(); for (int i = 0; i < c.length; i += 3) { assertEquals(c[i+1] == 0 ? 5 : 7, c[i+2], "z shall follow the polygon plane."); @@ -252,7 +306,7 @@ public class ToPrimitiveTest { assertTrue(result instanceof MeshPrimitive.Triangles); final MeshPrimitive primitive = (MeshPrimitive) result; assertEquals(3, primitive.getPositions().getDimension()); - assertEquals(6, primitive.getPositions().getLength(), "Expected 2 triangles."); + assertEquals(6, primitive.getIndex().getLength(), "Expected 2 triangles."); assertEquals(2 * Math.sqrt(5), trianglesArea3D(primitive), 1E-9); } @@ -273,7 +327,7 @@ public class ToPrimitiveTest { final Geometry result = new GeometryProcessor().toPrimitive(f.createPolyhedralSurface(crs, faces), null, null); assertTrue(result instanceof MeshPrimitive.Triangles); final MeshPrimitive primitive = (MeshPrimitive) result; - assertEquals(30, primitive.getPositions().getLength(), "Expected 2 triangles per face."); + assertEquals(30, primitive.getIndex().getLength(), "Expected 2 triangles per face."); assertEquals(20, trianglesArea3D(primitive), 1E-9); } @@ -304,7 +358,7 @@ public class ToPrimitiveTest { assertTrue(result instanceof MeshPrimitive.Triangles); final MeshPrimitive primitive = (MeshPrimitive) result; assertEquals(crs, primitive.getCoordinateReferenceSystem()); - assertEquals(36, primitive.getPositions().getLength(), "Expected 2 triangles per face."); + assertEquals(36, primitive.getIndex().getLength(), "Expected 2 triangles per face."); assertEquals(24, trianglesArea3D(primitive), 1E-9); } @@ -316,7 +370,7 @@ public class ToPrimitiveTest { final CoordinateReferenceSystem crs = CommonCRS.WGS84.geographic3D(); final Polyhedron hollow = GeometryFactory.DEFAULT.createPolyhedron(cube(crs, 0,0,0, 4,4,4), List.of(cube(crs, 1,1,1, 3,3,3))); final MeshPrimitive primitive = (MeshPrimitive) new GeometryProcessor().toPrimitive(hollow, null, null); - assertEquals(72, primitive.getPositions().getLength(), "Expected 2 triangles per face of both shells."); + assertEquals(72, primitive.getIndex().getLength(), "Expected 2 triangles per face of both shells."); assertEquals(96 + 24, trianglesArea3D(primitive), 1E-9); } @@ -349,7 +403,7 @@ public class ToPrimitiveTest { assertTrue(result instanceof MeshPrimitive.Triangles); final MeshPrimitive primitive = (MeshPrimitive) result; assertEquals(3, primitive.getCoordinateReferenceSystem().getCoordinateSystem().getDimension()); - assertEquals(36, primitive.getPositions().getLength(), "Expected 2 triangles per face."); + assertEquals(36, primitive.getIndex().getLength(), "Expected 2 triangles per face."); assertEquals(2*4 + 4*6, trianglesArea3D(primitive), 1E-9); } @@ -369,7 +423,7 @@ public class ToPrimitiveTest { final LineString line = f.createLineString(f.createDataPoints(NDArrays.of(crs, new double[] {0,0, 2,0, 2,1}))); final Geometry walls = new GeometryProcessor().toPrimitive(f.createPrism(line, range, null), null, null); assertTrue(walls instanceof MeshPrimitive.Triangles); - assertEquals(12, ((MeshPrimitive) walls).getPositions().getLength(), "Expected 2 triangles per segment."); + assertEquals(12, ((MeshPrimitive) walls).getIndex().getLength(), "Expected 2 triangles per segment."); assertEquals(2*3 + 1*3, trianglesArea3D((MeshPrimitive) walls), 1E-9); } @@ -382,7 +436,7 @@ public class ToPrimitiveTest { * Returns the sum of the areas of the triangles of the given non-indexed three-dimensional mesh. */ private static double trianglesArea3D(final MeshPrimitive primitive) { - final double[] c = primitive.getPositions().toArrayDouble(); + final double[] c = triangleCoordinates(primitive); double area = 0; for (int i = 0; i + 8 < c.length; i += 9) { final double ux = c[i+3] - c[i], uy = c[i+4] - c[i+1], uz = c[i+5] - c[i+2]; @@ -402,10 +456,28 @@ public class ToPrimitiveTest { } /** - * Returns the sum of the areas of the triangles of the given non-indexed mesh. + * Returns the coordinates of the vertices of the triangles of the given mesh, one triangle after the other. + * If the mesh is indexed, the coordinates of a vertex are repeated for each triangle using it. */ - private static double trianglesArea(final MeshPrimitive primitive) { + private static double[] triangleCoordinates(final MeshPrimitive primitive) { final double[] c = primitive.getPositions().toArrayDouble(); + if (primitive.getIndex() == null) { + return c; + } + final int dim = primitive.getPositions().getDimension(); + final int[] index = primitive.getIndex().toArrayInt(); + final double[] expanded = new double[index.length * dim]; + for (int i = 0; i < index.length; i++) { + System.arraycopy(c, index[i] * dim, expanded, i * dim, dim); + } + return expanded; + } + + /** + * Returns the sum of the areas of the triangles of the given mesh. + */ + private static double trianglesArea(final MeshPrimitive primitive) { + final double[] c = triangleCoordinates(primitive); double area = 0; for (int i = 0; i + 5 < c.length; i += 6) { area += Math.abs((c[i+2] - c[i]) * (c[i+5] - c[i+1]) - (c[i+4] - c[i]) * (c[i+3] - c[i+1])) / 2; 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 new file mode 100644 index 0000000000..1b8d9802e5 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/HoleToExteriorTest.java @@ -0,0 +1,178 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.sis.geometries.operation.triangulate; + +import java.util.Arrays; +import java.util.LinkedHashMap; +import java.util.List; +import java.util.Map; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.DataPointsType; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.curve.LinearRing; +import org.apache.sis.geometries.surface.Polygon; +import org.apache.sis.maths.Array; +import org.apache.sis.maths.NDArrays; +import org.apache.sis.referencing.CommonCRS; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; + + +/** + * Tests {@link HoleToExterior}. + * + * @author Johann Sorel (Geomatys) + */ +public class HoleToExteriorTest { + + /** + * Name of the attribute holding {@code x + 10*y} at each position, for checking that it is kept. + */ + private static final String ATT_TEST = "TEST"; + + /** + * Creates a ring of the given positions in the given reference system, with a {@value #ATT_TEST} attribute. + */ + private static LinearRing ring(final CoordinateReferenceSystem crs, final double... coordinates) { + final int dimension = crs.getCoordinateSystem().getDimension(); + final double[] values = new double[coordinates.length / dimension]; + for (int i = 0; i < values.length; i++) { + values[i] = coordinates[i * dimension] + 10 * coordinates[i * dimension + 1]; + } + final Map<String,Array> attributes = new LinkedHashMap<>(); + attributes.put(DataPointsType.ATT_POSITION, NDArrays.of(crs, coordinates)); + attributes.put(ATT_TEST, NDArrays.of(1, values)); + return GeometryFactory.DEFAULT.createLinearRing(GeometryFactory.DEFAULT.createDataPoints(attributes)); + } + + /** + * Returns the area enclosed by the given ring in the plane of the two given axes. + * A cut followed in both directions does not change the area. + */ + private static double area(final DataPoints ring, final int h0, final int h1) { + double area = 0; + for (int i = 1; i < ring.size(); i++) { + final double[] p = ring.getPosition(i-1).toArrayDouble(); + final double[] q = ring.getPosition(i).toArrayDouble(); + area += p[h0] * q[h1] - q[h0] * p[h1]; + } + return Math.abs(area) / 2; + } + + /** + * Verifies that the {@value #ATT_TEST} attribute of each position of the given ring is {@code x + 10*y}. + */ + private static void assertAttribute(final DataPoints ring) { + for (int i = 0; i < ring.size(); i++) { + final double[] p = ring.getPosition(i).toArrayDouble(); + assertEquals(p[0] + 10 * p[1], ring.getAttribute(i, ATT_TEST).get(0), 1E-9); + } + } + + /** + * A polygon without hole is returned unchanged. + */ + @Test + public void testWithoutHole() { + final Polygon polygon = GeometryFactory.DEFAULT.createPolygon( + ring(CommonCRS.WGS84.normalizedGeographic(), 0,0, 2,0, 2,2, 0,2, 0,0), null); + assertSame(polygon, HoleToExterior.apply(polygon)); + } + + /** + * Verifies that each position of the given ring is a position of the given polygon. + */ + private static void assertNoAddedPosition(final Polygon polygon, final DataPoints ring) { + final DataPoints positions = polygon.getDataPoints(); + next: for (int i = 0; i < ring.size(); i++) { + final double[] p = ring.getPosition(i).toArrayDouble(); + for (int j = 0; j < positions.size(); j++) { + if (Arrays.equals(p, positions.getPosition(j).toArrayDouble())) { + continue next; + } + } + fail("Added position: " + Arrays.toString(p)); + } + } + + /** + * The holes are joined to the exterior ring, which keeps the area of the polygon and the attributes + * of its positions. The cuts join vertices of the polygon, so no position is added. The second hole + * turns in the same direction as the exterior ring, so it shall be reversed. + */ + @Test + public void testWithHoles() { + 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), + ring(crs, 6,6, 8,6, 8,8, 6,8, 6,6))); + final Polygon joined = HoleToExterior.apply(polygon); + assertTrue(joined.getInteriorRings().isEmpty()); + final DataPoints ring = joined.getExteriorRing().getDataPoints(); + assertEquals(5 + 2 * (5 + 1), ring.size(), "Each hole adds its positions, its first position again and the ring vertex again."); + assertArrayEquals(ring.getPosition(0).toArrayDouble(), ring.getPosition(ring.size() - 1).toArrayDouble(), "The ring shall stay closed."); + assertEquals(100 - 4 - 4, area(ring, 0, 1), 1E-9); + assertNoAddedPosition(polygon, ring); + assertAttribute(ring); + } + + /** + * 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. + */ + @Test + public void testHiddenVertex() { + final CoordinateReferenceSystem crs = CommonCRS.WGS84.normalizedGeographic(); + //the ray from the hole vertex (4,5) reaches the segment from (10,0) to (10,10), but the spike + //of the exterior ring up to (7,4) hides the vertex (10,0) from the hole. + final Polygon polygon = GeometryFactory.DEFAULT.createPolygon( + ring(crs, 0,0, 7,0, 7,4, 8,0, 10,0, 10,10, 0,10, 0,0), + List.of(ring(crs, 2,4, 2,6, 4,5, 2,4))); + final DataPoints ring = HoleToExterior.apply(polygon).getExteriorRing().getDataPoints(); + assertEquals(8 + 4 + 1, ring.size(), "The hole adds its positions, its first position again and the ring vertex again."); + assertEquals(area(polygon.getExteriorRing().getDataPoints(), 0, 1) - 2, area(ring, 0, 1), 1E-9); + assertNoAddedPosition(polygon, ring); + boolean found = false; + for (int i = 1; i < ring.size(); i++) { + final double[] p = ring.getPosition(i - 1).toArrayDouble(); + final double[] q = ring.getPosition(i).toArrayDouble(); + if (Arrays.equals(p, new double[] {7, 4}) && Arrays.equals(q, new double[] {4, 5})) { + found = true; + } + } + assertTrue(found, "The cut shall start from the reflex vertex (7,4)."); + } + + /** + * The hole of a vertical polygon is joined in the plane of the polygon. + */ + @Test + public void testVerticalPolygon() { + final CoordinateReferenceSystem crs = CommonCRS.WGS84.geographic3D(); + final Polygon polygon = GeometryFactory.DEFAULT.createPolygon(ring(crs, 0,0,0, 6,0,0, 6,0,4, 0,0,4, 0,0,0), + List.of(ring(crs, 2,0,1, 2,0,3, 4,0,3, 4,0,1, 2,0,1))); + final DataPoints ring = HoleToExterior.apply(polygon).getExteriorRing().getDataPoints(); + assertEquals(24 - 4, area(ring, 0, 2), 1E-9); + for (int i = 0; i < ring.size(); i++) { + assertEquals(0, ring.getPosition(i).get(1), "y shall stay in the wall plane."); + } + assertAttribute(ring); + } +}
