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 13233d6dc2529986d135a574371a350c73b24d0b Author: jsorel <[email protected]> AuthorDate: Wed Oct 21 13:36:45 2026 +0200 feat(Geometry): add Centerline, StraightSkeleton and TerrainExtruder operations --- .../sis/geometries/internal/shared/Planar.java | 262 +++++++ .../operation/extrude/TerrainExtruder.java | 838 +++++++++++++++++++++ .../operation/spatialanalysis2d/Centerline.java | 305 ++++++++ .../spatialanalysis2d/StraightSkeleton.java | 448 +++++++++++ .../operation/extrude/TerrainExtruderTest.java | 338 +++++++++ .../spatialanalysis2d/CenterlineTest.java | 135 ++++ .../spatialanalysis2d/StraightSkeletonTest.java | 150 ++++ 7 files changed, 2476 insertions(+) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Planar.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Planar.java new file mode 100644 index 0000000000..7ef7816082 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Planar.java @@ -0,0 +1,262 @@ +/* + * 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.internal.shared; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.List; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.apache.sis.geometries.DataPointsType; +import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.curve.LinearRing; +import org.apache.sis.geometries.surface.Polygon; +import org.apache.sis.maths.Array; +import org.apache.sis.referencing.CRS; + + +/** + * Helper methods for computations in the (x,y) plane on coordinates stored in primitive arrays, + * two values per point. + * + * @author Johann Sorel (Geomatys) + */ +public final class Planar { + + private Planar() { + } + + /** + * Returns the (x,y) coordinates of the rings of a polygon, without consecutive duplicated points + * and without closing point. The exterior ring is counterclockwise and the holes are clockwise, + * so the polygon interior is always on the left. Degenerated holes are omitted, and the list is + * empty if the exterior ring is degenerated. + * + * @param polygon the polygon. + * @return exterior ring followed by the interior rings. + */ + public static List<double[]> rings(final Polygon polygon) { + final List<double[]> rings = new ArrayList<>(); + for (int i = -1, n = polygon.getNumInteriorRing(); i < n; i++) { + final double[] ring = ring((i < 0) ? polygon.getExteriorRing() : polygon.getInteriorRingN(i), i < 0); + if (ring.length >= 6) rings.add(ring); + else if (i < 0) break; + } + return rings; + } + + private static double[] ring(final LinearRing ring, final boolean counterclockwise) { + final Array positions = ring.getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION); + final int dim = positions.getDimension(); + final double[] coords = positions.toArrayDouble(); + final double[] xy = new double[coords.length / dim * 2]; + for (int i = 0; i < xy.length / 2; i++) { + xy[i*2] = coords[i*dim]; + xy[i*2+1] = coords[i*dim+1]; + } + final int n = removeDuplicates(xy, xy.length); + double area = 0; + for (int i = 0, j = n - 2; i < n; j = i, i += 2) { + area += xy[j] * xy[i+1] - xy[i] * xy[j+1]; + } + if ((area > 0) != counterclockwise) { + for (int i = 0, j = n - 2; i < j; i += 2, j -= 2) { + final double x = xy[i], y = xy[i+1]; + xy[i] = xy[j]; xy[i+1] = xy[j+1]; + xy[j] = x; xy[j+1] = y; + } + } + return Arrays.copyOf(xy, n); + } + + /** + * Removes consecutive duplicated points in place, including last points equal to the first one. + * + * @param xy the coordinates, two values per point. + * @param length number of valid values in the array. + * @return number of valid values after the removal. + */ + public static int removeDuplicates(final double[] xy, final int length) { + int n = 0; + for (int i = 0; i < length; i += 2) { + if (n == 0 || xy[i] != xy[n-2] || xy[i+1] != xy[n-1]) { + xy[n++] = xy[i]; + xy[n++] = xy[i+1]; + } + } + while (n >= 4 && xy[0] == xy[n-2] && xy[1] == xy[n-1]) { + n -= 2; + } + return n; + } + + /** + * Returns whether a point is inside rings, using the even-odd rule. + * A point on a ring gets an arbitrary answer. + * + * @param rings the rings, without closing point. + * @param x point x. + * @param y point y. + * @return whether the point is inside. + */ + public static boolean inside(final List<double[]> rings, final double x, final double y) { + boolean inside = false; + for (double[] ring : rings) { + for (int i = 0, j = ring.length - 2; i < ring.length; j = i, i += 2) { + final double xi = ring[i], yi = ring[i+1], xj = ring[j], yj = ring[j+1]; + if ((yi > y) != (yj > y) && x < (xj - xi) * (y - yi) / (yj - yi) + xi) { + inside = !inside; + } + } + } + return inside; + } + + /** + * Keeps the triangles whose centroid is inside, or outside, the given rings. + * + * @param xy the points coordinates. + * @param triangles three point indices per triangle. + * @param rings the rings, see {@link #inside(List, double, double)}. + * @param inside {@code true} for keeping triangles inside the rings, {@code false} for outside. + * @return the kept triangles. + */ + public static int[] keep(final double[] xy, final int[] triangles, final List<double[]> rings, final boolean inside) { + int n = 0; + for (int i = 0; i < triangles.length; i += 3) { + final int a = triangles[i]*2, b = triangles[i+1]*2, c = triangles[i+2]*2; + if (inside(rings, (xy[a] + xy[b] + xy[c]) / 3, (xy[a+1] + xy[b+1] + xy[c+1]) / 3) == inside) { + System.arraycopy(triangles, i, triangles, n, 3); + n += 3; + } + } + return Arrays.copyOf(triangles, n); + } + + /** + * Returns the number of segments needed for no segment being longer than the given distance. + * + * @param length length to divide. + * @param densification maximum segment length. + * @return number of segments, at least one. + */ + public static int steps(final double length, final double densification) { + return Math.max(1, (int) Math.ceil(length / densification)); + } + + /** + * Returns a key for an edge which does not depend on its direction. + * + * @param a index of an edge end. + * @param b index of the other edge end. + * @return the edge key. + */ + public static long edgeKey(final int a, final int b) { + return (((long) Math.min(a, b)) << 32) | (Math.max(a, b) & 0xFFFFFFFFL); + } + + /** + * Returns the horizontal component of a CRS, or an undefined 2D CRS if none. + * + * @param crs the CRS, may be null. + * @return the 2D CRS. + */ + public static CoordinateReferenceSystem horizontal(CoordinateReferenceSystem crs) { + if (crs != null && crs.getCoordinateSystem().getDimension() != 2) { + crs = CRS.getComponentAt(crs, 0, 2); + } + return (crs != null) ? crs : Geometries.getUndefinedCRS(2); + } + + /** + * Exact (x,y) location used as key in hash maps. + * + * @param x point x. + * @param y point y. + */ + public record XY(double x, double y) { + /** + * Creates a key, replacing -0 by +0 which are not equal for records. + * + * @param x point x. + * @param y point y. + * @return the key. + */ + public static XY of(final double x, final double y) { + return new XY(x + 0.0, y + 0.0); + } + } + + /** + * Growable array of doubles. + */ + public static final class DoubleList { + /** The values, valid up to {@link #size}. */ + public double[] values = new double[64]; + /** Number of valid values. */ + public int size; + + /** + * Appends values. + * + * @param v the values to append. + */ + public void add(final double... v) { + if (size + v.length > values.length) values = Arrays.copyOf(values, Math.max(size * 2, size + v.length)); + System.arraycopy(v, 0, values, size, v.length); + size += v.length; + } + + /** + * Returns a copy of the valid values. + * + * @return the values. + */ + public double[] toArray() { + return Arrays.copyOf(values, size); + } + } + + /** + * Growable array of integers. + */ + public static final class IntList { + /** The values, valid up to {@link #size}. */ + public int[] values = new int[16]; + /** Number of valid values. */ + public int size; + + /** + * Appends values. + * + * @param v the values to append. + */ + public void add(final int... v) { + if (size + v.length > values.length) values = Arrays.copyOf(values, Math.max(size * 2, size + v.length)); + System.arraycopy(v, 0, values, size, v.length); + size += v.length; + } + + /** + * Returns a copy of the valid values. + * + * @return the values. + */ + public int[] toArray() { + return Arrays.copyOf(values, size); + } + } +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/extrude/TerrainExtruder.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/extrude/TerrainExtruder.java new file mode 100644 index 0000000000..a84936f1e3 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/extrude/TerrainExtruder.java @@ -0,0 +1,838 @@ +/* + * 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.extrude; + +import java.util.ArrayList; +import java.util.Arrays; +import java.util.Comparator; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; +import java.util.stream.IntStream; +import org.locationtech.jts.geom.Envelope; +import org.locationtech.jts.index.quadtree.Quadtree; +import org.apache.sis.geometries.DataPointsType; +import org.apache.sis.geometries.GeometryCollection; +import org.apache.sis.geometries.curve.LineString; +import org.apache.sis.geometries.internal.shared.Planar; +import org.apache.sis.geometries.internal.shared.Planar.DoubleList; +import org.apache.sis.geometries.internal.shared.Planar.IntList; +import org.apache.sis.geometries.internal.shared.Planar.XY; +import org.apache.sis.geometries.mesh.MeshPrimitive; +import org.apache.sis.geometries.mesh.MeshPrimitiveIndexes; +import org.apache.sis.geometries.operation.Intersection; +import org.apache.sis.geometries.operation.OperationException; +import org.apache.sis.geometries.operation.ProcessorUtils; +import org.apache.sis.geometries.operation.spatialanalysis2d.Centerline; +import org.apache.sis.geometries.operation.spatialanalysis2d.StraightSkeleton; +import org.apache.sis.geometries.operation.triangulate.delaunay.Delaunay; +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.util.ArgumentChecks; + + +/** + * Embeds flat areas, such as building footprints or roads, in a terrain mesh. + * + * <p>The terrain is a height field: a triangle mesh in which each (x,y) location has a single elevation. + * Each {@linkplain #appendFlatConstraint(Polygon, double, double, boolean, boolean) appended polygon} cuts + * the terrain along its rings and replaces the part inside by a flat top:</p> + * <ul> + * <li>with an absolute height, the top is horizontal at the given elevation;</li> + * <li>relative to the ground, the top has no sideways tilt but follows the terrain along the + * {@linkplain Centerline centerline} of the polygon, like a road climbing a mountain. + * Each point of the top is at the elevation of the ground at its nearest point on the centerline, + * plus the given height.</li> + * </ul> + * <p>The top is joined to the terrain either by vertical walls, or by moving the terrain vertices + * on the polygon rings to the elevation of the top, in which case the neighboring terrain triangles slope + * to the top. When polygons overlap, the last appended one wins: it cuts into the tops of the previous ones.</p> + * + * <p>The top can be replaced by a hipped roof, computed from the {@linkplain StraightSkeleton straight skeleton} + * of the polygon: every face is planar and slopes up from one edge of the polygon with the same pitch, + * from the top elevation at the eaves to the roof elevation at the highest ridge.</p> + * + * <p>Only the positions of the terrain are used, other attributes such as normals or texture coordinates + * are not in the {@linkplain #getResult() result}. The polygons must have the same horizontal CRS than + * the terrain. Parts of a polygon outside the terrain get a top but no wall.</p> + * + * <p>Limitation: where the rings of a polygon cross a wall of a previous polygon, the new walls on each side + * of the crossing stand on different vertices (the previous top and the ground). The surface is closed but + * has a vertical T-junction there.</p> + * + * @author Johann Sorel (Geomatys) + */ +public final class TerrainExtruder { + + /** + * The axes of the (x,y) plane. + */ + private static final int[] XY_AXES = {0, 1}; + + /** + * Ratio on an edge under which a point is considered at the edge end. + */ + private static final double SNAP = 1e-9; + + /** + * The terrain given at construction, for its CRS. + */ + private final MeshPrimitive.Triangles terrain; + + /** + * Vertex positions, three values per vertex. Vertices are never removed, + * unused vertices are dropped when building the result. + */ + private double[] xyz; + private int nbVertex; + + /** + * Surface triangles, counterclockwise when seen from above, three vertex indices per triangle. + * A removed triangle has -1 as first index. Walls are not stored, they are created in the result + * between the edges facing each other on both sides of a cut. + */ + private int[] triangles; + private int nbTriangle; + + /** + * Index of the surface triangles by their (x,y) envelope. + */ + private final Quadtree index = new Quadtree(); + + /** + * Maximum distance between points of the polygon rings and of the centerline, + * used for polygons relative to the ground. + */ + private double densification; + + /** + * Prepares the extrusion of polygons in the given terrain. + * The terrain is copied, it is not modified by this class. + * + * @param terrain the terrain, a 3D height field, not null + * @throws IllegalArgumentException if the terrain positions are not 3D + */ + public TerrainExtruder(final MeshPrimitive.Triangles terrain) { + ArgumentChecks.ensureNonNull("terrain", terrain); + final Array positions = terrain.getPositions(); + if (positions == null || positions.getDimension() != 3) { + throw new IllegalArgumentException("Terrain positions must be 3D"); + } + this.terrain = terrain; + this.xyz = positions.toArrayDouble(); + this.nbVertex = xyz.length / 3; + final Array idx = terrain.getIndex(); + final int[] indices = (idx != null) ? MeshPrimitiveIndexes.toTriangles(idx, terrain.getType()).toArrayInt() + : IntStream.range(0, nbVertex / 3 * 3).toArray(); + this.triangles = new int[Math.max(indices.length, 30)]; + double length = 0; + for (int i = 0; i + 2 < indices.length; i += 3) { + final int a = indices[i], b = indices[i+1], c = indices[i+2]; + addTriangle(a, b, c); + length += distance2D(a, b) + distance2D(b, c) + distance2D(c, a); + } + densification = (length > 0) ? length / (indices.length / 3 * 3) : 1; + } + + /** + * Returns the maximum distance between the points of the rings and of the centerline, + * used for polygons relative to the ground. + * + * @return the densification distance, the mean terrain edge length by default + */ + public double getDensification() { + return densification; + } + + /** + * Sets the maximum distance between the points of the rings and of the centerline, + * used for polygons relative to the ground. Smaller values follow the terrain more closely + * but create more triangles. + * + * @param densification the densification distance, strictly positive + */ + public void setDensification(final double densification) { + ArgumentChecks.ensureStrictlyPositive("densification", densification); + this.densification = densification; + } + + /** + * Appends a flat area, possibly covered by a hipped roof. The terrain is cut along the polygon rings + * and the part inside the polygon is replaced by the top. Tops of previously appended polygons + * are cut the same way. + * + * @param polygon the area footprint, only the (x,y) coordinates are used + * @param height the elevation of the top if absolute, or its height above the ground; + * with a roof, this is the elevation of the eaves + * @param roofHeight elevation of the highest roof ridge, in the same way than {@code height}. + * If equal to {@code height} the top is flat, otherwise it is a hipped roof + * @param relativeToGround whether the top follows the ground along the polygon centerline + * @param walls {@code true} to join the top to the terrain with vertical walls, + * {@code false} to move the terrain vertices on the rings to the top elevation + * @throws OperationException if the polygon CRS differs from the terrain one, or a triangulation + * or the straight skeleton computation failed + */ + public void appendFlatConstraint(final Polygon polygon, final double height, final double roofHeight, + final boolean relativeToGround, final boolean walls) throws OperationException + { + ArgumentChecks.ensureNonNull("polygon", polygon); + ProcessorUtils.ensureSameCRS2D(terrain, polygon); + final List<double[]> rings = Planar.rings(polygon); + if (!rings.isEmpty()) { + new Append(polygon, rings, height, roofHeight, relativeToGround, walls).run(); + } + } + + /** + * Returns the terrain with the appended flat areas. + * + * @return a new mesh, with positions only + */ + public MeshPrimitive.Triangles getResult() { + final IntList faces = new IntList(); + for (int i = 0; i < nbTriangle * 3; i += 3) { + if (triangles[i] >= 0) faces.add(triangles[i], triangles[i+1], triangles[i+2]); + } + /* + * A surface edge without neighbor is on the terrain border or on a cut made by a flat area. + * On a cut, another edge without neighbor goes in the reverse direction at the same location + * on the other side. Those two edges are joined by a vertical wall, with the same orientation. + */ + final long[] edges = new long[faces.size]; + for (int i = 0; i < faces.size; i += 3) { + for (int k = 0; k < 3; k++) { + edges[i+k] = ((long) faces.values[i+k] << 32) | faces.values[i + (k+1) % 3]; + } + } + Arrays.sort(edges); + final Map<List<XY>, int[]> open = new HashMap<>(); + for (long edge : edges) { + final int a = (int) (edge >>> 32), b = (int) edge; + if (Arrays.binarySearch(edges, ((long) b << 32) | a) < 0) { + final int[] other = open.remove(List.of(key(b), key(a))); + if (other != null) { + faces.add(b, a, other[1], b, other[1], other[0]); + } else { + open.put(List.of(key(a), key(b)), new int[] {a, b}); + } + } + } + //keep the triangles not degenerated in 3D, and only their vertices + final int[] remap = new int[nbVertex]; + Arrays.fill(remap, -1); + final DoubleList positions = new DoubleList(); + final IntList indices = new IntList(); + for (int i = 0; i < faces.size; i += 3) { + final int a = faces.values[i]*3, b = faces.values[i+1]*3, c = faces.values[i+2]*3; + final double ux = xyz[b] - xyz[a], uy = xyz[b+1] - xyz[a+1], uz = xyz[b+2] - xyz[a+2]; + final double vx = xyz[c] - xyz[a], vy = xyz[c+1] - xyz[a+1], vz = xyz[c+2] - xyz[a+2]; + if (uy*vz - uz*vy == 0 && uz*vx - ux*vz == 0 && ux*vy - uy*vx == 0) continue; + for (int k = 0; k < 3; k++) { + final int v = faces.values[i+k]; + if (remap[v] < 0) { + remap[v] = positions.size / 3; + positions.add(xyz[v*3], xyz[v*3+1], xyz[v*3+2]); + } + indices.add(remap[v]); + } + } + final MeshPrimitive.Triangles result = new MeshPrimitive.Triangles(); + result.setPositions(NDArrays.of(terrain.getCoordinateReferenceSystem(), positions.toArray())); + result.setIndex(NDArrays.ofUnsigned(1, indices.toArray())); + return result; + } + + // --------------------------------------------------------------------------------------------- + // Working mesh + // --------------------------------------------------------------------------------------------- + + private int addVertex(final double x, final double y, final double z) { + if (nbVertex * 3 == xyz.length) { + xyz = Arrays.copyOf(xyz, Math.max(30, xyz.length * 2)); + } + xyz[nbVertex*3] = x; + xyz[nbVertex*3+1] = y; + xyz[nbVertex*3+2] = z; + return nbVertex++; + } + + /** + * Adds a surface triangle, made counterclockwise. Flat triangles are ignored. + */ + private void addTriangle(final int a, int b, int c) { + final double side = side(a, b, xyz[c*3], xyz[c*3+1]); + if (side == 0) { + return; + } + if (side < 0) { + final int t = b; b = c; c = t; + } + if (nbTriangle * 3 == triangles.length) { + triangles = Arrays.copyOf(triangles, triangles.length * 2); + } + final int id = nbTriangle++; + triangles[id*3] = a; + triangles[id*3+1] = b; + triangles[id*3+2] = c; + index.insert(envelope(id), id); + } + + private void removeTriangle(final int id) { + index.remove(envelope(id), id); + triangles[id*3] = -1; + } + + private Envelope envelope(final int triangle) { + final int a = triangles[triangle*3], b = triangles[triangle*3+1], c = triangles[triangle*3+2]; + final Envelope env = new Envelope(xyz[a*3], xyz[b*3], xyz[a*3+1], xyz[b*3+1]); + env.expandToInclude(xyz[c*3], xyz[c*3+1]); + return env; + } + + /** + * Returns the identifiers of the alive surface triangles whose envelope intersects the given one. + */ + private int[] query(final Envelope env) { + final IntList result = new IntList(); + for (Object item : index.query(env)) { + final int id = (Integer) item; + if (triangles[id*3] >= 0 && env.intersects(envelope(id))) { + result.add(id); + } + } + return result.toArray(); + } + + /** + * Returns the robust side of a point relative to the line from vertex a to vertex b, positive on the left. + */ + private double side(final int a, final int b, final double x, final double y) { + return Maths.lineSide(xyz[a*3], xyz[a*3+1], xyz[b*3], xyz[b*3+1], x, y); + } + + /** + * Returns whether vertex u is before vertex v in (x,y) order. Edges are always processed from their first + * vertex in that order, so that edges at the same location on both sides of a cut give the same points. + */ + private boolean before(final int u, final int v) { + return xyz[u*3] < xyz[v*3] || (xyz[u*3] == xyz[v*3] && xyz[u*3+1] < xyz[v*3+1]); + } + + private double distance2D(final int a, final int b) { + return Math.hypot(xyz[b*3] - xyz[a*3], xyz[b*3+1] - xyz[a*3+1]); + } + + private double[] xy(final int v) { + return new double[] {xyz[v*3], xyz[v*3+1]}; + } + + private XY key(final int v) { + return XY.of(xyz[v*3], xyz[v*3+1]); + } + + /** + * Returns the elevation of the surface at the given location, or NaN if outside the terrain. + */ + private double groundZ(final double x, final double y) { + for (int id : query(new Envelope(x, x, y, y))) { + final int a = triangles[id*3], b = triangles[id*3+1], c = triangles[id*3+2]; + if (side(a, b, x, y) >= 0 && side(b, c, x, y) >= 0 && side(c, a, x, y) >= 0) { + return interpolate(a, b, c, x, y); + } + } + return Double.NaN; + } + + private double interpolate(final int a, final int b, final int c, final double x, final double y) { + final double[] w = Maths.getBarycentricValue2D(xyz[a*3], xyz[a*3+1], xyz[b*3], xyz[b*3+1], xyz[c*3], xyz[c*3+1], x, y); + return w[0] * xyz[a*3+2] + w[1] * xyz[b*3+2] + w[2] * xyz[c*3+2]; + } + + /** + * Receives the points where a segment meets a surface edge. + */ + @FunctionalInterface + private interface Crossing { + /** + * @param t ratio on the segment. + * @param s ratio on the edge from p to q. + * @param p first edge vertex in (x,y) order. + * @param q second edge vertex in (x,y) order. + */ + void at(double t, double s, int p, int q); + } + + /** + * Finds where a segment meets the surface edges. An edge overlapping the segment gives its two ends. + */ + private void crossings(final double x0, final double y0, final double x1, final double y1, final Crossing action) { + final double[] a0 = {x0, y0}, a1 = {x1, y1}; + final Set<Long> seen = new HashSet<>(); + for (int id : query(new Envelope(x0, x1, y0, y1))) { + for (int k = 0; k < 3; k++) { + int p = triangles[id*3+k], q = triangles[id*3+(k+1)%3]; + if (!seen.add(Planar.edgeKey(p, q))) continue; + if (!before(p, q)) { + final int t = p; p = q; q = t; + } + final Intersection.SegmentSegment r = Intersection.segmentSegment2D(a0, a1, xy(p), xy(q), XY_AXES); + if (r.type() != Intersection.SegmentSegment.Type.NONE) { + action.at(r.startOnFirst(), r.startOnSecond(), p, q); + } + if (r.type() == Intersection.SegmentSegment.Type.OVERLAP) { + action.at(r.endOnFirst(), r.endOnSecond(), p, q); + } + } + } + } + + // --------------------------------------------------------------------------------------------- + // Append operation + // --------------------------------------------------------------------------------------------- + + /** + * Vertex inserted on an edge of the terrain. + * + * @param s ratio on the edge, from its first vertex in (x,y) order. + * @param vertex the inserted vertex. + */ + private record Split(double s, int vertex) {} + + /** + * State of the insertion of one polygon. + */ + private final class Append { + + private final Polygon polygon; + private final List<double[]> rings; + private final double height, roofHeight; + private final boolean relative, withWalls; + private final Envelope extent = new Envelope(); + + /** Vertices inserted on the surface edges, by edge key. */ + private final Map<Long, List<Split>> splits = new HashMap<>(); + + /** Vertices inserted inside surface triangles, by triangle identifier. */ + private final Map<Integer, IntList> interiors = new HashMap<>(); + + /** The rings split where they meet the terrain edges, two values per point. */ + private double[] chain; + /** Index of the next point in the same ring, for each chain point. */ + private int[] next; + /** Chain point indices by location. */ + private final Map<XY, IntList> chainAt = new HashMap<>(); + /** Locations of the ring vertices moved on a terrain vertex or edge. */ + private final Map<XY, XY> snapped = new HashMap<>(); + + /** Centerline sampled on the ground, (x,y,z) for each polyline, for relative tops. */ + private final List<double[]> profiles = new ArrayList<>(); + /** Elevation of the ground used when there is no centerline. */ + private double fallbackGround; + + Append(final Polygon polygon, final List<double[]> rings, final double height, final double roofHeight, + final boolean relative, final boolean withWalls) + { + this.polygon = polygon; + this.rings = rings; + this.height = height; + this.roofHeight = roofHeight; + this.relative = relative; + this.withWalls = withWalls; + final double[] exterior = rings.get(0); + for (int i = 0; i < exterior.length; i += 2) { + extent.expandToInclude(exterior[i], exterior[i+1]); + } + } + + void run() throws OperationException { + if (relative) { + //sampled before the surface is modified + buildProfiles(); + } + buildChain(); + cutSurface(); + buildTop(); + } + + // ----------------------------------------------------------------------------------------- + // Relative elevation + // ----------------------------------------------------------------------------------------- + + /** + * Computes the centerline of the polygon and samples the ground elevation along it, + * at each crossing with a surface edge and at the densification distance. + */ + private void buildProfiles() throws OperationException { + final GeometryCollection<LineString> lines = new Centerline(densification).compute(polygon); + for (int l = 0; l < lines.getNumGeometries(); l++) { + final Array positions = lines.getGeometryN(l).getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION); + final int dim = positions.getDimension(); + final double[] c = positions.toArrayDouble(); + final int n = c.length / dim; + final DoubleList profile = new DoubleList(); + for (int i = 0; i < n; i++) { + final double x0 = c[i*dim], y0 = c[i*dim+1]; + final double x1 = c[Math.min(i+1, n-1)*dim], y1 = c[Math.min(i+1, n-1)*dim+1]; + final DoubleList ratios = new DoubleList(); + ratios.add(0); + final int steps = Planar.steps(Math.hypot(x1 - x0, y1 - y0), densification); + for (int k = 1; k < steps; k++) ratios.add((double) k / steps); + if (i < n-1) crossings(x0, y0, x1, y1, (t, s, p, q) -> {if (t > 0 && t < 1) ratios.add(t);}); + final double[] sorted = ratios.toArray(); + Arrays.sort(sorted); + for (double t : sorted) { + final double x = x0 + (x1 - x0) * t, y = y0 + (y1 - y0) * t; + profile.add(x, y, groundZ(x, y)); + } + } + final double[] p = profile.toArray(); + if (fillMissing(p)) { + profiles.add(p); + } + } + //used if there is no centerline or the centerline is outside the terrain + double sum = 0; + int count = 0; + for (double[] ring : rings) { + for (int i = 0; i < ring.length; i += 2) { + final double z = groundZ(ring[i], ring[i+1]); + if (!Double.isNaN(z)) { + sum += z; + count++; + } + } + } + fallbackGround = (count > 0) ? sum / count : 0; + } + + /** + * Replaces the NaN elevations of a profile by the nearest known elevation along it. + * + * @return false if all elevations are NaN + */ + private boolean fillMissing(final double[] p) { + int known = 2; + while (known < p.length && Double.isNaN(p[known])) known += 3; + if (known >= p.length) { + return false; + } + for (int i = 2; i < p.length; i += 3) { + if (Double.isNaN(p[i])) p[i] = (i < known) ? p[known] : p[i-3]; + } + return true; + } + + /** + * Returns the elevation of the top at the given location. + */ + private double topZ(final double x, final double y) { + if (!relative) { + return height; + } + double best = Double.POSITIVE_INFINITY; + double z = fallbackGround; + for (double[] p : profiles) { + final int n = p.length / 3; + for (int i = 0; i < Math.max(1, n - 1); i++) { + final int j = Math.min(i + 1, n - 1); + final double x0 = p[i*3], y0 = p[i*3+1], dx = p[j*3] - x0, dy = p[j*3+1] - y0; + final double len2 = dx*dx + dy*dy; + final double t = (len2 == 0) ? 0 : Maths.clamp(((x - x0) * dx + (y - y0) * dy) / len2, 0, 1); + final double ex = x - x0 - dx * t, ey = y - y0 - dy * t, d = ex*ex + ey*ey; + if (d < best) { + best = d; + z = p[i*3+2] + (p[j*3+2] - p[i*3+2]) * t; + } + } + } + return z + height; + } + + // ----------------------------------------------------------------------------------------- + // Rings split on the terrain + // ----------------------------------------------------------------------------------------- + + /** + * Splits the rings where they cross the surface edges, inserting the crossing points in the edges, + * and inserts the other ring points in the surface. + */ + private void buildChain() { + final DoubleList points = new DoubleList(); + final IntList links = new IntList(); + for (double[] ring : rings) { + final DoubleList ringPoints = new DoubleList(); + for (int i = 0; i < ring.length; i += 2) { + final double x0 = ring[i], y0 = ring[i+1]; + final double x1 = ring[(i+2) % ring.length], y1 = ring[(i+3) % ring.length]; + final List<double[]> events = new ArrayList<>(); + final double[] start = locate(x0, y0); + snapped.put(XY.of(x0, y0), XY.of(start[0], start[1])); + events.add(new double[] {0, start[0], start[1]}); + final int steps = relative ? Planar.steps(Math.hypot(x1 - x0, y1 - y0), densification) : 1; + for (int k = 1; k < steps; k++) { + final double t = (double) k / steps; + final double[] p = locate(x0 + (x1 - x0) * t, y0 + (y1 - y0) * t); + events.add(new double[] {t, p[0], p[1]}); + } + //the ends of the segment are inserted by locate(…) + crossings(x0, y0, x1, y1, (t, s, p, q) -> { + if (t > 0 && t < 1) { + final int v = vertexOnEdge(p, q, s); + events.add(new double[] {t, xyz[v*3], xyz[v*3+1]}); + } + }); + events.sort(Comparator.comparingDouble((e) -> e[0])); + for (double[] e : events) ringPoints.add(e[1], e[2]); + } + ringPoints.size = Planar.removeDuplicates(ringPoints.values, ringPoints.size); + final int offset = points.size / 2, n = ringPoints.size / 2; + points.add(ringPoints.toArray()); + for (int i = 0; i < n; i++) links.add(offset + (i + 1) % n); + } + chain = points.toArray(); + next = links.toArray(); + for (int k = 0; k < next.length; k++) { + chainAt.computeIfAbsent(location(k), (key) -> new IntList()).add(k); + } + } + + private XY location(final int chainPoint) { + return XY.of(chain[chainPoint*2], chain[chainPoint*2+1]); + } + + /** + * Returns the vertex at the given ratio on the edge from p to q, with p before q in (x,y) order, + * inserting it in the edge if needed. + */ + private int vertexOnEdge(final int p, final int q, final double s) { + if (s <= SNAP) return p; + if (s >= 1 - SNAP) return q; + final List<Split> list = splits.computeIfAbsent(Planar.edgeKey(p, q), (key) -> new ArrayList<>()); + for (Split split : list) { + if (Math.abs(split.s() - s) <= SNAP) { + return split.vertex(); + } + } + final int vertex = addVertex(xyz[p*3] + (xyz[q*3] - xyz[p*3]) * s, + xyz[p*3+1] + (xyz[q*3+1] - xyz[p*3+1]) * s, + xyz[p*3+2] + (xyz[q*3+2] - xyz[p*3+2]) * s); + list.add(new Split(s, vertex)); + return vertex; + } + + /** + * Inserts a ring point in the triangles containing it. + * + * @return the point location, moved on the vertex at that location if the point is on an edge + */ + private double[] locate(final double x, final double y) { + final double[] result = {x, y}; + for (int id : query(new Envelope(x, x, y, y))) { + final int a = triangles[id*3], b = triangles[id*3+1], c = triangles[id*3+2]; + final double sab = side(a, b, x, y), sbc = side(b, c, x, y), sca = side(c, a, x, y); + if (sab < 0 || sbc < 0 || sca < 0) { + continue; + } + if (sab > 0 && sbc > 0 && sca > 0) { + interiors.computeIfAbsent(id, (key) -> new IntList()).add(addVertex(x, y, interpolate(a, b, c, x, y))); + continue; + } + final int u = (sab == 0) ? a : (sbc == 0) ? b : c; + final int v = (sab == 0) ? b : (sbc == 0) ? c : a; + final int p = before(u, v) ? u : v, q = (p == u) ? v : u; + final double s = Maths.projectionRatio(xyz[p*3], xyz[p*3+1], xyz[q*3], xyz[q*3+1], x, y); + final int vertex = vertexOnEdge(p, q, Maths.clamp(s, 0, 1)); + result[0] = xyz[vertex*3]; + result[1] = xyz[vertex*3+1]; + } + return result; + } + + /** + * Returns the vertices inserted in the edge from p to q, ordered from p to q. + */ + private IntList along(final int p, final int q) { + final IntList result = new IntList(); + final List<Split> list = splits.get(Planar.edgeKey(p, q)); + if (list != null) { + final Comparator<Split> order = Comparator.comparingDouble(Split::s); + list.sort(before(p, q) ? order : order.reversed()); + list.forEach((split) -> result.add(split.vertex())); + } + return result; + } + + // ----------------------------------------------------------------------------------------- + // Surface cut + // ----------------------------------------------------------------------------------------- + + /** + * Removes the surface inside the polygon. Triangles crossed by the rings are triangulated again + * with the ring segments as constraints, keeping the parts outside the polygon. + */ + private void cutSurface() throws OperationException { + for (int id : query(extent)) { + final int a = triangles[id*3], b = triangles[id*3+1], c = triangles[id*3+2]; + /* + * A ring segment entering the triangle interior leaves it on an edge, where a vertex was + * inserted, or at a vertex, in which case it entered on an edge. Triangles without inserted + * vertices are therefore entirely on one side of the rings. + */ + if (interiors.containsKey(id) || splits.containsKey(Planar.edgeKey(a, b)) + || splits.containsKey(Planar.edgeKey(b, c)) || splits.containsKey(Planar.edgeKey(c, a))) { + cutTriangle(id); + } else if (Planar.inside(rings, (xyz[a*3] + xyz[b*3] + xyz[c*3]) / 3, + (xyz[a*3+1] + xyz[b*3+1] + xyz[c*3+1]) / 3)) { + removeTriangle(id); + } + } + } + + private void cutTriangle(final int id) throws OperationException { + final int a = triangles[id*3], b = triangles[id*3+1], c = triangles[id*3+2]; + final IntList vertices = new IntList(); + final Map<XY,Integer> local = new HashMap<>(); + final IntList constraints = new IntList(); + + //triangle edges, split at the inserted vertices + final int[] corners = {a, b, c}; + for (int k = 0; k < 3; k++) { + final int p = corners[k], q = corners[(k+1) % 3]; + int previous = addLocal(p, vertices, local); + final IntList inner = along(p, q); + inner.add(q); + for (int i = 0; i < inner.size; i++) { + final int current = addLocal(inner.values[i], vertices, local); + constraints.add(previous, current); + previous = current; + } + } + final IntList inside = interiors.get(id); + for (int i = 0; inside != null && i < inside.size; i++) { + addLocal(inside.values[i], vertices, local); + } + + //ring segments inside the triangle + for (int i = 0, n = vertices.size; i < n; i++) { + final IntList points = chainAt.get(key(vertices.values[i])); + for (int j = 0; points != null && j < points.size; j++) { + final int k = points.values[j]; + final Integer other = local.get(location(next[k])); + final double mx = (chain[k*2] + chain[next[k]*2]) / 2; + final double my = (chain[k*2+1] + chain[next[k]*2+1]) / 2; + if (other != null && side(a, b, mx, my) > 0 && side(b, c, mx, my) > 0 && side(c, a, mx, my) > 0) { + constraints.add(i, other); + } + } + } + + final double[] coords = new double[vertices.size * 2]; + for (int i = 0; i < vertices.size; i++) { + coords[i*2] = xyz[vertices.values[i]*3]; + coords[i*2+1] = xyz[vertices.values[i]*3+1]; + } + final int[] parts = Planar.keep(coords, Delaunay.triangulate(coords, constraints.toArray()), rings, false); + removeTriangle(id); + for (int i = 0; i < parts.length; i += 3) { + addTriangle(vertices.values[parts[i]], vertices.values[parts[i+1]], vertices.values[parts[i+2]]); + } + } + + private int addLocal(final int vertex, final IntList vertices, final Map<XY,Integer> local) { + return local.computeIfAbsent(key(vertex), (k) -> { + vertices.add(vertex); + return vertices.size - 1; + }); + } + + // ----------------------------------------------------------------------------------------- + // Top + // ----------------------------------------------------------------------------------------- + + /** + * Triangulates the top on the chain points. For a hipped roof, the straight skeleton nodes are added + * with their time and its arcs as constraints, so that each triangle is inside a single planar face. + * The ring vertices are the skeleton nodes of time 0, they are already chain points. + */ + private void buildTop() throws OperationException { + final Map<XY,Integer> local = new HashMap<>(); + final DoubleList points = new DoubleList(); // x, y, time + final IntList constraints = new IntList(); + for (int k = 0; k < next.length; k++) { + constraints.add(addTop(location(k), 0, local, points), addTop(location(next[k]), 0, local, points)); + } + double roofScale = 0; + if (roofHeight != height) { + final StraightSkeleton skeleton = new StraightSkeleton(polygon); + final double[] nodes = skeleton.getNodes(); + final int[] arcs = skeleton.getArcs(); + for (int i = 0; i < arcs.length; i++) { + final int n = arcs[i] * 3; + final XY p = XY.of(nodes[n], nodes[n+1]); + arcs[i] = addTop(snapped.getOrDefault(p, p), nodes[n+2], local, points); + } + for (int i = 0; i < arcs.length; i += 2) { + if (arcs[i] != arcs[i+1]) constraints.add(arcs[i], arcs[i+1]); + } + if (skeleton.getMaximumTime() > 0) { + roofScale = (roofHeight - height) / skeleton.getMaximumTime(); + } + } + final int nbLocal = points.size / 3; + final double[] xy = new double[nbLocal * 2]; + final double[] z = new double[nbLocal]; + final int[] topVertex = new int[nbLocal]; + for (int l = 0; l < nbLocal; l++) { + xy[l*2] = points.values[l*3]; + xy[l*2+1] = points.values[l*3+1]; + z[l] = topZ(xy[l*2], xy[l*2+1]) + roofScale * points.values[l*3+2]; + topVertex[l] = -1; + } + if (!withWalls) { + //reuse the surface vertices on the rings, moved to the top elevation + for (int id : query(extent)) { + for (int k = 0; k < 3; k++) { + final int v = triangles[id*3+k]; + final Integer l = local.get(key(v)); + if (l != null) { + xyz[v*3+2] = z[l]; + if (topVertex[l] < 0) topVertex[l] = v; + } + } + } + } + for (int l = 0; l < nbLocal; l++) { + if (topVertex[l] < 0) topVertex[l] = addVertex(xy[l*2], xy[l*2+1], z[l]); + } + final int[] tops = Planar.keep(xy, Delaunay.triangulate(xy, constraints.toArray()), rings, true); + for (int i = 0; i < tops.length; i += 3) { + addTriangle(topVertex[tops[i]], topVertex[tops[i+1]], topVertex[tops[i+2]]); + } + } + + private int addTop(final XY p, final double time, final Map<XY,Integer> local, final DoubleList points) { + return local.computeIfAbsent(p, (k) -> { + points.add(k.x(), k.y(), time); + return points.size / 3 - 1; + }); + } + } +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/spatialanalysis2d/Centerline.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/spatialanalysis2d/Centerline.java new file mode 100644 index 0000000000..4e238e3681 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/spatialanalysis2d/Centerline.java @@ -0,0 +1,305 @@ +/* + * 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.spatialanalysis2d; + +import java.util.ArrayList; +import java.util.BitSet; +import java.util.HashMap; +import java.util.List; +import java.util.Map; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.apache.sis.geometries.GeometryCollection; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.curve.LineString; +import org.apache.sis.geometries.internal.shared.Planar; +import org.apache.sis.geometries.internal.shared.Planar.DoubleList; +import org.apache.sis.geometries.internal.shared.Planar.IntList; +import org.apache.sis.geometries.internal.shared.Planar.XY; +import org.apache.sis.geometries.operation.OperationException; +import org.apache.sis.geometries.operation.triangulate.delaunay.Delaunay; +import org.apache.sis.geometries.surface.Polygon; +import org.apache.sis.util.ArgumentChecks; + + +/** + * Computes the centerline of a polygon, the lines running along the middle of its elongated parts, + * for example the axis of a road. + * + * <p>The centerline is computed with the chordal axis transform. The rings are densified so that no segment + * is longer than the densification distance, then the polygon is triangulated with a constrained Delaunay + * triangulation. The edges of the triangles which are not on a ring are the chords of the polygon. + * A triangle with two chords contributes the segment joining their middles, and a triangle with three chords + * contributes the segments joining their middles to its centroid. Triangles with a single chord are at the ends + * of the polygon and contribute nothing.</p> + * + * <p>Branches ending inside the polygon which are shorter than the {@linkplain #getMinimumBranchLength() + * minimum branch length} are removed, they are caused by the bumps of the rings and not by the general shape. + * The remaining segments are merged in polylines, broken at the nodes where branches meet.</p> + * + * <p>Only the first two coordinates are used. The result is in the horizontal component of the polygon CRS.</p> + * + * @author Johann Sorel (Geomatys) + */ +public final class Centerline { + + /** + * Maximum length of a ring segment. + */ + private final double densification; + + /** + * Branches with a free end shorter than this length are removed. + */ + private double minimumBranchLength; + + /** + * Creates a centerline extractor. + * + * @param densification maximum length of a ring segment, also used as default minimum branch length. + * Smaller values give a smoother centerline but more triangles. + */ + public Centerline(final double densification) { + ArgumentChecks.ensureStrictlyPositive("densification", densification); + this.densification = densification; + this.minimumBranchLength = 2 * densification; + } + + /** + * Returns the length under which a branch with a free end is removed. + * + * @return the minimum branch length, twice the densification distance by default + */ + public double getMinimumBranchLength() { + return minimumBranchLength; + } + + /** + * Sets the length under which a branch with a free end is removed. + * + * @param length the minimum branch length, zero to keep all branches + */ + public void setMinimumBranchLength(final double length) { + ArgumentChecks.ensurePositive("length", length); + this.minimumBranchLength = length; + } + + /** + * Computes the centerline of a polygon. + * + * @param polygon the polygon, not null + * @return the polylines of the centerline, empty if the polygon has no elongated part + * @throws OperationException if the triangulation of the polygon failed + */ + public GeometryCollection<LineString> compute(final Polygon polygon) throws OperationException { + ArgumentChecks.ensureNonNull("polygon", polygon); + final CoordinateReferenceSystem crs = Planar.horizontal(polygon.getCoordinateReferenceSystem()); + final GeometryFactory factory = GeometryFactory.DEFAULT; + final List<double[]> lines = computeXY(polygon); + final LineString[] geometries = new LineString[lines.size()]; + for (int i = 0; i < geometries.length; i++) { + geometries[i] = factory.createLineString(factory.createDataPoints(crs, lines.get(i))); + } + return factory.createMultiLineString(crs, geometries); + } + + /** + * Computes the centerline polylines as (x,y) coordinates, two values per point. + */ + private List<double[]> computeXY(final Polygon polygon) throws OperationException { + final List<double[]> rings = Planar.rings(polygon); + + //densified rings, with unique points since rings may touch + final Map<XY, Integer> unique = new HashMap<>(); + final DoubleList points = new DoubleList(); + final IntList segments = new IntList(); + for (double[] ring : rings) { + final IntList ids = new IntList(); + for (int i = 0; i < ring.length; i += 2) { + final double x0 = ring[i], y0 = ring[i+1]; + final double x1 = ring[(i+2) % ring.length], y1 = ring[(i+3) % ring.length]; + final int steps = Planar.steps(Math.hypot(x1 - x0, y1 - y0), densification); + for (int k = 0; k < steps; k++) { + final double r = (double) k / steps; + ids.add(unique.computeIfAbsent(XY.of(x0 + (x1 - x0) * r, y0 + (y1 - y0) * r), (p) -> { + points.add(p.x(), p.y()); + return points.size / 2 - 1; + })); + } + } + for (int j = 0; j < ids.size; j++) { + segments.add(ids.values[j], ids.values[(j+1) % ids.size]); + } + } + final double[] xy = points.toArray(); + final int[] triangles = Planar.keep(xy, Delaunay.triangulate(xy, segments.toArray()), rings, true); + + //chords are the edges shared by two triangles, edges on the rings have a single triangle inside + final Map<Long,Integer> edgeCount = new HashMap<>(); + for (int i = 0; i < triangles.length; i += 3) { + for (int k = 0; k < 3; k++) { + edgeCount.merge(Planar.edgeKey(triangles[i+k], triangles[i+(k+1)%3]), 1, Integer::sum); + } + } + + //build the graph of the chordal axis + final Graph graph = new Graph(); + final Map<Long,Integer> chordNodes = new HashMap<>(); + final int[] chords = new int[3]; + for (int i = 0; i < triangles.length; i += 3) { + int nb = 0; + for (int k = 0; k < 3; k++) { + final int a = triangles[i+k], b = triangles[i+(k+1)%3]; + final long key = Planar.edgeKey(a, b); + if (edgeCount.get(key) == 2) { + chords[nb++] = chordNodes.computeIfAbsent(key, + (e) -> graph.addNode((xy[a*2] + xy[b*2]) / 2, (xy[a*2+1] + xy[b*2+1]) / 2)); + } + } + if (nb == 2) { + graph.addEdge(chords[0], chords[1]); + } else if (nb == 3) { + final int a = triangles[i]*2, b = triangles[i+1]*2, c = triangles[i+2]*2; + final int center = graph.addNode((xy[a] + xy[b] + xy[c]) / 3, (xy[a+1] + xy[b+1] + xy[c+1]) / 3); + for (int k = 0; k < 3; k++) { + graph.addEdge(center, chords[k]); + } + } + } + graph.prune(minimumBranchLength); + return graph.polylines(); + } + + /** + * Undirected graph of the centerline segments. + */ + private static final class Graph { + /** Node coordinates, two values per node. */ + private final DoubleList nodes = new DoubleList(); + /** Edge ids of each node. */ + private final List<IntList> incident = new ArrayList<>(); + /** Node ids of each edge, two values per edge. */ + private final IntList edges = new IntList(); + /** Removed edges. */ + private final BitSet removed = new BitSet(); + + int addNode(final double x, final double y) { + nodes.add(x, y); + incident.add(new IntList()); + return incident.size() - 1; + } + + void addEdge(final int a, final int b) { + incident.get(a).add(edges.size / 2); + incident.get(b).add(edges.size / 2); + edges.add(a, b); + } + + private int other(final int edge, final int node) { + final int a = edges.values[edge*2]; + return (a == node) ? edges.values[edge*2+1] : a; + } + + private int degree(final int node) { + int degree = 0; + final IntList list = incident.get(node); + for (int i = 0; i < list.size; i++) { + if (!removed.get(list.values[i])) degree++; + } + return degree; + } + + /** + * Returns an incident edge of the node which is not removed and not the excluded edge, or -1. + */ + private int nextEdge(final int node, final int excluded) { + final IntList list = incident.get(node); + for (int i = 0; i < list.size; i++) { + final int e = list.values[i]; + if (e != excluded && !removed.get(e)) return e; + } + return -1; + } + + private double length(final int edge) { + final int a = edges.values[edge*2]*2, b = edges.values[edge*2+1]*2; + return Math.hypot(nodes.values[a] - nodes.values[b], nodes.values[a+1] - nodes.values[b+1]); + } + + /** + * Removes the branches with a free end which are shorter than the given length + * and attached to a node where at least two other branches meet. + */ + void prune(final double minLength) { + boolean changed = true; + while (changed) { + changed = false; + for (int node = 0, n = incident.size(); node < n; node++) { + if (degree(node) != 1) continue; + final BitSet path = new BitSet(); + double length = 0; + int current = node; + for (int edge = nextEdge(node, -1); edge >= 0; edge = nextEdge(current, edge)) { + path.set(edge); + length += length(edge); + current = other(edge, current); + if (degree(current) != 2 || length >= minLength) break; + } + if (length < minLength && degree(current) >= 3) { + removed.or(path); + changed = true; + } + } + } + } + + /** + * Merges the remaining edges in polylines broken at the nodes where the degree is not 2. + * Edges remaining after the walks from those nodes are closed loops. + */ + List<double[]> polylines() { + final List<double[]> lines = new ArrayList<>(); + final BitSet visited = (BitSet) removed.clone(); + for (int node = 0, n = incident.size(); node < n; node++) { + if (degree(node) == 2) continue; + final IntList list = incident.get(node); + for (int i = 0; i < list.size; i++) { + if (!visited.get(list.values[i])) { + lines.add(walk(node, list.values[i], visited)); + } + } + } + for (int e = visited.nextClearBit(0); e < edges.size / 2; e = visited.nextClearBit(e + 1)) { + lines.add(walk(edges.values[e*2], e, visited)); + } + return lines; + } + + private double[] walk(final int start, int edge, final BitSet visited) { + final DoubleList line = new DoubleList(); + line.add(nodes.values[start*2], nodes.values[start*2+1]); + int current = start; + while (edge >= 0 && !visited.get(edge)) { + visited.set(edge); + current = other(edge, current); + line.add(nodes.values[current*2], nodes.values[current*2+1]); + if (current == start || degree(current) != 2) break; + edge = nextEdge(current, edge); + } + return line.toArray(); + } + } +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/spatialanalysis2d/StraightSkeleton.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/spatialanalysis2d/StraightSkeleton.java new file mode 100644 index 0000000000..af7b3d50ce --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/spatialanalysis2d/StraightSkeleton.java @@ -0,0 +1,448 @@ +/* + * 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.spatialanalysis2d; + +import java.util.ArrayList; +import java.util.HashSet; +import java.util.List; +import java.util.Set; +import org.apache.sis.geometries.internal.shared.Planar; +import org.apache.sis.geometries.internal.shared.Planar.DoubleList; +import org.apache.sis.geometries.internal.shared.Planar.IntList; +import org.apache.sis.geometries.operation.OperationException; +import org.apache.sis.geometries.surface.Polygon; +import org.apache.sis.util.ArgumentChecks; + + +/** + * Computes the straight skeleton of a polygon, which may have holes. + * + * <p>The straight skeleton is traced by the vertices of the polygon while all its edges move inward, + * parallel to themselves, at the same speed. The time at which a skeleton node is reached is its distance + * to the lines of the edges which created it. Seen as heights, the skeleton gives a hipped roof where every + * face is planar and slopes up from one edge of the polygon with the same pitch.</p> + * + * <p>The computation simulates the moving wavefront. At each step all possible events are searched: + * an edge collapsing to a point (edge event) or a reflex vertex reaching an opposite edge (split event), + * and the earliest one is processed. This costs O(n²) per event, so O(n³) for the whole polygon, which is + * fine for footprints of buildings but not for polygons with thousands of vertices. + * Parallel opposite edges meeting at the same time, frequent in buildings, are handled by + * joining the collapsed parts of the wavefront with ridge arcs.</p> + * + * <p>Only the first two coordinates are used.</p> + * + * @author Johann Sorel (Geomatys) + */ +public final class StraightSkeleton { + + /** + * Unit direction and line offset of each input edge. The inward normal is (-dy, dx) + * and the edge line at time t is the points q with normal·q = offset + t. + */ + private double[] dx, dy, offset; + + /** + * Skeleton nodes, three values per node: x, y, time. + */ + private final DoubleList nodes = new DoubleList(); + + /** + * Skeleton arcs, two node indices per arc, and their keys for avoiding duplicates. + */ + private final IntList arcs = new IntList(); + private final Set<Long> arcKeys = new HashSet<>(); + + /** + * Vertices of the wavefront which may still be alive. + */ + private final List<Vertex> vertices = new ArrayList<>(); + + /** + * Distance tolerance relative to the polygon size, current time of the simulation and largest node time. + */ + private double epsilon, now, maxTime; + + /** + * Computes the straight skeleton of the given polygon. + * + * @param polygon the polygon, not null + * @throws OperationException if the computation did not converge + */ + public StraightSkeleton(final Polygon polygon) throws OperationException { + ArgumentChecks.ensureNonNull("polygon", polygon); + final List<double[]> rings = new ArrayList<>(); + for (double[] ring : Planar.rings(polygon)) { + ring = removeCollinear(ring); + if (ring.length >= 6) rings.add(ring); + else if (rings.isEmpty()) return; // Degenerated exterior ring. + } + if (!rings.isEmpty()) { + initialize(rings); + run(); + } + } + + /** + * Returns the skeleton nodes, three values per node: x, y and the time at which the wavefront + * reached the node. The vertices of the polygon come first, with time 0. + * + * @return the nodes, a copy + */ + public double[] getNodes() { + return nodes.toArray(); + } + + /** + * Returns the skeleton arcs, as pairs of indices in the {@linkplain #getNodes() nodes}. + * + * @return the arcs, a copy + */ + public int[] getArcs() { + return arcs.toArray(); + } + + /** + * Returns the largest time of the nodes, which is the half-width of the widest part of the polygon. + * + * @return the largest node time + */ + public double getMaximumTime() { + return maxTime; + } + + /** + * Creates the input edges, the polygon nodes and the initial wavefront. + * The rings are oriented with the polygon interior on their left. + */ + private void initialize(final List<double[]> rings) { + final double[] exterior = rings.get(0); + double minX = Double.POSITIVE_INFINITY, minY = Double.POSITIVE_INFINITY; + double maxX = Double.NEGATIVE_INFINITY, maxY = Double.NEGATIVE_INFINITY; + for (int i = 0; i < exterior.length; i += 2) { + minX = Math.min(minX, exterior[i]); maxX = Math.max(maxX, exterior[i]); + minY = Math.min(minY, exterior[i+1]); maxY = Math.max(maxY, exterior[i+1]); + } + epsilon = 1e-9 * Math.max(Math.hypot(maxX - minX, maxY - minY), Double.MIN_NORMAL); + final int nbEdge = rings.stream().mapToInt((ring) -> ring.length / 2).sum(); + dx = new double[nbEdge]; + dy = new double[nbEdge]; + offset = new double[nbEdge]; + int first = 0; + for (double[] ring : rings) { + final int n = ring.length / 2; + for (int i = 0; i < n; i++) { + final double x0 = ring[i*2], y0 = ring[i*2+1]; + final double ux = ring[(i+1) % n * 2] - x0, uy = ring[(i+1) % n * 2 + 1] - y0, len = Math.hypot(ux, uy); + dx[first + i] = ux / len; + dy[first + i] = uy / len; + offset[first + i] = dx[first + i] * y0 - dy[first + i] * x0; + } + final Vertex[] loop = new Vertex[n]; + for (int i = 0; i < n; i++) { + loop[i] = create(ring[i*2], ring[i*2+1], first + (i + n - 1) % n, first + i, addNode(ring[i*2], ring[i*2+1])); + } + for (int i = 0; i < n; i++) { + link(loop[i], loop[(i+1) % n]); + } + first += n; + } + } + + /** + * Searches and processes the events until the wavefront vanished. + */ + private void run() throws OperationException { + final int limit = 20 * vertices.size() * vertices.size() + 100; + for (int iteration = 0; ; iteration++) { + if (iteration > limit) { + throw new OperationException("Straight skeleton computation did not converge."); + } + vertices.removeIf((v) -> !v.alive); + if (vertices.isEmpty()) { + return; + } + //search the earliest event, edge events first on equal times + double best = Double.POSITIVE_INFINITY; + Vertex bestA = null, bestB = null; + boolean bestSplit = false; + for (Vertex a : vertices) { + final double t = edgeEvent(a, a.next); + if (t < best - epsilon) { + best = t; bestA = a; bestB = a.next; + } + } + for (Vertex v : vertices) { + if (!v.reflex || v.antiparallel) continue; + for (Vertex p : vertices) { + if (p == v || p.next == v) continue; + final double t = splitEvent(v, p, p.next); + if (t < best - epsilon) { + best = t; bestA = v; bestB = p; bestSplit = true; + } + } + } + if (bestA == null) { + //no more event: remaining loops are degenerated, join them + for (Vertex v : vertices) { + if (v.alive) collapse(v); + } + return; + } + now = Math.max(now, best); + if (bestSplit) { + processSplit(bestA, bestB); + } else { + resolve(merge(bestA, bestB, (bestA.x(now) + bestB.x(now)) / 2, (bestA.y(now) + bestB.y(now)) / 2)); + } + } + } + + /** + * Returns the time when vertices a and b, consecutive on the wavefront, meet. Infinity if never. + */ + private double edgeEvent(final Vertex a, final Vertex b) { + if (a == b || a.antiparallel || b.antiparallel) return Double.POSITIVE_INFINITY; + final int e = a.out; + final double gap = dx[e] * (b.x(now) - a.x(now)) + dy[e] * (b.y(now) - a.y(now)); + if (gap <= epsilon) return now; + final double approach = dx[e] * (a.vx - b.vx) + dy[e] * (a.vy - b.vy); + return (approach > 0) ? now + gap / approach : Double.POSITIVE_INFINITY; + } + + /** + * Returns the time when the reflex vertex v reaches the wavefront edge from p to q. Infinity if never. + */ + private double splitEvent(final Vertex v, final Vertex p, final Vertex q) { + final int e = p.out; + if (e == v.in || e == v.out) return Double.POSITIVE_INFINITY; + final double vx = v.x(now), vy = v.y(now), nx = -dy[e], ny = dx[e]; + //distance of v in front of the moving edge line, and speed at which they approach + final double ahead = nx * vx + ny * vy - (offset[e] + now); + final double closing = 1 - (nx * v.vx + ny * v.vy); + if (ahead < -epsilon || closing <= 0) return Double.POSITIVE_INFINITY; + final double tau = Math.max(0, ahead) / closing, t = now + tau; + //the hit point must be on the edge at that time + final double sh = dx[e] * (vx + v.vx * tau) + dy[e] * (vy + v.vy * tau); + final double sp = dx[e] * p.x(t) + dy[e] * p.y(t); + final double sq = dx[e] * q.x(t) + dy[e] * q.y(t); + return (sh >= sp - epsilon && sh <= sq + epsilon) ? t : Double.POSITIVE_INFINITY; + } + + /** + * The reflex vertex v reaches the wavefront edge from p to p.next: the wavefront is split + * in two loops, or two loops are merged if the edge is in another loop. + */ + private void processSplit(final Vertex v, final Vertex p) { + final double x = v.x(now), y = v.y(now); + final int node = addNode(x, y); + addArc(v.node, node); + v.alive = false; + final Vertex w1 = create(x, y, v.in, p.out, node); + final Vertex w2 = create(x, y, p.out, v.out, node); + final Vertex q = p.next; + link(v.prev, w1); + link(w1, q); + link(p, w2); + link(w2, v.next); + resolve(w1); + resolve(w2); + } + + /** + * Replaces the consecutive vertices from first to last by a single vertex at the given location, + * joined to them by arcs. If those vertices are the whole loop, the loop is removed. + * + * @return the new vertex, or null if the loop has been removed + */ + private Vertex merge(final Vertex first, final Vertex last, final double x, final double y) { + final int node = addNode(x, y); + for (Vertex v = first; ; v = v.next) { + addArc(v.node, node); + v.alive = false; + if (v == last) break; + } + if (last.next == first) { + return null; + } + final Vertex m = create(x, y, first.in, last.out, node); + link(first.prev, m); + link(m, last.next); + return m; + } + + /** + * Cleans the loop containing the given vertex after an event: loops without area are joined by arcs + * and removed, and degenerated vertices between two opposite parallel edges are moved along + * the collapsed part of the wavefront up to the nearest of their neighbors, and merged with it. + */ + private void resolve(Vertex start) { + while (start != null && start.alive) { + final List<Vertex> loop = loop(start); + if (isFlat(loop)) { + collapse(start); + return; + } + final Vertex w = loop.stream().filter((v) -> v.antiparallel).findFirst().orElse(null); + if (w == null) { + return; + } + final Vertex p = w.prev, n = w.next; + final double px = p.x(now), py = p.y(now), nx = n.x(now), ny = n.y(now); + final double distP = dx[w.out] * (px - w.x) + dy[w.out] * (py - w.y); + final double distN = dx[w.out] * (nx - w.x) + dy[w.out] * (ny - w.y); + if (Math.abs(distP - distN) <= epsilon || p.prev == n) { + start = merge(p, n, (px + nx) / 2, (py + ny) / 2); + } else if (distN < distP) { + start = merge(w, n, nx, ny); + } else { + start = merge(p, w, px, py); + } + } + } + + /** + * Removes the loop containing the given vertex, joining its vertices by arcs at their current location. + */ + private void collapse(final Vertex start) { + final List<Vertex> loop = loop(start); + final int[] current = new int[loop.size()]; + for (int i = 0; i < current.length; i++) { + final Vertex v = loop.get(i); + current[i] = addNode(v.x(now), v.y(now)); + addArc(v.node, current[i]); + v.alive = false; + } + for (int i = 0; i < current.length; i++) { + addArc(current[i], current[(i + 1) % current.length]); + } + } + + private static List<Vertex> loop(final Vertex start) { + final List<Vertex> loop = new ArrayList<>(); + Vertex v = start; + do { + loop.add(v); + v = v.next; + } while (v != start); + return loop; + } + + /** + * Returns whether the area of the loop is negligible compared to its perimeter. + */ + private boolean isFlat(final List<Vertex> loop) { + double area = 0, perimeter = 0; + for (int i = 0, n = loop.size(); i < n; i++) { + final Vertex a = loop.get(i), b = loop.get((i + 1) % n); + final double ax = a.x(now), ay = a.y(now), bx = b.x(now), by = b.y(now); + area += ax * by - bx * ay; + perimeter += Math.hypot(bx - ax, by - ay); + } + return Math.abs(area) / 2 <= epsilon * perimeter; + } + + private Vertex create(final double x, final double y, final int in, final int out, final int node) { + final Vertex w = new Vertex(); + w.x = x; + w.y = y; + w.t = now; + w.in = in; + w.out = out; + w.node = node; + //velocity v such as v·n_in = 1 and v·n_out = 1 + final double nix = -dy[in], niy = dx[in], nox = -dy[out], noy = dx[out]; + final double dot = nix * nox + niy * noy; + w.antiparallel = (1 + dot) < 1e-12; + if (!w.antiparallel) { + w.vx = (nix + nox) / (1 + dot); + w.vy = (niy + noy) / (1 + dot); + } + //turning right while the interior is on the left + w.reflex = (dx[in] * dy[out] - dy[in] * dx[out]) < 0; + vertices.add(w); + return w; + } + + private static void link(final Vertex a, final Vertex b) { + a.next = b; + b.prev = a; + } + + /** + * Returns the node at the given location and current time, creating it if none is closer than the tolerance. + */ + private int addNode(final double x, final double y) { + for (int i = 0; i < nodes.size; i += 3) { + if (Math.abs(nodes.values[i] - x) <= epsilon && Math.abs(nodes.values[i+1] - y) <= epsilon) { + return i / 3; + } + } + nodes.add(x, y, now); + maxTime = Math.max(maxTime, now); + return nodes.size / 3 - 1; + } + + private void addArc(final int a, final int b) { + if (a != b && arcKeys.add(Planar.edgeKey(a, b))) { + arcs.add(a, b); + } + } + + /** + * Removes the vertices aligned with their neighbors. + */ + private static double[] removeCollinear(final double[] ring) { + final DoubleList result = new DoubleList(); + final int n = ring.length / 2; + for (int i = 0; i < n; i++) { + final int p = (i + n - 1) % n * 2, c = i * 2, q = (i + 1) % n * 2; + if ((ring[c] - ring[p]) * (ring[q+1] - ring[c+1]) - (ring[c+1] - ring[p+1]) * (ring[q] - ring[c]) != 0) { + result.add(ring[c], ring[c+1]); + } + } + return result.toArray(); + } + + /** + * A vertex of the moving wavefront. + */ + private static final class Vertex { + /** Position at the creation time. */ + double x, y, t; + /** Velocity. */ + double vx, vy; + /** Incoming and outgoing input edges. */ + int in, out; + /** Skeleton node where the vertex started. */ + int node; + /** Neighbors on the wavefront. */ + Vertex prev, next; + /** Whether the vertex turns right, it may split an opposite edge. */ + boolean reflex; + /** Whether the incoming and outgoing edges are opposite, the vertex is on a collapsed part. */ + boolean antiparallel; + boolean alive = true; + + double x(final double time) { + return x + vx * (time - t); + } + + double y(final double time) { + return y + vy * (time - t); + } + } +} diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/extrude/TerrainExtruderTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/extrude/TerrainExtruderTest.java new file mode 100644 index 0000000000..abaa39b025 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/extrude/TerrainExtruderTest.java @@ -0,0 +1,338 @@ +/* + * 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.extrude; + +import java.util.ArrayList; +import java.util.HashMap; +import java.util.HashSet; +import java.util.List; +import java.util.Map; +import java.util.Set; +import java.util.function.DoubleBinaryOperator; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.curve.LinearRing; +import org.apache.sis.geometries.mesh.MeshPrimitive; +import org.apache.sis.geometries.surface.Polygon; +import org.apache.sis.maths.Maths; +import org.apache.sis.maths.NDArrays; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; + + +/** + * Tests {@link TerrainExtruder}. + * + * <p>The terrain is a 10×10 square made of two triangles, sloped along x with z = x. + * Its diagonal goes from (0,0) to (10,10).</p> + * <pre> + * 10 ■-----------■ + * | ⟋ | + * | ⟋ | + * | ⟋ | + * | ⟋ | + * | ⟋ | + * 0 ■-----------■ + * 0 10 + * </pre> + * + * @author Johann Sorel (Geomatys) + */ +public class TerrainExtruderTest { + + private static final CoordinateReferenceSystem CRS3D = Geometries.getUndefinedCRS(3); + + /** + * Footprint crossing the terrain diagonal at (3,3) and (5,5). + */ + private static final double[] FOOTPRINT = {2,3, 5,3, 5,6, 2,6}; + + /** + * Absolute top joined to the terrain by walls. + */ + @Test + public void testAbsoluteWithWalls() { + final Mesh mesh = extrude(20, 20, polygon(FOOTPRINT)); + for (double[] p : new double[][] {{2,3},{5,3},{5,6},{2,6},{3,3},{5,5}}) { + assertTrue(mesh.hasVertex(p[0], p[1], 20), "Missing top vertex at " + p[0] + " " + p[1]); + } + mesh.assertElevations((x, y) -> inside(FOOTPRINT, x, y) ? 20 : x); + assertTrue(mesh.countVertical() > 0, "Walls expected."); + mesh.assertWatertight(); + } + + /** + * Absolute top joined to the terrain by moving the terrain vertices on the rings. + */ + @Test + public void testAbsoluteFlattened() { + final TerrainExtruder extruder = new TerrainExtruder(terrain()); + extruder.appendFlatConstraint(polygon(FOOTPRINT), 20, 20, false, false); + final Mesh mesh = new Mesh(extruder.getResult()); + assertEquals(0, mesh.countVertical(), "No wall expected."); + final Set<List<Double>> xy = new HashSet<>(); + for (int i = 0; i < mesh.positions.length; i += 3) { + assertTrue(xy.add(List.of(mesh.positions[i], mesh.positions[i+1])), "Duplicated vertex location."); + if (onRing(FOOTPRINT, mesh.positions[i], mesh.positions[i+1])) { + assertEquals(20, mesh.positions[i+2], 1E-9); + } + } + mesh.assertWatertight(); + } + + /** + * Top relative to the ground, along the slope: flat cross-section, climbing with the terrain. + */ + @Test + public void testRelative() { + final TerrainExtruder extruder = new TerrainExtruder(terrain()); + extruder.setDensification(0.5); + extruder.appendFlatConstraint(polygon(new double[] {1,4, 9,4, 9,6, 1,6}), 2, 2, true, true); + final Mesh mesh = new Mesh(extruder.getResult()); + //top vertices on the long sides, by x + final Map<Double,Double> south = new HashMap<>(), north = new HashMap<>(); + for (int i = 0; i < mesh.positions.length; i += 3) { + final double x = mesh.positions[i], y = mesh.positions[i+1], z = mesh.positions[i+2]; + if (x < 3 || x > 7 || z <= x + 1) continue; // Skip ground vertices and road ends. + if (y == 4) south.put(x, z); + if (y == 6) north.put(x, z); + assertEquals(x + 2, z, 1E-6, "Top shall be 2 above the ground at the centerline."); + } + assertFalse(south.isEmpty()); + south.forEach((x, z) -> { + if (north.containsKey(x)) assertEquals(z, north.get(x), 1E-6, "Cross-section shall be flat."); + }); + mesh.assertWatertight(); + } + + /** + * Two overlapping footprints, the last appended one wins. + * Not watertight: where the second rings cross a wall of the first footprint, the new walls on each + * side have different bottoms (top of the first footprint or ground), leaving a vertical T-junction. + */ + @Test + public void testOverlap() { + final double[] first = {2,2, 6,2, 6,6, 2,6}; + final double[] second = {4,4, 8,4, 8,8, 4,8}; + final TerrainExtruder extruder = new TerrainExtruder(terrain()); + extruder.appendFlatConstraint(polygon(first), 15, 15, false, true); + extruder.appendFlatConstraint(polygon(second), 25, 25, false, true); + final Mesh mesh = new Mesh(extruder.getResult()); + mesh.assertElevations((x, y) -> inside(second, x, y) ? 25 : inside(first, x, y) ? 15 : x); + for (int t = 0; t < mesh.index.length / 3; t++) { + final double[] c = mesh.centroid(t); + assertFalse(mesh.isVertical(t) && inside(second, c[0], c[1]) && !onRing(second, c[0], c[1]), + "Wall of the first footprint inside the second one."); + } + } + + /** + * Footprint with a courtyard, the terrain stays in the courtyard. + */ + @Test + public void testHole() { + final double[] exterior = {2,2, 8,2, 8,8, 2,8}; + final double[] hole = {4,4, 4,6, 6,6, 6,4}; + final Mesh mesh = extrude(20, 20, polygon(exterior, hole)); + mesh.assertElevations((x, y) -> inside(exterior, x, y) && !inside(hole, x, y) ? 20 : x); + int holeWalls = 0; + for (int t = 0; t < mesh.index.length / 3; t++) { + final double[] c = mesh.centroid(t); + if (mesh.isVertical(t) && onRing(hole, c[0], c[1])) holeWalls++; + } + assertTrue(holeWalls > 0, "Walls expected around the courtyard."); + mesh.assertWatertight(); + } + + /** + * Hipped roof on a rectangle: eaves at 20, ridge from (3,5) to (7,5) at 23. + * Every roof face is planar, its elevation grows with the distance to the footprint edges. + */ + @Test + public void testRoof() { + final double[] footprint = {2,4, 8,4, 8,6, 2,6}; + final Mesh mesh = extrude(20, 23, polygon(footprint)); + assertTrue(mesh.hasVertex(3, 5, 23), "Missing ridge start."); + assertTrue(mesh.hasVertex(7, 5, 23), "Missing ridge end."); + for (double[] p : new double[][] {{2,4},{8,4},{8,6},{2,6},{4,4},{6,6}}) { + assertTrue(mesh.hasVertex(p[0], p[1], 20), "Missing eave vertex at " + p[0] + " " + p[1]); + } + //distance to the nearest edge of the rectangle + mesh.assertElevations((x, y) -> inside(footprint, x, y) + ? 20 + 3 * Math.min(Math.min(x - 2, 8 - x), Math.min(y - 4, 6 - y)) : x); + mesh.assertWatertight(); + } + + /** + * Extrudes a single polygon with walls in the terrain. + */ + private static Mesh extrude(final double height, final double roofHeight, final Polygon polygon) { + final TerrainExtruder extruder = new TerrainExtruder(terrain()); + extruder.appendFlatConstraint(polygon, height, roofHeight, false, true); + return new Mesh(extruder.getResult()); + } + + /** + * Creates the terrain, z = x. + */ + private static MeshPrimitive.Triangles terrain() { + final MeshPrimitive.Triangles terrain = new MeshPrimitive.Triangles(); + terrain.setPositions(NDArrays.of(CRS3D, new double[] {0,0,0, 10,0,10, 10,10,10, 0,10,0})); + terrain.setIndex(NDArrays.ofUnsigned(1, new int[] {0,1,2, 0,2,3})); + return terrain; + } + + /** + * Creates a polygon from unclosed (x,y) rings, the exterior one first. + */ + private static Polygon polygon(final double[] exterior, final double[]... holes) { + final List<LinearRing> interiors = new ArrayList<>(); + for (double[] hole : holes) { + interiors.add(ring(hole)); + } + return GeometryFactory.DEFAULT.createPolygon(ring(exterior), interiors); + } + + private static LinearRing ring(final double[] xy) { + final double[] xyz = new double[xy.length / 2 * 3 + 3]; + for (int i = 0; i < xyz.length / 3; i++) { + xyz[i*3] = xy[i*2 % xy.length]; + xyz[i*3+1] = xy[(i*2 + 1) % xy.length]; + } + return GeometryFactory.DEFAULT.createLinearRing(GeometryFactory.DEFAULT.createDataPoints(CRS3D, xyz)); + } + + /** + * Even-odd point in ring test. + */ + private static boolean inside(final double[] ring, final double x, final double y) { + boolean inside = false; + for (int i = 0, j = ring.length - 2; i < ring.length; j = i, i += 2) { + if ((ring[i+1] > y) != (ring[j+1] > y) && x < (ring[j] - ring[i]) * (y - ring[i+1]) / (ring[j+1] - ring[i+1]) + ring[i]) { + inside = !inside; + } + } + return inside; + } + + /** + * Whether the point is on a segment of the ring, with a small tolerance. + */ + private static boolean onRing(final double[] ring, final double x, final double y) { + for (int i = 0, j = ring.length - 2; i < ring.length; j = i, i += 2) { + if (Maths.distanceSquare(ring[j], ring[j+1], ring[i], ring[i+1], x, y) < 1E-12) { + return true; + } + } + return false; + } + + /** + * Result mesh with helper methods. + */ + private static final class Mesh { + final double[] positions; + final int[] index; + + Mesh(final MeshPrimitive.Triangles mesh) { + positions = mesh.getPositions().toArrayDouble(); + index = mesh.getIndex().toArrayInt(); + } + + double[] vertex(final int triangle, final int k) { + final int v = index[triangle*3 + k]; + return new double[] {positions[v*3], positions[v*3+1], positions[v*3+2]}; + } + + double[] centroid(final int triangle) { + final double[] c = new double[3]; + for (int k = 0; k < 3; k++) { + final double[] p = vertex(triangle, k); + for (int i = 0; i < 3; i++) c[i] += p[i] / 3; + } + return c; + } + + boolean isVertical(final int triangle) { + final double[] a = vertex(triangle, 0), b = vertex(triangle, 1), c = vertex(triangle, 2); + return (b[0] - a[0]) * (c[1] - a[1]) - (b[1] - a[1]) * (c[0] - a[0]) == 0; + } + + int countVertical() { + int count = 0; + for (int t = 0; t < index.length / 3; t++) { + if (isVertical(t)) count++; + } + return count; + } + + boolean hasVertex(final double x, final double y, final double z) { + for (int i = 0; i < positions.length; i += 3) { + if (Math.abs(positions[i] - x) < 1E-9 && Math.abs(positions[i+1] - y) < 1E-9 + && Math.abs(positions[i+2] - z) < 1E-9) { + return true; + } + } + return false; + } + + /** + * Verifies that the centroid of each non-vertical triangle is at the expected elevation. + * Faces are planar, so this also verifies that no face crosses an elevation change. + */ + void assertElevations(final DoubleBinaryOperator expected) { + for (int t = 0; t < index.length / 3; t++) { + if (!isVertical(t)) { + final double[] c = centroid(t); + assertEquals(expected.applyAsDouble(c[0], c[1]), c[2], 1E-9, "Wrong elevation at " + c[0] + " " + c[1]); + } + } + } + + /** + * Verifies that each directed edge is used once and its reverse once, except on the terrain border. + */ + void assertWatertight() { + final Map<List<Integer>, Integer> edges = new HashMap<>(); + for (int i = 0; i < index.length; i += 3) { + for (int k = 0; k < 3; k++) { + edges.merge(List.of(index[i+k], index[i+(k+1)%3]), 1, Integer::sum); + } + } + edges.forEach((edge, count) -> { + final int a = edge.get(0), b = edge.get(1); + assertEquals(1, count, "Edge used more than once in the same direction."); + if (!onBorder(a, b)) { + assertTrue(edges.containsKey(List.of(b, a)), "Edge without neighbor at " + + positions[a*3] + " " + positions[a*3+1] + " " + positions[a*3+2] + " -> " + + positions[b*3] + " " + positions[b*3+1] + " " + positions[b*3+2]); + } + }); + } + + private boolean onBorder(final int a, final int b) { + for (int axis = 0; axis < 2; axis++) { + final double va = positions[a*3+axis], vb = positions[b*3+axis]; + if (va == vb && (va == 0 || va == 10)) return true; + } + return false; + } + } +} diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/spatialanalysis2d/CenterlineTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/spatialanalysis2d/CenterlineTest.java new file mode 100644 index 0000000000..126d725c53 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/spatialanalysis2d/CenterlineTest.java @@ -0,0 +1,135 @@ +/* + * 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.spatialanalysis2d; + +import java.util.ArrayList; +import java.util.List; +import org.apache.sis.geometries.DataPointsType; +import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryCollection; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.curve.LineString; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; + + +/** + * Tests {@link Centerline}. + * + * @author Johann Sorel (Geomatys) + */ +public class CenterlineTest { + + /** + * A long rectangle gives a single polyline along its middle. + * <pre> + * 10 ■----------------------------------■ + * | ============================ | + * 0 ■----------------------------------■ + * 0 100 + * </pre> + */ + @Test + public void testRectangle() { + final List<double[]> lines = compute(0,0, 100,0, 100,10, 0,10, 0,0); + assertEquals(1, lines.size()); + final double[] line = lines.get(0); + double minX = Double.POSITIVE_INFINITY, maxX = Double.NEGATIVE_INFINITY; + for (int i = 0; i < line.length; i += 2) { + if (line[i] > 10 && line[i] < 90) { + //ends may bend toward a corner + assertEquals(5, line[i+1], 1E-9, "Centerline shall be in the middle of the rectangle."); + } + minX = Math.min(minX, line[i]); + maxX = Math.max(maxX, line[i]); + } + assertTrue(minX < 10, "Centerline shall reach the left end."); + assertTrue(maxX > 90, "Centerline shall reach the right end."); + } + + /** + * A L shape gives a single polyline with a bend. + * <pre> + * 100 ■--■ + * | | + * 10 ■--------■ | + * 0 ■-----------■ + * 0 90 100 + * </pre> + */ + @Test + public void testLShape() { + final List<double[]> lines = compute(0,0, 100,0, 100,100, 90,100, 90,10, 0,10, 0,0); + assertEquals(1, lines.size()); + assertTrue(reaches(lines, (x, y) -> x < 50 && Math.abs(y - 5) < 2), "Centerline shall follow the horizontal branch."); + assertTrue(reaches(lines, (x, y) -> y > 50 && Math.abs(x - 95) < 2), "Centerline shall follow the vertical branch."); + } + + /** + * A T shape gives branches reaching the three ends. + * <pre> + * 100 ■-----------------■ + * 90 ■------■ ■------■ + * | | + * 0 ■---■ + * 0 45 55 100 + * </pre> + */ + @Test + public void testTShape() { + final List<double[]> lines = compute(45,0, 55,0, 55,90, 100,90, 100,100, 0,100, 0,90, 45,90, 45,0); + assertTrue(lines.size() >= 3, "Three branches expected."); + assertTrue(reaches(lines, (x, y) -> x < 15), "Centerline shall reach the left end."); + assertTrue(reaches(lines, (x, y) -> x > 85), "Centerline shall reach the right end."); + assertTrue(reaches(lines, (x, y) -> y < 15), "Centerline shall reach the bottom end."); + } + + private interface Condition { + boolean test(double x, double y); + } + + /** + * Returns whether a point of the lines satisfies the condition. + */ + private static boolean reaches(final List<double[]> lines, final Condition condition) { + for (double[] line : lines) { + for (int i = 0; i < line.length; i += 2) { + if (condition.test(line[i], line[i+1])) return true; + } + } + return false; + } + + /** + * Computes the centerline of a polygon without hole, with a densification of 2 and branches of 10, + * and returns the coordinates of each polyline. + */ + private static List<double[]> compute(final double... exterior) { + final GeometryFactory factory = GeometryFactory.DEFAULT; + final Centerline centerline = new Centerline(2); + centerline.setMinimumBranchLength(10); + final GeometryCollection<LineString> result = centerline.compute(factory.createPolygon( + factory.createLinearRing(factory.createDataPoints(Geometries.getUndefinedCRS(2), exterior)), List.of())); + final List<double[]> lines = new ArrayList<>(); + for (int i = 0; i < result.getNumGeometries(); i++) { + lines.add(result.getGeometryN(i).getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION).toArrayDouble()); + } + return lines; + } +} diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/spatialanalysis2d/StraightSkeletonTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/spatialanalysis2d/StraightSkeletonTest.java new file mode 100644 index 0000000000..38956cc0af --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/spatialanalysis2d/StraightSkeletonTest.java @@ -0,0 +1,150 @@ +/* + * 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.spatialanalysis2d; + +import java.util.ArrayList; +import java.util.List; +import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.curve.LinearRing; +import org.apache.sis.maths.Maths; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; + + +/** + * Tests {@link StraightSkeleton}. + * + * @author Johann Sorel (Geomatys) + */ +public class StraightSkeletonTest { + + /** + * A rectangle gives two hips at each end joined by a ridge. + * <pre> + * 4 ■-------------------■ + * | ⟍ ⟋ | + * 2 | ■---------■ | + * | ⟋ ⟍ | + * 0 ■-------------------■ + * 0 2 8 10 + * </pre> + */ + @Test + public void testRectangle() { + assertSkeleton(skeleton(new double[] {0,0, 10,0, 10,4, 0,4, 0,0}), 6, 5, 2, 2,2, 8,2); + } + + /** + * A square gives a pyramid, all hips meeting at the center. + */ + @Test + public void testSquare() { + assertSkeleton(skeleton(new double[] {0,0, 4,0, 4,4, 0,4, 0,0}), 5, 4, 2, 2,2); + } + + /** + * A L shape, the reflex vertex creates a valley up to the junction of the two ridges. + * <pre> + * 10 ■---■ + * | | + * 4 | ■-----------■ + * | | + * 0 ■---------------■ + * 0 4 10 + * </pre> + */ + @Test + public void testLShape() { + assertSkeleton(skeleton(new double[] {0,0, 10,0, 10,4, 4,4, 4,10, 0,10, 0,0}), 9, 8, 2, 2,2, 8,2, 2,8); + } + + /** + * A square with an off-center hole, the widest part is on the right of the hole. + * <pre> + * 10 ■-------------------■ + * | | + * 6 | ■-----■ | + * | | | | + * 3 | ■-----■ | + * 0 ■-------------------■ + * 0 2 5 10 + * </pre> + */ + @Test + public void testHole() { + final StraightSkeleton skeleton = skeleton(new double[] {0,0, 10,0, 10,10, 0,10, 0,0}, + new double[] {2,3, 5,3, 5,6, 2,6, 2,3}); + assertSkeleton(skeleton, -1, -1, 2.5); + final double[] nodes = skeleton.getNodes(); + for (int i = 0; i < nodes.length; i += 3) { + final double x = nodes[i], y = nodes[i+1]; + if (nodes[i+2] > 0) { + assertTrue(x > 0 && x < 10 && y > 0 && y < 10, "Node outside the exterior ring."); + assertFalse(x >= 2 && x <= 5 && y >= 3 && y <= 6, "Node inside the hole."); + } + } + } + + /** + * Verifies the number of nodes and arcs (ignored if negative), the maximum time, the nodes at the + * maximum time given as (x,y) pairs, and that no two arcs cross each other. + */ + private static void assertSkeleton(final StraightSkeleton skeleton, final int nbNodes, final int nbArcs, + final double maxTime, final double... topNodes) + { + final double[] n = skeleton.getNodes(); + final int[] arcs = skeleton.getArcs(); + if (nbNodes >= 0) assertEquals(nbNodes, n.length / 3); + if (nbArcs >= 0) assertEquals(nbArcs, arcs.length / 2); + assertEquals(maxTime, skeleton.getMaximumTime(), 1E-9); +next: for (int k = 0; k < topNodes.length; k += 2) { + for (int i = 0; i < n.length; i += 3) { + if (Math.abs(n[i] - topNodes[k]) < 1E-9 && Math.abs(n[i+1] - topNodes[k+1]) < 1E-9) { + assertEquals(maxTime, n[i+2], 1E-9); + continue next; + } + } + fail("Missing node at " + topNodes[k] + " " + topNodes[k+1]); + } + for (int i = 0; i < arcs.length; i += 2) { + for (int j = i + 2; j < arcs.length; j += 2) { + final int a = arcs[i]*3, b = arcs[i+1]*3, c = arcs[j]*3, d = arcs[j+1]*3; + if (a == c || a == d || b == c || b == d) continue; + final boolean cross = Maths.lineSide(n[c], n[c+1], n[d], n[d+1], n[a], n[a+1]) + * Maths.lineSide(n[c], n[c+1], n[d], n[d+1], n[b], n[b+1]) < 0 + && Maths.lineSide(n[a], n[a+1], n[b], n[b+1], n[c], n[c+1]) + * Maths.lineSide(n[a], n[a+1], n[b], n[b+1], n[d], n[d+1]) < 0; + assertFalse(cross, "Crossing arcs."); + } + } + } + + private static StraightSkeleton skeleton(final double[] exterior, final double[]... holes) { + final List<LinearRing> interiors = new ArrayList<>(); + for (double[] hole : holes) { + interiors.add(ring(hole)); + } + return new StraightSkeleton(GeometryFactory.DEFAULT.createPolygon(ring(exterior), interiors)); + } + + private static LinearRing ring(final double[] coordinates) { + return GeometryFactory.DEFAULT.createLinearRing(GeometryFactory.DEFAULT.createDataPoints(Geometries.getUndefinedCRS(2), coordinates)); + } +}
