diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/SpatialIndexFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/SpatialIndexFunctions.java index 05ce51e864..9124b1eb1f 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/function/SpatialIndexFunctions.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/SpatialIndexFunctions.java @@ -35,6 +35,7 @@ import org.locationtech.jts.index.strtree.Boundable; import org.locationtech.jts.index.strtree.GeometryItemDistance; import org.locationtech.jts.index.strtree.STRtree; +//import org.locationtech.jts.triangulatepoly.VertexSequencePackedRtree; public class SpatialIndexFunctions @@ -249,4 +250,20 @@ public static Geometry monotoneChains(Geometry geom) { } return geom.getFactory().buildGeometry(lines); } + + /* + public static Geometry sprTreeBounds(Geometry geom) + { + Coordinate[] pts = geom.getCoordinates(); + VertexSequencePackedRtree index = new VertexSequencePackedRtree(pts); + Envelope[] bounds = index.getBounds(); + Geometry[] polys = new Geometry[bounds.length]; + int i = 0; + for (Envelope env : bounds) { + polys[i++] = geom.getFactory().toGeometry(env); + } + return geom.getFactory().createGeometryCollection(polys); + } + */ + } diff --git a/modules/app/src/main/java/org/locationtech/jtstest/function/TriangulatePolyFunctions.java b/modules/app/src/main/java/org/locationtech/jtstest/function/TriangulatePolyFunctions.java new file mode 100644 index 0000000000..566cf5a13d --- /dev/null +++ b/modules/app/src/main/java/org/locationtech/jtstest/function/TriangulatePolyFunctions.java @@ -0,0 +1,32 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jtstest.function; + +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.triangulate.polygon.ConstrainedDelaunayTriangulator; +import org.locationtech.jts.triangulate.polygon.PolygonTriangulator; + + +public class TriangulatePolyFunctions +{ + public static Geometry triangulate(Geometry geom) + { + return PolygonTriangulator.triangulate(geom); + } + + public static Geometry constrainedDelaunay(Geometry geom) + { + return ConstrainedDelaunayTriangulator.triangulate(geom); + } + + +} diff --git a/modules/app/src/main/java/org/locationtech/jtstest/geomfunction/GeometryFunctionRegistry.java b/modules/app/src/main/java/org/locationtech/jtstest/geomfunction/GeometryFunctionRegistry.java index 4183997a57..098cca12d4 100644 --- a/modules/app/src/main/java/org/locationtech/jtstest/geomfunction/GeometryFunctionRegistry.java +++ b/modules/app/src/main/java/org/locationtech/jtstest/geomfunction/GeometryFunctionRegistry.java @@ -65,6 +65,7 @@ import org.locationtech.jtstest.function.SpatialIndexFunctions; import org.locationtech.jtstest.function.SpatialPredicateFunctions; import org.locationtech.jtstest.function.TriangleFunctions; +import org.locationtech.jtstest.function.TriangulatePolyFunctions; import org.locationtech.jtstest.function.TriangulationFunctions; import org.locationtech.jtstest.function.ValidationFunctions; import org.locationtech.jtstest.function.WriterFunctions; @@ -131,6 +132,7 @@ public static GeometryFunctionRegistry createTestBuilderRegistry() funcRegistry.add(SnappingFunctions.class); funcRegistry.add(SortingFunctions.class); funcRegistry.add(TriangulationFunctions.class); + funcRegistry.add(TriangulatePolyFunctions.class); funcRegistry.add(TriangleFunctions.class); funcRegistry.add(ValidationFunctions.class); funcRegistry.add(WriterFunctions.class); diff --git a/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java b/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java index baa4fae5fe..8cdaa1d16f 100644 --- a/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java +++ b/modules/core/src/main/java/org/locationtech/jts/geom/Triangle.java @@ -33,12 +33,9 @@ public class Triangle * Note: this implementation is not robust for angles very close to 90 * degrees. * - * @param a - * a vertex of the triangle - * @param b - * a vertex of the triangle - * @param c - * a vertex of the triangle + * @param a a vertex of the triangle + * @param b a vertex of the triangle + * @param c a vertex of the triangle * @return true if the triangle is acute */ public static boolean isAcute(Coordinate a, Coordinate b, Coordinate c) @@ -51,6 +48,38 @@ public static boolean isAcute(Coordinate a, Coordinate b, Coordinate c) return false; return true; } + + /** + * Tests whether a triangle is oriented counter-clockwise. + * + * @param a a vertex of the triangle + * @param b a vertex of the triangle + * @param c a vertex of the triangle + * @return true if the triangle orientation is counter-clockwise + */ + public static boolean isCCW(Coordinate a, Coordinate b, Coordinate c) + { + return Orientation.COUNTERCLOCKWISE == Orientation.index(a, b, c); + } + + /** + * Tests whether a triangle intersects a point. + * + * @param a a vertex of the triangle + * @param b a vertex of the triangle + * @param c a vertex of the triangle + * @param p the point to test + * @return true if the triangle intersects the point + */ + public static boolean intersects(Coordinate a, Coordinate b, Coordinate c, Coordinate p) + { + int exteriorIndex = isCCW(a, b, c) ? + Orientation.CLOCKWISE : Orientation.COUNTERCLOCKWISE; + if (exteriorIndex == Orientation.index(a, b, p)) return false; + if (exteriorIndex == Orientation.index(b, c, p)) return false; + if (exteriorIndex == Orientation.index(c, a, p)) return false; + return true; + } /** * Computes the line which is the perpendicular bisector of the line segment @@ -437,8 +466,8 @@ public static double interpolateZ(Coordinate p, Coordinate v0, Coordinate v1, double u = (-c * dx + a * dy) / det; double z = v0.getZ() + t * (v1.getZ() - v0.getZ()) + u * (v2.getZ() - v0.getZ()); return z; - } - + } + /** * The coordinates of the vertices of the triangle */ @@ -487,9 +516,18 @@ public Coordinate inCentre() */ public boolean isAcute() { - return isAcute(this.p0, this.p1, this.p2); + return isAcute(p0, p1, p2); } + /** + * Tests whether this triangle is oriented counter-clockwise. + * + * @return true if the triangle orientation is counter-clockwise + */ + public boolean isCCW() { + return isCCW(p0, p1, p2); + } + /** * Computes the circumcentre of this triangle. The circumcentre is the centre * of the circumcircle, the smallest circle which encloses the triangle. It is @@ -507,7 +545,7 @@ public boolean isAcute() */ public Coordinate circumcentre() { - return circumcentre(this.p0, this.p1, this.p2); + return circumcentre(p0, p1, p2); } /** @@ -522,7 +560,7 @@ public Coordinate circumcentre() */ public Coordinate centroid() { - return centroid(this.p0, this.p1, this.p2); + return centroid(p0, p1, p2); } /** @@ -532,7 +570,7 @@ public Coordinate centroid() */ public double longestSideLength() { - return longestSideLength(this.p0, this.p1, this.p2); + return longestSideLength(p0, p1, p2); } /** @@ -545,7 +583,7 @@ public double longestSideLength() */ public double area() { - return area(this.p0, this.p1, this.p2); + return area(p0, p1, p2); } /** @@ -563,7 +601,7 @@ public double area() */ public double signedArea() { - return signedArea(this.p0, this.p1, this.p2); + return signedArea(p0, p1, p2); } /** @@ -574,7 +612,7 @@ public double signedArea() */ public double area3D() { - return area3D(this.p0, this.p1, this.p2); + return area3D(p0, p1, p2); } /** diff --git a/modules/core/src/main/java/org/locationtech/jts/math/MathUtil.java b/modules/core/src/main/java/org/locationtech/jts/math/MathUtil.java index 2a7f84275a..63bca1b672 100644 --- a/modules/core/src/main/java/org/locationtech/jts/math/MathUtil.java +++ b/modules/core/src/main/java/org/locationtech/jts/math/MathUtil.java @@ -47,6 +47,31 @@ public static int clamp(int x, int min, int max) return x; } + /** + * Clamps an integer to a given maximum limit. + * + * @param x the value to clamp + * @param max the maximum value + * @return the clamped value + */ + public static int clampMax(int x, int max) + { + if (x > max) return max; + return x; + } + + /** + * Computes the ceiling function of the dividend of two integers. + * + * @param num the numerator + * @param denom the denominator + * @return the ceiling of num / denom + */ + public static int ceil(int num, int denom) { + int div = num / denom; + return div * denom >= num ? div : div + 1; + } + private static final double LOG_10 = Math.log(10); /** diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/ConstrainedDelaunayTriangulator.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/ConstrainedDelaunayTriangulator.java new file mode 100644 index 0000000000..8872a2a2c1 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/ConstrainedDelaunayTriangulator.java @@ -0,0 +1,95 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import java.util.ArrayList; +import java.util.List; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.Polygon; +import org.locationtech.jts.geom.util.PolygonExtracter; +import org.locationtech.jts.triangulate.tri.Tri; +import org.locationtech.jts.triangulate.tri.TriangulationBuilder; + +/** + * Computes the Constrained Delaunay Triangulation of polygons. + * The Constrained Delaunay Triangulation of a polygon is a set of triangles + * covering the polygon, with the maximum total interior angle over all + * possible triangulations. It provides the "best quality" triangulation + * of the polygon. + *

+ * Holes are supported. + */ +public class ConstrainedDelaunayTriangulator { + + /** + * Computes the Constrained Delaunay Triangulation of each polygon element in a geometry. + * + * @param geom the input geometry + * @return a GeometryCollection of the computed triangle polygons + */ + public static Geometry triangulate(Geometry geom) { + ConstrainedDelaunayTriangulator cdt = new ConstrainedDelaunayTriangulator(geom); + return cdt.compute(); + } + + private final GeometryFactory geomFact; + private final Geometry inputGeom; + + /** + * Constructs a new Constrained Delaunay triangulator. + * + * @param inputGeom the input geometry + */ + public ConstrainedDelaunayTriangulator(Geometry inputGeom) { + geomFact = new GeometryFactory(); + this.inputGeom = inputGeom; + } + + private Geometry compute() { + List polys = PolygonExtracter.getPolygons(inputGeom); + List triList = new ArrayList(); + for (Polygon poly : polys) { + List polyTriList = triangulatePolygon(poly); + triList.addAll(polyTriList); + } + return Tri.toGeometry(triList, geomFact); + } + + /** + * Computes the triangulation of a single polygon + * and returns it as a list of {@link Tri}s. + * + * @param poly the input polygon + * @return list of Tris forming the triangulation + */ + List triangulatePolygon(Polygon poly) { + /** + * Normalize to ensure that shell and holes have canonical orientation. + * + * TODO: perhaps better to just correct orientation of rings? + */ + Polygon polyNorm = (Polygon) poly.norm(); + Coordinate[] polyShell = PolygonHoleJoiner.join(polyNorm); + List triList = PolygonEarClipper.triangulate(polyShell); + + //long start = System.currentTimeMillis(); + TriangulationBuilder.build(triList); + TriDelaunayImprover.improve(triList); + //System.out.println("swap used: " + (System.currentTimeMillis() - start) + " milliseconds"); + + return triList; + } + +} diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonEarClipper.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonEarClipper.java new file mode 100644 index 0000000000..df6dff7f79 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonEarClipper.java @@ -0,0 +1,400 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import java.util.ArrayList; +import java.util.List; + +import org.locationtech.jts.algorithm.Angle; +import org.locationtech.jts.algorithm.Orientation; +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.CoordinateList; +import org.locationtech.jts.geom.Envelope; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.Polygon; +import org.locationtech.jts.geom.Triangle; +import org.locationtech.jts.triangulate.tri.Tri; + +/** + * Triangulates a polygon using the Ear-Clipping technique. + * The polygon is provided as a closed list of contiguous vertices + * defining its boundary. + * The vertices must have clockwise orientation. + *

+ * The polygon boundary must not self-cross, + * but may self-touch at points or along an edge. + * It may contain repeated points, which are treated as a single vertex. + * By default every vertex is triangulated, + * including ones which are "flat" (the adjacent segments are collinear). + * These can be removed by setting {@link #setSkipFlatCorners(boolean)} + *

+ * The polygon representation does not allow holes. + * Polygons with holes can be triangulated by preparing them + * with {@link PolygonHoleJoiner}. + * + * @author Martin Davis + * + */ +class PolygonEarClipper { + + private static final int NO_VERTEX_INDEX = -1; + + /** + * Triangulates a polygon via ear-clipping. + * + * @param polyShell the vertices of the polygon + * @return a list of the Tris + */ + public static List triangulate(Coordinate[] polyShell) { + PolygonEarClipper clipper = new PolygonEarClipper(polyShell); + return clipper.compute(); + } + + private boolean isFlatCornersSkipped = false; + + /** + * The polygon vertices are provided in CW orientation. + * Thus for convex interior angles + * the vertices forming the angle are in CW orientation. + */ + private final Coordinate[] vertex; + + private final int[] vertexNext; + private int vertexSize; + // first available vertex index + private int vertexFirst; + + // indices for current corner + private int[] cornerIndex; + + /** + * Indexing vertices improves ear intersection testing performance a lot. + * The polyShell vertices are contiguous, so are suitable for an SPRtree. + */ + private VertexSequencePackedRtree vertexCoordIndex; + + /** + * Creates a new ear-clipper instance. + * + * @param polyShell the polygon vertices to process + */ + public PolygonEarClipper(Coordinate[] polyShell) { + this.vertex = polyShell; + + // init working storage + vertexSize = vertex.length - 1; + vertexNext = createNextLinks(vertexSize); + vertexFirst = 0; + + vertexCoordIndex = new VertexSequencePackedRtree(vertex); + } + + private static int[] createNextLinks(int size) { + int[] next = new int[size]; + for (int i = 0; i < size; i++) { + next[i] = i + 1; + } + next[size - 1] = 0; + return next; + } + + /** + * Sets whether flat corners formed by collinear adjacent line segments + * are included in the triangulation. + * Skipping flat corners reduces the number of triangles in the output. + * However, it produces a triangulation which does not include + * all input vertices. This may be undesirable for downstream processes + * (such as computing a Constrained Delaunay Triangulation for + * purposes of computing the medial axis). + *

+ * The default is to include all vertices in the result triangulation. + * This still produces a valid triangulation, with no zero-area triangles. + *

+ * Note that repeated vertices are always skipped. + * + * @param isFlatCornersSkipped whether to skip collinear vertices + */ + public void setSkipFlatCorners(boolean isFlatCornersSkipped) { + this.isFlatCornersSkipped = isFlatCornersSkipped; + } + + public List compute() { + List triList = new ArrayList(); + + /** + * Count scanned corners, to catch infinite loops + * (which indicate an algorithm bug) + */ + int cornerScanCount = 0; + + initCornerIndex(); + Coordinate[] corner = new Coordinate[3]; + fetchCorner(corner); + + /** + * Scan continuously around vertex ring, + * until all ears have been found. + */ + while (true) { + /** + * Non-convex corner- remove if flat, or skip + * (a concave corner will turn into a convex corner + * after enough ears are removed) + */ + if (! isConvex(corner)) { + // remove the corner if it is flat or a repeated point + boolean isCornerRemoved = hasRepeatedPoint(corner) + || (isFlatCornersSkipped && isFlat(corner)); + if (isCornerRemoved) { + removeCorner(); + } + cornerScanCount++; + if ( cornerScanCount > 2 * vertexSize ) { + throw new IllegalStateException("Unable to find a convex corner"); + } + } + /** + * Convex corner - check if it is a valid ear + */ + else if ( isValidEar(cornerIndex[1], corner) ) { + triList.add(Tri.create(corner)); + removeCorner(); + cornerScanCount = 0; + } + if ( cornerScanCount > 2 * vertexSize ) { + //System.out.println(toGeometry()); + throw new IllegalStateException("Unable to find a valid ear"); + } + + //--- done when all corners are processed and removed + if ( vertexSize < 3 ) { + return triList; + } + + /** + * Skip to next corner. + * This is done even after an ear is removed, + * since that creates fewer skinny triangles. + */ + nextCorner(corner); + } + } + + private boolean isValidEar(int cornerIndex, Coordinate[] corner) { + int intApexIndex = findIntersectingVertex(cornerIndex, corner); + //--- no intersections found + if (intApexIndex == NO_VERTEX_INDEX) + return true; + //--- check for duplicate corner apex vertex + if ( vertex[intApexIndex].equals2D(corner[1]) ) { + //--- a duplicate corner vertex requires a full scan + return isValidEarScan(cornerIndex, corner); + } + return false; + } + + /** + * Finds another vertex intersecting the corner triangle, if any. + * Uses the vertex spatial index for efficiency. + *

+ * Also finds any vertex which is a duplicate of the corner apex vertex, + * which then requires a full scan of the vertices to confirm ear is valid. + * This is usually a rare situation, so has little impact on performance. + * + * @param cornerIndex the index of the corner apex vertex + * @param corner the corner vertices + * @return the index of an intersecting or duplicate vertex, or {@link #NO_VERTEX_INDEX} if none + */ + private int findIntersectingVertex(int cornerIndex, Coordinate[] corner) { + Envelope cornerEnv = envelope(corner); + int[] result = vertexCoordIndex.query(cornerEnv); + + int dupApexIndex = NO_VERTEX_INDEX; + //--- check for duplicate vertices + for (int i = 0; i < result.length; i++) { + int vertIndex = result[i]; + + if (vertIndex == cornerIndex + || vertIndex == vertex.length - 1 + || isRemoved(vertIndex)) + continue; + + Coordinate v = vertex[vertIndex]; + /** + * If another vertex at the corner is found, + * need to do a full scan to check the incident segments. + * This happens when the polygon ring self-intersects, + * usually due to hold joining. + * But only report this if no properly intersecting vertex is found, + * for efficiency. + */ + if ( v.equals2D(corner[1]) ) { + dupApexIndex = vertIndex; + } + //--- don't need to check other corner vertices + else if ( v.equals2D(corner[0]) || v.equals2D(corner[2]) ) { + continue; + } + //--- this is a properly intersecting vertex + else if (Triangle.intersects(corner[0], corner[1], corner[2], v) ) + return vertIndex; + } + if (dupApexIndex != NO_VERTEX_INDEX) { + return dupApexIndex; + } + return NO_VERTEX_INDEX; + } + + /** + * Scan all vertices in current ring to check if any are duplicates + * of the corner apex vertex, and if so whether the corner ear + * intersects the adjacent segments and thus is invalid. + * + * @param cornerIndex the index of the corner apex + * @param corner the corner vertices + * @return true if the corner ia a valid ear + */ + private boolean isValidEarScan(int cornerIndex, Coordinate[] corner) { + double cornerAngle = Angle.angleBetweenOriented(corner[0], corner[1], corner[2]); + + int currIndex = nextIndex(vertexFirst); + int prevIndex = vertexFirst; + Coordinate vPrev = vertex[prevIndex]; + for (int i = 0; i < vertexSize; i++) { + Coordinate v = vertex[currIndex]; + /** + * Because of hole-joining vertices can occur more than once. + * If vertex is same as corner[1], + * check whether either adjacent edge lies inside the ear corner. + * If so the ear is invalid. + */ + if ( currIndex != cornerIndex + && v.equals2D(corner[1]) ) { + Coordinate vNext = vertex[nextIndex(currIndex)]; + + //TODO: for robustness use segment orientation instead + double aOut = Angle.angleBetweenOriented(corner[0], corner[1], vNext); + double aIn = Angle.angleBetweenOriented(corner[0], corner[1], vPrev); + if ( aOut > 0 && aOut < cornerAngle ) { + return false; + } + if ( aIn > 0 && aIn < cornerAngle ) { + return false; + } + if ( aOut == 0 && aIn == cornerAngle ) { + return false; + } + } + + //--- move to next vertex + vPrev = v; + prevIndex = currIndex; + currIndex = nextIndex(currIndex); + } + return true; + } + + private static Envelope envelope(Coordinate[] corner) { + Envelope cornerEnv = new Envelope(corner[0], corner[1]); + cornerEnv.expandToInclude(corner[2]); + return cornerEnv; + } + + /** + * Remove the corner apex vertex and update the candidate corner location. + */ + private void removeCorner() { + int cornerApexIndex = cornerIndex[1]; + if ( vertexFirst == cornerApexIndex) { + vertexFirst = vertexNext[cornerApexIndex]; + } + vertexNext[cornerIndex[0]] = vertexNext[cornerApexIndex]; + vertexCoordIndex.remove(cornerApexIndex); + vertexNext[cornerApexIndex] = NO_VERTEX_INDEX; + vertexSize--; + //-- adjust following corner indexes + cornerIndex[1] = nextIndex(cornerIndex[0]); + cornerIndex[2] = nextIndex(cornerIndex[1]); + } + + private boolean isRemoved(int vertexIndex) { + return NO_VERTEX_INDEX == vertexNext[vertexIndex]; + } + + private void initCornerIndex() { + cornerIndex = new int[3]; + cornerIndex[0] = 0; + cornerIndex[1] = 1; + cornerIndex[2] = 2; + } + + /** + * Fetch the corner vertices from the indices. + * + * @param corner an array for the corner vertices + */ + private void fetchCorner(Coordinate[] cornerVertex) { + cornerVertex[0] = vertex[cornerIndex[0]]; + cornerVertex[1] = vertex[cornerIndex[1]]; + cornerVertex[2] = vertex[cornerIndex[2]]; + } + + /** + * Move to next corner. + */ + private void nextCorner(Coordinate[] cornerVertex) { + if ( vertexSize < 3 ) { + return; + } + cornerIndex[0] = nextIndex(cornerIndex[0]); + cornerIndex[1] = nextIndex(cornerIndex[0]); + cornerIndex[2] = nextIndex(cornerIndex[1]); + fetchCorner(cornerVertex); + } + + /** + * Get the index of the next available shell coordinate starting from the given + * index. + * + * @param index candidate position + * @return index of the next available shell coordinate + */ + private int nextIndex(int index) { + return vertexNext[index]; + } + + private static boolean isConvex(Coordinate[] pts) { + return Orientation.CLOCKWISE == Orientation.index(pts[0], pts[1], pts[2]); + } + + private static boolean isFlat(Coordinate[] pts) { + return Orientation.COLLINEAR == Orientation.index(pts[0], pts[1], pts[2]); + } + + private static boolean hasRepeatedPoint(Coordinate[] pts) { + return pts[1].equals2D(pts[0]) || pts[1].equals2D(pts[2]); + } + + public Polygon toGeometry() { + GeometryFactory fact = new GeometryFactory(); + CoordinateList coordList = new CoordinateList(); + int index = vertexFirst; + for (int i = 0; i < vertexSize; i++) { + Coordinate v = vertex[index]; + index = nextIndex(index); + // if (i < shellCoordAvailable.length && shellCoordAvailable.get(i)) + coordList.add(v, true); + } + coordList.closeRing(); + return fact.createPolygon(fact.createLinearRing(coordList.toCoordinateArray())); + } +} \ No newline at end of file diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonHoleJoiner.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonHoleJoiner.java new file mode 100755 index 0000000000..5aebe3f437 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonHoleJoiner.java @@ -0,0 +1,334 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import java.util.ArrayList; +import java.util.Collections; +import java.util.Comparator; +import java.util.HashMap; +import java.util.List; +import java.util.TreeSet; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Envelope; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.LinearRing; +import org.locationtech.jts.geom.Polygon; +import org.locationtech.jts.noding.BasicSegmentString; +import org.locationtech.jts.noding.MCIndexSegmentSetMutualIntersector; +import org.locationtech.jts.noding.SegmentIntersectionDetector; +import org.locationtech.jts.noding.SegmentSetMutualIntersector; +import org.locationtech.jts.noding.SegmentString; +import org.locationtech.jts.noding.SegmentStringUtil; + +/** + * Transforms a polygon with holes into a single self-touching ring + * by connecting holes to the exterior shell or to another hole. + * The holes are added from the lowest upwards. + * As the resulting shell develops, a hole might be added to what was + * originally another hole. + */ +class PolygonHoleJoiner { + + public static Polygon joinAsPolygon(Polygon inputPolygon) { + return inputPolygon.getFactory().createPolygon(join(inputPolygon)); + } + + public static Coordinate[] join(Polygon inputPolygon) { + PolygonHoleJoiner joiner = new PolygonHoleJoiner(inputPolygon); + return joiner.compute(); + } + + private static final double EPS = 1.0E-4; + + private List shellCoords; + // orderedCoords is a copy of shellCoords for sort purposes + private TreeSet orderedCoords; + // Key: starting end of the cut; Value: list of the other end of the cut + private HashMap> cutMap; + private SegmentSetMutualIntersector polygonIntersector; + + private Polygon inputPolygon; + + public PolygonHoleJoiner(Polygon inputPolygon) { + this.inputPolygon = inputPolygon; + polygonIntersector = createPolygonIntersector(inputPolygon); + } + + /** + * Computes the joined ring. + * + * @return the points in the joined ring + */ + public Coordinate[] compute() { + //--- copy the input polygon shell coords + shellCoords = ringCoordinates(inputPolygon.getExteriorRing()); + if ( inputPolygon.getNumInteriorRing() != 0 ) { + joinHoles(); + } + return shellCoords.toArray(new Coordinate[0]); + } + + private static List ringCoordinates(LinearRing ring) { + Coordinate[] coords = ring.getCoordinates(); + List coordList = new ArrayList(); + for (Coordinate p : coords) { + coordList.add(p); + } + return coordList; + } + + private void joinHoles() { + orderedCoords = new TreeSet(); + orderedCoords.addAll(shellCoords); + cutMap = new HashMap>(); + List orderedHoles = sortHoles(inputPolygon); + for (int i = 0; i < orderedHoles.size(); i++) { + joinHole(orderedHoles.get(i)); + } + } + + /** + * Joins a single hole to the current shellRing. + * + * @param hole the hole to join + */ + private void joinHole(LinearRing hole) { + /** + * 1) Get a list of HoleVertex Index. + * 2) Get a list of ShellVertex. + * 3) Get the pair that has the shortest distance between them. + * This pair is the endpoints of the cut + * 4) The selected ShellVertex may occurs multiple times in + * shellCoords[], so find the proper one and add the hole after it. + */ + final Coordinate[] holeCoords = hole.getCoordinates(); + List holeLeftVerticesIndex = getLeftMostVertex(hole); + Coordinate holeCoord = holeCoords[holeLeftVerticesIndex.get(0)]; + List shellCoordsList = getLeftShellVertex(holeCoord); + Coordinate shellCoord = shellCoordsList.get(0); + int shortestHoleVertexIndex = 0; + //--- pick the shell-hole vertex pair that gives the shortest distance + if ( Math.abs(shellCoord.x - holeCoord.x) < EPS ) { + double shortest = Double.MAX_VALUE; + for (int i = 0; i < holeLeftVerticesIndex.size(); i++) { + for (int j = 0; j < shellCoordsList.size(); j++) { + double currLength = Math.abs(shellCoordsList.get(j).y - holeCoords[holeLeftVerticesIndex.get(i)].y); + if ( currLength < shortest ) { + shortest = currLength; + shortestHoleVertexIndex = i; + shellCoord = shellCoordsList.get(j); + } + } + } + } + int shellVertexIndex = getShellCoordIndex(shellCoord, + holeCoords[holeLeftVerticesIndex.get(shortestHoleVertexIndex)]); + addHoleToShell(shellVertexIndex, holeCoords, holeLeftVerticesIndex.get(shortestHoleVertexIndex)); + } + + /** + * Get the ith shellvertex in shellCoords[] that the current should add after + * + * @param shellVertex Coordinate of the shell vertex + * @param holeVertex Coordinate of the hole vertex + * @return the ith shellvertex + */ + private int getShellCoordIndex(Coordinate shellVertex, Coordinate holeVertex) { + int numSkip = 0; + ArrayList newValueList = new ArrayList(); + newValueList.add(holeVertex); + if ( cutMap.containsKey(shellVertex) ) { + for (Coordinate coord : cutMap.get(shellVertex)) { + if ( coord.y < holeVertex.y ) { + numSkip++; + } + } + cutMap.get(shellVertex).add(holeVertex); + } else { + cutMap.put(shellVertex, newValueList); + } + if ( !cutMap.containsKey(holeVertex) ) { + cutMap.put(holeVertex, new ArrayList(newValueList)); + } + return getShellCoordIndexSkip(shellVertex, numSkip); + } + + /** + * Find the index of the coordinate in ShellCoords ArrayList, + * skipping over some number of matches + * + * @param coord + * @return + */ + private int getShellCoordIndexSkip(Coordinate coord, int numSkip) { + for (int i = 0; i < shellCoords.size(); i++) { + if ( shellCoords.get(i).equals2D(coord, EPS) ) { + if ( numSkip == 0 ) + return i; + numSkip--; + } + } + throw new IllegalStateException("Vertex is not in shellcoords"); + } + + /** + * Gets a list of shell vertices that could be used to join with the hole. + * This list contains only one item if the chosen vertex does not share the same + * x value with holeCoord + * + * @param holeCoord the hole coordinates + * @return a list of candidate join vertices + */ + private List getLeftShellVertex(Coordinate holeCoord) { + ArrayList list = new ArrayList(); + Coordinate closest = orderedCoords.higher(holeCoord); + while (closest.x == holeCoord.x) { + closest = orderedCoords.higher(closest); + } + do { + closest = orderedCoords.lower(closest); + } while (!isJoinable(holeCoord, closest) && !closest.equals(orderedCoords.first())); + list.add(closest); + if ( closest.x != holeCoord.x ) + return list; + double chosenX = closest.x; + list.clear(); + while (chosenX == closest.x) { + list.add(closest); + closest = orderedCoords.lower(closest); + if ( closest == null ) + return list; + } + return list; + } + + /** + * Determine if a line segment between a hole vertex + * and a shell vertex lies inside the input polygon. + * + * @param holeCoord a hole coordinate + * @param shellCoord a shell coordinate + * @return true if the line lies inside the polygon + */ + private boolean isJoinable(Coordinate holeCoord, Coordinate shellCoord) { + /** + * Since the line runs between a hole and the shell, + * it is inside the polygon if it does not cross the polygon boundary. + */ + boolean isJoinable = ! crossesPolygon(holeCoord, shellCoord); + /* + //--- slow code for testing only + LineString join = geomFact.createLineString(new Coordinate[] { holeCoord, shellCoord }); + boolean isJoinableSlow = inputPolygon.covers(join) + if (isJoinableSlow != isJoinable) { + System.out.println(WKTWriter.toLineString(holeCoord, shellCoord)); + } + //Assert.isTrue(isJoinableSlow == isJoinable); + */ + return isJoinable; + } + + /** + * Tests whether a line segment crosses the polygon boundary. + * + * @param p0 a vertex + * @param p1 a vertex + * @return true if the line segment crosses the polygon boundary + */ + private boolean crossesPolygon(Coordinate p0, Coordinate p1) { + SegmentString segString = new BasicSegmentString( + new Coordinate[] { p0, p1 }, null); + List segStrings = new ArrayList(); + segStrings.add(segString); + + SegmentIntersectionDetector segInt = new SegmentIntersectionDetector(); + segInt.setFindProper(true); + polygonIntersector.process(segStrings, segInt); + + return segInt.hasProperIntersection(); + } + + /** + * Add hole at proper position in shell coordinate list. + * Also adds hole points to ordered coordinates. + * + * @param shellVertexIndex + * @param holeCoords + * @param holeVertexIndex + */ + private void addHoleToShell(int shellVertexIndex, Coordinate[] holeCoords, int holeVertexIndex) { + List newCoords = new ArrayList(); + newCoords.add(new Coordinate(shellCoords.get(shellVertexIndex))); + final int nPts = holeCoords.length - 1; + int i = holeVertexIndex; + do { + newCoords.add(new Coordinate(holeCoords[i])); + i = (i + 1) % nPts; + } while (i != holeVertexIndex); + newCoords.add(new Coordinate(holeCoords[holeVertexIndex])); + shellCoords.addAll(shellVertexIndex, newCoords); + orderedCoords.addAll(newCoords); + } + + /** + * Sort the holes by minimum X, minimum Y. + * + * @param poly polygon that contains the holes + * @return a list of ordered hole geometry + */ + private List sortHoles(final Polygon poly) { + List holes = new ArrayList(); + for (int i = 0; i < poly.getNumInteriorRing(); i++) { + holes.add(poly.getInteriorRingN(i)); + } + Collections.sort(holes, new EnvelopeComparator()); + return holes; + } + + /** + * Gets a list of indices of the leftmost vertices in a ring. + * + * @param geom the hole ring + * @return index of the left most vertex + */ + private List getLeftMostVertex(LinearRing ring) { + Coordinate[] coords = ring.getCoordinates(); + ArrayList list = new ArrayList(); + double minX = ring.getEnvelopeInternal().getMinX(); + for (int i = 0; i < coords.length; i++) { + if ( Math.abs(coords[i].x - minX) < EPS ) { + list.add(i); + } + } + return list; + } + + private SegmentSetMutualIntersector createPolygonIntersector(Polygon polygon) { + List polySegStrings = SegmentStringUtil.extractSegmentStrings(polygon); + return new MCIndexSegmentSetMutualIntersector(polySegStrings); + } + + /** + * + * @author mdavis + * + */ + private static class EnvelopeComparator implements Comparator { + public int compare(Geometry o1, Geometry o2) { + Envelope e1 = o1.getEnvelopeInternal(); + Envelope e2 = o2.getEnvelopeInternal(); + return e1.compareTo(e2); + } + } + +} diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonTriangulator.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonTriangulator.java new file mode 100755 index 0000000000..2da07e8a8e --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/PolygonTriangulator.java @@ -0,0 +1,99 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import java.util.ArrayList; +import java.util.List; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.Polygon; +import org.locationtech.jts.geom.util.PolygonExtracter; +import org.locationtech.jts.triangulate.tri.Tri; + +/** + * Computes a triangulation of each polygon in a {@link Geometry}. + * A polygon triangulation is a non-overlapping set of triangles which + * cover the polygon and have the same vertices as the polygon. + * The priority is on performance rather than triangulation quality, + * so that the output may contain many narrow triangles. + *

+ * Holes are handled by joining them to the shell to form a + * (self-touching) polygon shell with no holes. + * Although invalid, this can be triangulated effectively. + *

+ * For better-quality triangulation use {@link ConstrainedDelaunayTriangulator}. + * + * @see ConstrainedDelaunayTriangulator + * + * @author Martin Davis + * + */ +public class PolygonTriangulator { + + /** + * Computes a triangulation of each polygon in a geometry. + * + * @param geom a geometry containing polygons + * @return a GeometryCollection containing the triangle polygons + */ + public static Geometry triangulate(Geometry geom) { + PolygonTriangulator clipper = new PolygonTriangulator(geom); + return clipper.compute(); + } + + private final GeometryFactory geomFact; + private final Geometry inputGeom; + + /** + * Constructs a new triangulator. + * + * @param inputGeom the input geometry + */ + public PolygonTriangulator(Geometry inputGeom) { + geomFact = new GeometryFactory(); + this.inputGeom = inputGeom; + } + + private Geometry compute() { + @SuppressWarnings("unchecked") + List polys = PolygonExtracter.getPolygons(inputGeom); + List triList = new ArrayList(); + for (Polygon poly : polys) { + List polyTriList = triangulatePolygon(poly); + triList.addAll(polyTriList); + } + return Tri.toGeometry(triList, geomFact); + } + + /** + * Computes the triangulation of a single polygon + * + * @return GeometryCollection of triangular polygons + */ + private List triangulatePolygon(Polygon poly) { + /** + * Normalize to ensure that shell and holes have canonical orientation. + * + * TODO: perhaps better to just correct orientation of rings? + */ + Polygon polyNorm = (Polygon) poly.norm(); + Coordinate[] polyShell = PolygonHoleJoiner.join(polyNorm); + + List triList = PolygonEarClipper.triangulate(polyShell); + //Tri.validate(triList); + + return triList; + } + +} diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/TriDelaunayImprover.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/TriDelaunayImprover.java new file mode 100755 index 0000000000..c57d03e278 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/TriDelaunayImprover.java @@ -0,0 +1,186 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import java.util.List; + +import org.locationtech.jts.algorithm.Orientation; +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.triangulate.quadedge.TrianglePredicate; +import org.locationtech.jts.triangulate.tri.Tri; +import org.locationtech.jts.triangulate.tri.TriangulationBuilder; + +/** + * Improves the quality of a triangulation of {@link Tri}s via + * iterated Delaunay flipping. + * This produces a Constrained Delaunay Triangulation + * with the constraints being the boundary of the input triangulation. + * + * @author mdavis + */ +class TriDelaunayImprover { + + /** + * Improves the quality of a triangulation of {@link Tri}s via + * iterated Delaunay flipping. + * The Tris are assumed to be linked into a Triangulation + * (e.g. via {@link TriangulationBuilder}). + * + * @param triList the list of Tris to flip. + */ + public static void improve(List triList) { + TriDelaunayImprover improver = new TriDelaunayImprover(triList); + improver.improve(); + } + + private static int MAX_ITERATION = 200; + private List triList; + + private TriDelaunayImprover(List triList) { + this.triList = triList; + } + + private void improve() { + for (int i = 0; i < MAX_ITERATION; i++) { + int improveCount = improveScan(triList); + //System.out.println("improve #" + i + " - count = " + improveCount); + if ( improveCount == 0 ) { + return; + } + } + } + + /** + * Improves a triangulation by examining pairs of adjacent triangles + * (forming a quadrilateral) and testing if flipping the diagonal of + * the quadrilateral would produce two new triangles with larger minimum + * interior angles. + * + * @return the number of flips that were made + */ + private int improveScan(List triList) { + int improveCount = 0; + for (int i = 0; i < triList.size() - 1; i++) { + Tri tri = triList.get(i); + for (int j = 0; j < 3; j++) { + //Tri neighb = tri.getAdjacent(j); + //tri.validateAdjacent(j); + if ( improveNonDelaunay(tri, j) ) { + // TODO: improve performance by only rescanning tris adjacent to flips? + improveCount++; + } + } + } + return improveCount; + } + + /** + * Does a flip of the common edge of two Tris if the Delaunay condition is not met. + * + * @param tri0 a Tri + * @param tri1 a Tri + * @return true if the triangles were flipped + */ + private boolean improveNonDelaunay(Tri tri, int index) { + if ( tri == null ) { + return false; + } + Tri tri1 = tri.getAdjacent(index); + if ( tri1 == null ) { + return false; + } + //tri0.validate(); + //tri1.validate(); + + + int index1 = tri1.getIndex(tri); + + Coordinate adj0 = tri.getCoordinate(index); + Coordinate adj1 = tri.getCoordinate(Tri.next(index)); + Coordinate opp0 = tri.getCoordinate(Tri.oppVertex(index)); + Coordinate opp1 = tri1.getCoordinate(Tri.oppVertex(index1)); + + /** + * The candidate new edge is opp0 - opp1. + * Check if it is inside the quadrilateral formed by the two triangles. + * This is the case if the quadrilateral is convex. + */ + if ( ! isConvex(adj0, adj1, opp0, opp1) ) { + return false; + } + + /** + * The candidate edge is inside the quadrilateral. Check to see if the flipping + * criteria is met. The flipping criteria is to flip if the two triangles are + * not Delaunay (i.e. one of the opposite vertices is in the circumcircle of the + * other triangle). + */ + if ( ! isDelaunay(adj0, adj1, opp0, opp1) ) { + tri.flip(index); + return true; + } + return false; + } + + /** + * Tests if the quadrilateral formed by two adjacent triangles is convex. + * opp0-adj0-adj1 and opp1-adj1-adj0 are the triangle corners + * and hence are known to be convex. + * The quadrilateral is convex if the other corners opp0-adj0-opp1 + * and opp1-adj1-opp0 have the same orientation (since at least one must be convex). + * + * @param adj0 adjacent edge vertex 0 + * @param adj1 adjacent edge vertex 1 + * @param opp0 corner vertex of triangle 0 + * @param opp1 corner vertex of triangle 1 + * @return true if the quadrilateral is convex + */ + private static boolean isConvex(Coordinate adj0, Coordinate adj1, Coordinate opp0, Coordinate opp1) { + int dir0 = Orientation.index(opp0, adj0, opp1); + int dir1 = Orientation.index(opp1, adj1, opp0); + boolean isConvex = dir0 == dir1; + return isConvex; + } + + /** + * Tests if either of a pair of adjacent triangles satisfy the Delaunay condition. + * The triangles are opp0-adj0-adj1 and opp1-adj1-adj0. + * The Delaunay condition is not met if one opposite vertex + * lies is in the circumcircle of the other triangle. + * + * @param adj0 adjacent edge vertex 0 + * @param adj1 adjacent edge vertex 1 + * @param opp0 corner vertex of triangle 0 + * @param opp1 corner vertex of triangle 1 + * @return true if the triangles are Delaunay + */ + private static boolean isDelaunay(Coordinate adj0, Coordinate adj1, Coordinate opp0, Coordinate opp1) { + if (isInCircle(adj0, adj1, opp0, opp1)) return false; + if (isInCircle(adj1, adj0, opp1, opp0)) return false; + return true; + } + + /** + * Tests whether a point p is in the circumcircle of a triangle abc + * (oriented clockwise). + * @param a a vertex of the triangle + * @param b a vertex of the triangle + * @param c a vertex of the triangle + * @param p the point + * + * @return true if the point is in the circumcircle + */ + private static boolean isInCircle(Coordinate a, Coordinate b, Coordinate c, Coordinate p) { + return TrianglePredicate.isInCircleRobust(a, c, b, p); + } + +} \ No newline at end of file diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/VertexSequencePackedRtree.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/VertexSequencePackedRtree.java new file mode 100644 index 0000000000..f1fd182925 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/VertexSequencePackedRtree.java @@ -0,0 +1,256 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Envelope; +import org.locationtech.jts.math.MathUtil; +import org.locationtech.jts.util.IntArrayList; + +/** + * A semi-static spatial index for points which occur + * in a spatially-coherent sequence. + * In particular, this is suitable for indexing the vertices + * of a {@link LineString} or {@link Polygon} ring. + *

+ * The index is constructed in a batch fashion on a given sequence of coordinates. + * Coordinates can be removed via the {@link #remove(int)} method. + *

+ * Note that this index queries only the individual points + * of the input coordinate sequence, + * not any line segments which might be lie between them. + * + * @author Martin Davis + * + */ +class VertexSequencePackedRtree { + + /** + * Number of items/nodes in a parent node. + * Determined empirically. Performance is not too sensitive to this. + */ + private static final int NODE_CAPACITY = 16; + + private Coordinate[] items; + private int[] levelOffset; + private int nodeCapacity = NODE_CAPACITY; + private Envelope[] bounds; + + /** + * Creates a new tree over the given sequence of coordinates. + * The sequence should be spatially coherent to provide query performance. + * + * @param pts a sequence of points + */ + public VertexSequencePackedRtree(Coordinate[] pts) { + this.items = pts; + build(); + } + + public Envelope[] getBounds() { + return bounds.clone(); + } + + private void build() { + levelOffset = computeLevelOffsets(); + bounds = createBounds(); + } + + /** + * Computes the level offsets. + * This is the position in the bounds array of each level. + * + * The levelOffsets array includes a sentinel value of offset[0] = 0. + * The top level is always of size 1, + * and so also indicates the total number of bounds. + * + * @return the level offsets + */ + private int[] computeLevelOffsets() { + IntArrayList offsets = new IntArrayList(); + offsets.add(0); + int levelSize = items.length; + int currOffset = 0; + do { + levelSize = levelNodeCount(levelSize); + currOffset += levelSize; + offsets.add(currOffset); + } while (levelSize > 1); + return offsets.toArray(); + } + + private int levelNodeCount(int numNodes) { + return MathUtil.ceil(numNodes, nodeCapacity); + } + + private Envelope[] createBounds() { + int boundsSize = levelOffset[levelOffset.length - 1] + 1; + Envelope[] bounds = new Envelope[boundsSize]; + fillItemBounds(bounds); + + for (int lvl = 1; lvl < levelOffset.length; lvl++) { + fillLevelBounds(lvl, bounds); + } + return bounds; + } + + private void fillLevelBounds(int lvl, Envelope[] bounds) { + int levelStart = levelOffset[lvl - 1]; + int levelEnd = levelOffset[lvl]; + int nodeStart = levelStart; + int levelBoundIndex = levelOffset[lvl]; + do { + int nodeEnd = MathUtil.clampMax(nodeStart + nodeCapacity, levelEnd); + bounds[levelBoundIndex++] = computeNodeEnvelope(bounds, nodeStart, nodeEnd); + nodeStart = nodeEnd; + } + while (nodeStart < levelEnd); + } + + private void fillItemBounds(Envelope[] bounds) { + int nodeStart = 0; + int boundIndex = 0; + do { + int nodeEnd = MathUtil.clampMax(nodeStart + nodeCapacity, items.length); + bounds[boundIndex++] = computeItemEnvelope(items, nodeStart, nodeEnd); + nodeStart = nodeEnd; + } + while (nodeStart < items.length); + } + + private static Envelope computeNodeEnvelope(Envelope[] bounds, int start, int end) { + Envelope env = new Envelope(); + for (int i = start; i < end; i++) { + env.expandToInclude(bounds[i]); + } + return env; + } + + private static Envelope computeItemEnvelope(Coordinate[] items, int start, int end) { + Envelope env = new Envelope(); + for (int i = start; i < end; i++) { + env.expandToInclude(items[i]); + } + return env; + } + + //------------------------ + + /** + * Queries the index to find all items which intersect an extent. + * The query result is a list of the indices of input coordinates + * which intersect the extent. + * + * @param queryEnv the query extent + * @return an array of the indices of the input coordinates + */ + public int[] query(Envelope queryEnv) { + IntArrayList resultList = new IntArrayList(); + int level = levelOffset.length - 1; + queryNode(queryEnv, level, 0, resultList); + int[] result = resultList.toArray(); + return result; + } + + private void queryNode(Envelope queryEnv, int level, int nodeIndex, IntArrayList resultList) { + int boundsIndex = levelOffset[level] + nodeIndex; + Envelope nodeEnv = bounds[boundsIndex]; + //--- node is empty + if (nodeEnv == null) + return; + if (! queryEnv.intersects(nodeEnv)) + return; + + int childNodeIndex = nodeIndex * nodeCapacity; + if (level == 0) { + queryItemRange(queryEnv, childNodeIndex, resultList); + } + else { + queryNodeRange(queryEnv, level - 1, childNodeIndex, resultList); + } + } + + private void queryNodeRange(Envelope queryEnv, int level, int nodeStartIndex, IntArrayList resultList) { + int levelMax = levelSize(level); + for (int i = 0; i < nodeCapacity; i++) { + int index = nodeStartIndex + i; + if (index >= levelMax) + return; + queryNode(queryEnv, level, index, resultList); + } + } + + private int levelSize(int level) { + return levelOffset[level + 1] - levelOffset[level]; + } + + private void queryItemRange(Envelope queryEnv, int itemIndex, IntArrayList resultList) { + for (int i = 0; i < nodeCapacity; i++) { + int index = itemIndex + i; + if (index >= items.length) + return; + Coordinate p = items[index]; + if (p != null + && queryEnv.contains(p)) + resultList.add(index); + } + } + + //------------------------ + + /** + * Removes the input item at the given index from the spatial index. + * + * @param index the index of the item in the input + */ + public void remove(int index) { + items[index] = null; + + //--- prune the item parent node if all its items are removed + int nodeIndex = index / nodeCapacity; + if (! isItemsNodeEmpty(nodeIndex)) + return; + + bounds[nodeIndex] = null; + + if (levelOffset.length <= 2) + return; + + //-- prune the node parent if all children removed + int nodeLevelIndex = nodeIndex / nodeCapacity; + if (! isNodeEmpty(1, nodeLevelIndex)) + return; + int nodeIndex1 = levelOffset[1] + nodeLevelIndex; + bounds[nodeIndex1] = null; + + //TODO: propagate removal up the tree nodes? + } + + private boolean isNodeEmpty(int level, int index) { + int start = index * nodeCapacity; + int end = MathUtil.clampMax(start + nodeCapacity, levelOffset[level]); + for (int i = start; i < end; i++) { + if (bounds[i] != null) return false; + } + return true; + } + + private boolean isItemsNodeEmpty(int nodeIndex) { + int start = nodeIndex * nodeCapacity; + int end = MathUtil.clampMax(start + nodeCapacity, items.length); + for (int i = start; i < end; i++) { + if (items[i] != null) return false; + } + return true; + } + +} diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/package-info.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/package-info.java new file mode 100644 index 0000000000..43d687cc57 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/polygon/package-info.java @@ -0,0 +1,19 @@ +/* + * Copyright (c) 2021 Martin Davis + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ + +/** + * Classes for triangulating polygons. + * {@link ConstrainedDelaunayTriangulator} can be used to provide high-quality + * near-Delaunay triangulations of polygonal geometry. + * The {@link PolygonTriangulator} produces lower-quality but faster triangulations. + */ +package org.locationtech.jts.triangulate.polygon; \ No newline at end of file diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/Tri.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/Tri.java new file mode 100755 index 0000000000..337905e815 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/Tri.java @@ -0,0 +1,549 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.tri; + +import java.util.List; + +import org.locationtech.jts.algorithm.Orientation; +import org.locationtech.jts.algorithm.RobustLineIntersector; +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.geom.GeometryFactory; +import org.locationtech.jts.geom.Polygon; +import org.locationtech.jts.io.WKTWriter; +import org.locationtech.jts.util.Assert; + +/** + * A memory-efficient representation of a triangle in a triangulation. + * Contains three vertices, and links to adjacent Tris for each edge. + * Tris are constructed independently, and if needed linked + * into a triangulation using {@link TriangulationBuilder}. + * + * @author mdavis + * + */ +public class Tri { + + /** + * Creates a {@link GeometryCollection} of {@link Polygon}s + * representing the triangles in a list. + * + * @param triList a list of Tris + * @param geomFact the GeometryFactory to use + * @return the polygons for the triangles + */ + public static Geometry toGeometry(List triList, GeometryFactory geomFact) { + Geometry[] geoms = new Geometry[triList.size()]; + for (int i = 0; i < triList.size(); i++) { + geoms[i] = triList.get(i).toPolygon(geomFact); + } + return geomFact.createGeometryCollection(geoms); + } + + /** + * Validates a list of Tris. + * + * @param triList the tris to validate + */ + public static void validate(List triList) { + for (Tri tri : triList) { + tri.validate(); + } + } + + /** + * Creates a triangle with the given vertices. + * The vertices should be oriented clockwise. + * + * @param p0 the first triangle vertex + * @param p1 the second triangle vertex + * @param p2 the third triangle vertex + * @return the created triangle + */ + public static Tri create(Coordinate p0, Coordinate p1, Coordinate p2) { + return new Tri(p0, p1, p2); + } + + /** + * Creates a triangle from an array with three vertex coordinates. + * The vertices should be oriented clockwise. + * + * @param pts the array of vertex coordinates + * @return the created triangle + */ + public static Tri create(Coordinate[] pts) { + return new Tri(pts[0], pts[1], pts[2]); + } + + private Coordinate p0; + private Coordinate p1; + private Coordinate p2; + + /** + * triN is the adjacent triangle across the edge pN - pNN. + * pNN is the next vertex CW from pN. + */ + private Tri tri0; + private Tri tri1; + private Tri tri2; + + /** + * Creates a triangle with the given vertices. + * The vertices should be oriented clockwise. + * + * @param p0 the first triangle vertex + * @param p1 the second triangle vertex + * @param p2 the third triangle vertex + */ + public Tri(Coordinate p0, Coordinate p1, Coordinate p2) { + this.p0 = p0; + this.p1 = p1; + this.p2 = p2; + //Assert.isTrue( Orientation.CLOCKWISE != Orientation.index(p0, p1, p2), "Tri is not oriented correctly"); + } + + /** + * Sets the adjacent triangles. + * The vertices of the adjacent triangles are + * assumed to match the appropriate vertices in this triangle. + * + * @param tri0 the triangle adjacent to edge 0 + * @param tri1 the triangle adjacent to edge 1 + * @param tri2 the triangle adjacent to edge 2 + */ + public void setAdjacent(Tri tri0, Tri tri1, Tri tri2) { + this.tri0 = tri0; + this.tri1 = tri1; + this.tri2 = tri2; + } + + /** + * Sets the triangle adjacent to the edge originating + * at a given vertex. + * The vertices of the adjacent triangles are + * assumed to match the appropriate vertices in this triangle. + * + * @param pt the edge start point + * @param tri the adjacent triangle + */ + public void setAdjacent(Coordinate pt, Tri tri) { + int index = getIndex(pt); + setTri(index, tri); + // TODO: validate that tri is adjacent at the edge specified + } + + /** + * Sets the triangle adjacent to an edge. + * The vertices of the adjacent triangle are + * assumed to match the appropriate vertices in this triangle. + * + * @param edgeIndex the edge triangle is adjacent to + * @param tri the adjacent triangle + */ + public void setTri(int edgeIndex, Tri tri) { + switch (edgeIndex) { + case 0: tri0 = tri; return; + case 1: tri1 = tri; return; + case 2: tri2 = tri; return; + } + Assert.shouldNeverReachHere(); + } + + private void setCoordinates(Coordinate p0, Coordinate p1, Coordinate p2) { + this.p0 = p0; + this.p1 = p1; + this.p2 = p2; + //Assert.isTrue( Orientation.CLOCKWISE != Orientation.index(p0, p1, p2), "Tri is not oriented correctly"); + } + + /** + * Spits a triangle by a point located inside the triangle. + * Creates the three new resulting triangles with adjacent links + * set correctly. + * Returns the new triangle whose 0'th vertex is the splitting point. + * + * @param p the point to insert + * @return the new triangle whose 0'th vertex is p + */ + public Tri split(Coordinate p) { + Tri tt0 = new Tri(p, p0, p1); + Tri tt1 = new Tri(p, p1, p2); + Tri tt2 = new Tri(p, p2, p0); + tt0.setAdjacent(tt2, tri0, tt1); + tt1.setAdjacent(tt0, tri1, tt2); + tt2.setAdjacent(tt1, tri2, tt0); + return tt0; + } + + /** + * Interchanges the vertices of this triangle and a neighbor + * so that their common edge + * becomes the the other diagonal of the quadrilateral they form. + * Neighbour triangle links are modified accordingly. + * + * @param index the index of the adjacent tri to flip with + */ + public void flip(int index) { + Tri tri = getAdjacent(index); + int index1 = tri.getIndex(this); + + Coordinate adj0 = getCoordinate(index); + Coordinate adj1 = getCoordinate(next(index)); + Coordinate opp0 = getCoordinate(oppVertex(index)); + Coordinate opp1 = tri.getCoordinate(oppVertex(index1)); + + flip(tri, index, index1, adj0, adj1, opp0, opp1); + } + + private void flip(Tri tri, int index0, int index1, Coordinate adj0, Coordinate adj1, Coordinate opp0, Coordinate opp1) { + //System.out.println("Flipping: " + this + " -> " + tri); + + //validate(); + //tri.validate(); + + this.setCoordinates(opp1, opp0, adj0); + tri.setCoordinates(opp0, opp1, adj1); + /** + * Order: 0: opp0-adj0 edge, 1: opp0-adj1 edge, + * 2: opp1-adj0 edge, 3: opp1-adj1 edge + */ + Tri[] adjacent = getAdjacentTris(tri, index0, index1); + this.setAdjacent(tri, adjacent[0], adjacent[2]); + //--- update the adjacent triangles with new adjacency + if ( adjacent[2] != null ) { + adjacent[2].replace(tri, this); + } + tri.setAdjacent(this, adjacent[3], adjacent[1]); + if ( adjacent[1] != null ) { + adjacent[1].replace(this, tri); + } + //validate(); + //tri.validate(); + } + + /** + * Replaces an adjacent triangle with a different one. + * + * @param triOld an adjacent triangle + * @param triNew the triangle to replace it with + */ + private void replace(Tri triOld, Tri triNew) { + if ( tri0 != null && tri0 == triOld ) { + tri0 = triNew; + } else if ( tri1 != null && tri1 == triOld ) { + tri1 = triNew; + } else if ( tri2 != null && tri2 == triOld ) { + tri2 = triNew; + } + } + + /** + * Gets the triangles adjacent to the quadrilateral + * formed by this triangle and an adjacent one. + * The triangles are returned in the following order: + *

+ * Order: 0: opp0-adj0 edge, 1: opp0-adj1 edge, + * 2: opp1-adj0 edge, 3: opp1-adj1 edge + * + * @param tri1 an adjacent triangle + * @param index the index of the common edge in this triangle + * @param index1 the index of the common edge in the adjacent triangle + * @return + */ + private Tri[] getAdjacentTris(Tri triAdj, int index, int indexAdj) { + Tri[] adj = new Tri[4]; + adj[0] = getAdjacent(prev(index)); + adj[1] = getAdjacent(next(index)); + adj[2] = triAdj.getAdjacent(next(indexAdj)); + adj[3] = triAdj.getAdjacent(prev(indexAdj)); + return adj; + } + + /** + * Validates that a tri is correct. + * Currently just checks that orientation is CW. + * + * @throw IllegalArgumentException if tri is not valid + */ + public void validate() { + if ( Orientation.CLOCKWISE != Orientation.index(p0, p1, p2) ) { + throw new IllegalArgumentException("Tri is not oriented correctly"); + } + + validateAdjacent(0); + validateAdjacent(1); + validateAdjacent(2); + } + + /** + * Validates that the vertices of an adjacent linked triangle are correct. + * + * @param index the index of the adjacent triangle + */ + public void validateAdjacent(int index) { + Tri tri = getAdjacent(index); + if (tri == null) return; + + assert(this.isAdjacent(tri)); + assert(tri.isAdjacent(this)); + + Coordinate e0 = getCoordinate(index); + Coordinate e1 = getCoordinate(next(index)); + int indexNeighbor = tri.getIndex(this); + Coordinate n0 = tri.getCoordinate(indexNeighbor); + Coordinate n1 = tri.getCoordinate(next(indexNeighbor)); + Assert.isTrue(e0.equals2D(n1), "Edge coord not equal"); + Assert.isTrue(e1.equals2D(n0), "Edge coord not equal"); + + //--- check that no edges cross + RobustLineIntersector li = new RobustLineIntersector(); + for (int i = 0; i < 3; i++) { + for (int j = 0; j < 3; j++) { + Coordinate p00 = getCoordinate(i); + Coordinate p01 = getCoordinate(next(i)); + Coordinate p10 = tri.getCoordinate(j); + Coordinate p11 = tri.getCoordinate(next(j)); + li.computeIntersection(p00, p01, p10, p11); + assert(! li.isProper()); + } + } + } + + /** + * Gets the start and end vertex of the edge adjacent to another triangle. + * + * @param neighbor + * @return + */ + /* + //TODO: define when needed + public Coordinate[] getEdge(Tri neighbor) { + int index = getIndex(neighbor); + int next = next(index); + + Coordinate e0 = getCoordinate(index); + Coordinate e1 = getCoordinate(next); + assert (neighbor.hasCoordinate(e0)); + assert (neighbor.hasCoordinate(e1)); + int iN = neighbor.getIndex(e0); + int iNPrev = prev(iN); + assert (neighbor.getIndex(e1) == iNPrev); + + return new Coordinate[] { getCoordinate(index), getCoordinate(next) }; + } + + public Coordinate getEdgeStart(int i) { + return getCoordinate(i); + } + + public Coordinate getEdgeEnd(int i) { + return getCoordinate(next(i)); + } + + public boolean hasCoordinate(Coordinate v) { + if ( p0.equals(v) || p1.equals(v) || p2.equals(v) ) { + return true; + } + return false; + } + */ + + /** + * Gets the coordinate for a vertex. + * This is the start vertex of the edge. + * + * @param index the vertex (edge) index + * @return the vertex coordinate + */ + public Coordinate getCoordinate(int index) { + if ( index == 0 ) { + return p0; + } + if ( index == 1 ) { + return p1; + } + return p2; + } + + /** + * Gets the index of the triangle vertex which has a given coordinate (if any). + * This is also the index of the edge which originates at the vertex. + * + * @param p the coordinate to find + * @return the vertex index, or -1 if it is not in the triangle + */ + public int getIndex(Coordinate p) { + if ( p0.equals2D(p) ) + return 0; + if ( p1.equals2D(p) ) + return 1; + if ( p2.equals2D(p) ) + return 2; + return -1; + } + + /** + * Gets the edge index which a triangle is adjacent to (if any), + * based on the adjacent triangle link. + * + * @param tri the tri to find + * @return the index of the edge adjacent to the triangle, or -1 if not found + */ + public int getIndex(Tri tri) { + if ( tri0 == tri ) + return 0; + if ( tri1 == tri ) + return 1; + if ( tri2 == tri ) + return 2; + return -1; + } + + /** + * Gets the triangle adjacent to an edge. + * + * @param index the edge index + * @return the adjacent triangle (may be null) + */ + public Tri getAdjacent(int index) { + switch(index) { + case 0: return tri0; + case 1: return tri1; + case 2: return tri2; + } + Assert.shouldNeverReachHere(); + return null; + } + + /** + * Tests if there is an adjacent triangle to an edge. + * + * @param index the edge index + * @return true if there is a triangle adjacent to edge + */ + public boolean hasAdjacent(int index) { + return null != getAdjacent(index); + } + + /** + * Tests if a triangle is adjacent to some edge of this triangle. + * + * @param tri the triangle to test + * @return true if the triangle is adjacent + * @see getIndex(Tri) + */ + public boolean isAdjacent(Tri tri) { + return getIndex(tri) >= 0; + } + + /** + * Computes the number of triangle adjacent to this triangle. + * This is a number in the range [0,2]. + * + * @return the number of adjacent triangles + */ + public int numAdjacent() { + int num = 0; + if ( tri0 != null ) + num++; + if ( tri1 != null ) + num++; + if ( tri2 != null ) + num++; + return num; + } + + /** + * Computes the vertex or edge index which is the next one + * (clockwise) around the triangle. + * + * @param index the index + * @return the next index value + */ + public static int next(int index) { + switch (index) { + case 0: return 1; + case 1: return 2; + case 2: return 0; + } + return -1; + } + + /** + * Computes the vertex or edge index which is the previous one + * (counter-clockwise) around the triangle. + * + * @param index the index + * @return the previous index value + */ + public static int prev(int index) { + switch (index) { + case 0: return 2; + case 1: return 0; + case 2: return 1; + } + return -1; + } + + /** + * Gets the index of the vertex opposite an edge. + * + * @param edgeIndex the edge index + * @return the index of the opposite vertex + */ + public static int oppVertex(int edgeIndex) { + return prev(edgeIndex); + } + + /** + * Gets the index of the edge opposite a vertex. + * + * @param vertexIndex the index of the vertex + * @return the index of the opposite edge + */ + public static int oppEdge(int vertexIndex) { + return next(vertexIndex); + } + + /** + * Computes a coordinate for the midpoint of a triangle edge. + * + * @param edgeIndex the edge index + * @return the midpoint of the triangle edge + */ + public Coordinate midpoint(int edgeIndex) { + Coordinate p0 = getCoordinate(edgeIndex); + Coordinate p1 = getCoordinate(next(edgeIndex)); + double midX = (p0.getX() + p1.getX()) / 2; + double midY = (p0.getY() + p1.getY()) / 2; + return new Coordinate(midX, midY); + } + + /** + * Creates a {@link Polygon} representing this triangle. + * + * @param geomFact the geometry factory + * @return a polygon + */ + public Polygon toPolygon(GeometryFactory geomFact) { + return geomFact.createPolygon( + geomFact.createLinearRing(new Coordinate[] { p0.copy(), p1.copy(), p2.copy(), p0.copy() }), null); + } + + @Override + public String toString() { + return String.format("POLYGON ((%s, %s, %s, %s))", + WKTWriter.format(p0), WKTWriter.format(p1), WKTWriter.format(p2), + WKTWriter.format(p0)); + } + +} diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriEdge.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriEdge.java new file mode 100755 index 0000000000..e7fe3a42b0 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriEdge.java @@ -0,0 +1,68 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.tri; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.io.WKTWriter; + +/** + * Represents an edge in a {@link Tri}, + * to be used as a key for looking up Tris + * while building a triangulation. + * The edge value is normalized to allow lookup + * of adjacent triangles. + * + * @author mdavis + * + */ +class TriEdge { + public Coordinate p0; + public Coordinate p1; + + public TriEdge(Coordinate a, Coordinate b) { + p0 = a; + p1 = b; + normalize(); + } + + private void normalize() { + if ( p0.compareTo(p1) < 0 ) { + Coordinate tmp = p0; + p0 = p1; + p1 = tmp; + } + } + + @Override + public int hashCode() { + int result = 17; + result = 37 * result + Coordinate.hashCode(p0.x); + result = 37 * result + Coordinate.hashCode(p1.x); + result = 37 * result + Coordinate.hashCode(p0.y); + result = 37 * result + Coordinate.hashCode(p1.y); + return result; + } + + @Override + public boolean equals(Object arg) { + if ( !(arg instanceof TriEdge) ) + return false; + TriEdge other = (TriEdge) arg; + if ( p0.equals(other.p0) && p1.equals(other.p1) ) + return true; + return false; + } + + public String toString() { + return WKTWriter.toLineString(new Coordinate[] { p0, p1}); + } +} diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriangulationBuilder.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriangulationBuilder.java new file mode 100755 index 0000000000..a4828fc9d4 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/TriangulationBuilder.java @@ -0,0 +1,82 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.tri; + +import java.util.HashMap; +import java.util.List; + +import org.locationtech.jts.geom.Coordinate; + +/** + * Builds a triangulation from a set of {@link Tri}s + * by populating the links to adjacent triangles. + * + * @author mdavis + * + */ +public class TriangulationBuilder { + + /** + * Computes the triangulation of a set of {@link Tri}s. + * + * @param triList the list of Tris + */ + public static void build(List triList) { + new TriangulationBuilder(triList); + } + + private HashMap triMap; + + /** + * Computes the triangulation of a set of {@link Tri}s. + * + * @param triList the list of Tris + */ + private TriangulationBuilder(List triList) { + triMap = new HashMap(); + for (Tri tri : triList) { + add(tri); + } + } + + private Tri find(Coordinate p0, Coordinate p1) { + TriEdge e = new TriEdge(p0, p1); + return triMap.get(e); + } + + private void add(Tri tri) { + Coordinate p0 = tri.getCoordinate(0); + Coordinate p1 = tri.getCoordinate(1); + Coordinate p2 = tri.getCoordinate(2); + + // get adjacent triangles, if any + Tri n0 = find(p0, p1); + Tri n1 = find(p1, p2); + Tri n2 = find(p2, p0); + + tri.setAdjacent(n0, n1, n2); + addAdjacent(tri, n0, p0, p1); + addAdjacent(tri, n1, p1, p2); + addAdjacent(tri, n2, p2, p0); + } + + private void addAdjacent(Tri tri, Tri adj, Coordinate p0, Coordinate p1) { + /** + * If adjacent is null, this tri is first one to be recorded for edge + */ + if (adj == null) { + triMap.put(new TriEdge(p0, p1), tri); + return; + } + adj.setAdjacent(p1, tri); + } +} diff --git a/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/package-info.java b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/package-info.java new file mode 100644 index 0000000000..e6b47f6730 --- /dev/null +++ b/modules/core/src/main/java/org/locationtech/jts/triangulate/tri/package-info.java @@ -0,0 +1,18 @@ +/* + * Copyright (c) 2021 Martin Davis + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ + +/** + * Classes for representing a planar triangulation as a set of linked triangles. + * Triangles are represented by memory-efficient {@link Tri} objects. + * A set of triangles can be linked into a triangulation using {@link TriangulationBuilder}. + */ +package org.locationtech.jts.triangulate.tri; \ No newline at end of file diff --git a/modules/core/src/test/java/org/locationtech/jts/geom/TriangleTest.java b/modules/core/src/test/java/org/locationtech/jts/geom/TriangleTest.java index aa10d55340..938f46d516 100644 --- a/modules/core/src/test/java/org/locationtech/jts/geom/TriangleTest.java +++ b/modules/core/src/test/java/org/locationtech/jts/geom/TriangleTest.java @@ -15,12 +15,13 @@ import junit.framework.TestCase; import junit.textui.TestRunner; +import test.jts.GeometryTestCase; /** * @version 1.7 */ -public class TriangleTest extends TestCase +public class TriangleTest extends GeometryTestCase { private PrecisionModel precisionModel = new PrecisionModel(); @@ -155,7 +156,7 @@ public void testCentroid() throws Exception checkCentroid("POLYGON((10 10, 20 10, 15 20, 10 10))", new Coordinate( (10.0 + 20.0 + 15.0) / 3.0, (10.0 + 10.0 + 20.0) / 3.0)); } - + public void checkCentroid(String wkt, Coordinate expectedValue) throws Exception { @@ -222,5 +223,40 @@ public void checkLongestSideLength(String wkt, double expectedValue) //System.out.println("(Instance) longestSideLength = " + length); assertEquals(expectedValue, length, 0.00000001); } + + //=============================================================== + + public void testIsCCW() { + checkIsCCW("POLYGON ((30 90, 80 50, 20 20, 30 90))", false); + checkIsCCW("POLYGON ((90 90, 20 40, 10 10, 90 90))", true); + } + + public void checkIsCCW(String wkt, boolean expectedValue) + { + Coordinate[] pt = read(wkt).getCoordinates(); + boolean actual = Triangle.isCCW(pt[0], pt[1], pt[2]); + assertEquals(expectedValue, actual); + } + + //=============================================================== + + public void testIntersects() { + checkIntersects("POLYGON ((30 90, 80 50, 20 20, 30 90))", "POINT (70 20)", false); + // triangle vertex + checkIntersects("POLYGON ((30 90, 80 50, 20 20, 30 90))", "POINT (30 90)", true); + checkIntersects("POLYGON ((30 90, 80 50, 20 20, 30 90))", "POINT (40 40)", true); + + // on an edge + checkIntersects("POLYGON ((30 90, 70 50, 71.5 16.5, 30 90))", "POINT (50 70)", true); + } + + public void checkIntersects(String wktTri, String wktPt, boolean expectedValue) + { + Coordinate[] tri = read(wktTri).getCoordinates(); + Coordinate pt = read(wktPt).getCoordinate(); + + boolean actual = Triangle.intersects(tri[0], tri[1], tri[2], pt); + assertEquals(expectedValue, actual); + } } diff --git a/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/ConstrainedDelaunayTriangulatorTest.java b/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/ConstrainedDelaunayTriangulatorTest.java new file mode 100644 index 0000000000..fe3e220f67 --- /dev/null +++ b/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/ConstrainedDelaunayTriangulatorTest.java @@ -0,0 +1,73 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.triangulate.polygon.ConstrainedDelaunayTriangulator; + +import junit.textui.TestRunner; +import test.jts.GeometryTestCase; + +public class ConstrainedDelaunayTriangulatorTest extends GeometryTestCase { + + public static void main(String args[]) { + TestRunner.run(ConstrainedDelaunayTriangulatorTest.class); + } + + public ConstrainedDelaunayTriangulatorTest(String name) { + super(name); + } + + public void testQuad() { + checkTri("POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 90 10, 90 90)))"); + } + + public void testPent() { + checkTri("POLYGON ((10 10, 20 40, 90 90, 100 50, 90 10, 10 10))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 10, 10 10)), POLYGON ((90 90, 20 40, 100 50, 90 90)), POLYGON ((100 50, 20 40, 90 10, 100 50)))"); + } + + public void testHoleCW() { + checkTri("POLYGON ((10 90, 90 90, 90 20, 10 10, 10 90), (30 70, 80 70, 50 30, 30 70))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 10 90, 30 70, 10 10)), POLYGON ((10 10, 30 70, 50 30, 10 10)), POLYGON ((80 70, 30 70, 90 90, 80 70)), POLYGON ((10 90, 30 70, 90 90, 10 90)), POLYGON ((80 70, 90 90, 90 20, 80 70)), POLYGON ((90 20, 10 10, 50 30, 90 20)), POLYGON ((90 20, 50 30, 80 70, 90 20)))"); + } + + public void testMultiPolygon() { + checkTri("MULTIPOLYGON (((10 10, 20 50, 50 50, 40 20, 10 10)), ((20 60, 60 60, 90 20, 90 90, 20 60)), ((10 90, 10 70, 40 70, 50 90, 10 90)))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 20 50, 40 20, 10 10)), POLYGON ((50 50, 20 50, 40 20, 50 50)), POLYGON ((90 90, 90 20, 60 60, 90 90)), POLYGON ((90 90, 60 60, 20 60, 90 90)), POLYGON ((10 70, 10 90, 40 70, 10 70)), POLYGON ((50 90, 10 90, 40 70, 50 90)))"); + } + + public void testFail() { + checkTri( + "POLYGON ((110 170, 138 272, 145 286, 152 296, 160 307, 303 307, 314 301, 332 287, 343 278, 352 270, 385 99, 374 89, 359 79, 178 89, 167 91, 153 99, 146 107, 173 157, 182 163, 191 170, 199 176, 208 184, 218 194, 226 203, 198 252, 188 247, 182 239, 175 231, 167 223, 161 213, 156 203, 155 198, 110 170))" + ); + } + + private void checkTri(String wkt, String wktExpected) { + Geometry geom = read(wkt); + Geometry actual = ConstrainedDelaunayTriangulator.triangulate(geom); + Geometry expected = read(wktExpected); + checkEqual(expected, actual); + } + + /** + * Check union of result equals original geom + * @param wkt + */ + private void checkTri(String wkt) { + Geometry geom = read(wkt); + Geometry actual = ConstrainedDelaunayTriangulator.triangulate(geom); + Geometry actualUnion = actual.union(); + checkEqual(geom, actualUnion); + } +} diff --git a/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/PolygonTriangulatorTest.java b/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/PolygonTriangulatorTest.java new file mode 100644 index 0000000000..c77c663603 --- /dev/null +++ b/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/PolygonTriangulatorTest.java @@ -0,0 +1,83 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import org.locationtech.jts.geom.Geometry; +import org.locationtech.jts.triangulate.polygon.PolygonTriangulator; + +import junit.textui.TestRunner; +import test.jts.GeometryTestCase; + +public class PolygonTriangulatorTest extends GeometryTestCase { + + public static void main(String args[]) { + TestRunner.run(PolygonTriangulatorTest.class); + } + + public PolygonTriangulatorTest(String name) { + super(name); + } + + public void testQuad() { + checkTri("POLYGON ((10 10, 20 40, 90 90, 90 10, 10 10))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 90, 10 10)), POLYGON ((90 90, 90 10, 10 10, 90 90)))"); + } + + public void testPent() { + checkTri("POLYGON ((10 10, 20 40, 90 90, 100 50, 90 10, 10 10))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 20 40, 90 90, 10 10)), POLYGON ((90 90, 100 50, 90 10, 90 90)), POLYGON ((90 10, 10 10, 90 90, 90 10)))"); + } + + public void testHoleCW() { + checkTri("POLYGON ((10 90, 90 90, 90 20, 10 10, 10 90), (30 70, 80 70, 50 30, 30 70))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 10 90, 30 70, 10 10)), POLYGON ((80 70, 30 70, 10 90, 80 70)), POLYGON ((10 10, 30 70, 50 30, 10 10)), POLYGON ((80 70, 10 90, 90 90, 80 70)), POLYGON ((90 20, 10 10, 50 30, 90 20)), POLYGON ((80 70, 90 90, 90 20, 80 70)), POLYGON ((90 20, 50 30, 80 70, 90 20)))"); + } + + public void testTouchingHoles() { + checkTri("POLYGON ((10 90, 90 90, 90 10, 10 10, 10 90), (20 80, 50 70, 30 30, 20 80), (70 20, 50 70, 80 80, 70 20))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 10 90, 20 80, 10 10)), POLYGON ((50 70, 20 80, 10 90, 50 70)), POLYGON ((10 10, 20 80, 30 30, 10 10)), POLYGON ((30 30, 50 70, 70 20, 30 30)), POLYGON ((80 80, 50 70, 10 90, 80 80)), POLYGON ((90 10, 10 10, 30 30, 90 10)), POLYGON ((80 80, 10 90, 90 90, 80 80)), POLYGON ((90 10, 30 30, 70 20, 90 10)), POLYGON ((70 20, 80 80, 90 90, 70 20)), POLYGON ((90 90, 90 10, 70 20, 90 90)))"); + } + + public void testRepeatedPoints() { + checkTri("POLYGON ((71 195, 178 335, 178 335, 239 185, 380 210, 290 60, 110 70, 71 195))" + ,"GEOMETRYCOLLECTION (POLYGON ((71 195, 178 335, 239 185, 71 195)), POLYGON ((71 195, 239 185, 290 60, 71 195)), POLYGON ((71 195, 290 60, 110 70, 71 195)), POLYGON ((239 185, 380 210, 290 60, 239 185)))"); + } + + public void testMultiPolygon() { + checkTri("MULTIPOLYGON (((10 10, 20 50, 50 50, 40 20, 10 10)), ((20 60, 60 60, 90 20, 90 90, 20 60)), ((10 90, 10 70, 40 70, 50 90, 10 90)))" + ,"GEOMETRYCOLLECTION (POLYGON ((10 10, 20 50, 50 50, 10 10)), POLYGON ((50 50, 40 20, 10 10, 50 50)), POLYGON ((90 90, 90 20, 60 60, 90 90)), POLYGON ((60 60, 20 60, 90 90, 60 60)), POLYGON ((10 70, 10 90, 50 90, 10 70)), POLYGON ((50 90, 40 70, 10 70, 50 90)))"); + } + + public void testCeeShape() { + checkTri( + "POLYGON ((110 170, 138 272, 145 286, 152 296, 160 307, 303 307, 314 301, 332 287, 343 278, 352 270, 385 99, 374 89, 359 79, 178 89, 167 91, 153 99, 146 107, 173 157, 182 163, 191 170, 199 176, 208 184, 218 194, 226 203, 198 252, 188 247, 182 239, 175 231, 167 223, 161 213, 156 203, 155 198, 110 170))" + ); + } + + private void checkTri(String wkt, String wktExpected) { + Geometry geom = read(wkt); + Geometry actual = PolygonTriangulator.triangulate(geom); + Geometry expected = read(wktExpected); + checkEqual(expected, actual); + } + + /** + * Check union of result equals original geom + * @param wkt + */ + private void checkTri(String wkt) { + Geometry geom = read(wkt); + Geometry actual = PolygonTriangulator.triangulate(geom); + Geometry actualUnion = actual.union(); + checkEqual(geom, actualUnion); + } +} diff --git a/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/VertexSequencePackedRtreeTest.java b/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/VertexSequencePackedRtreeTest.java new file mode 100644 index 0000000000..d6e1ed237e --- /dev/null +++ b/modules/core/src/test/java/org/locationtech/jts/triangulate/polygon/VertexSequencePackedRtreeTest.java @@ -0,0 +1,88 @@ +/* + * Copyright (c) 2021 Martin Davis. + * + * All rights reserved. This program and the accompanying materials + * are made available under the terms of the Eclipse Public License 2.0 + * and Eclipse Distribution License v. 1.0 which accompanies this distribution. + * The Eclipse Public License is available at http://www.eclipse.org/legal/epl-v20.html + * and the Eclipse Distribution License is available at + * + * http://www.eclipse.org/org/documents/edl-v10.php. + */ +package org.locationtech.jts.triangulate.polygon; + +import org.locationtech.jts.geom.Coordinate; +import org.locationtech.jts.geom.Envelope; +import org.locationtech.jts.triangulate.polygon.VertexSequencePackedRtree; + +import junit.framework.TestCase; +import junit.textui.TestRunner; + +public class VertexSequencePackedRtreeTest extends TestCase { + public static void main(String args[]) { + TestRunner.run(VertexSequencePackedRtreeTest.class); + } + + public VertexSequencePackedRtreeTest(String name) { + super(name); + } + + public void test1() { + VertexSequencePackedRtree tree = createSPRtree(1,1); + checkQuery(tree, 1,1,4,4, result( 0 )); + } + + public void test2() { + VertexSequencePackedRtree tree = createSPRtree(0,0, 1,1); + checkQuery(tree, 1,1,4,4, result( 1 )); + } + + public void test6() { + VertexSequencePackedRtree tree = createSPRtree(0,0, 1,1, 2,2, 3,3, 4,4, 5,5); + checkQuery(tree, 2,2,4,4, result( 2,3,4 )); + checkQuery(tree, 0,0,0,0, result( 0 )); + } + + public void test10() { + VertexSequencePackedRtree tree = createSPRtree(0,0, 1,1, 2,2, 3,3, 4,4, 5,5, 6,6, 7,7, 8,8, 9,9, 10,10); + checkQuery(tree, 2,2,4,4, result( 2,3,4 )); + checkQuery(tree, 7,7,8,8, result( 7,8 )); + checkQuery(tree, 0,0,0,0, result( 0 )); + } + + public void test6WithDups() { + VertexSequencePackedRtree tree = createSPRtree(0,0, 1,1, 2,2, 3,3, 4,4, 5,5, 4,4, 3,3, 2,2, 1,1, 0,0); + checkQuery(tree, 2,2,4,4, result( 2,3,4, 6, 7, 8 )); + checkQuery(tree, 0,0,0,0, result( 0, 10 )); + } + + private void checkQuery(VertexSequencePackedRtree tree, + double xmin, double ymin, double xmax, double ymax, int[] expected) { + Envelope env = new Envelope(xmin, xmax, ymin, ymax); + int[] result = tree.query(env); + assertEquals(expected.length, result.length); + assertTrue( isEqualResult(expected, result) ); + } + + private boolean isEqualResult(int[] expected, int[] result) { + for (int i = 0; i < result.length; i++) { + if (expected[i] != result[i]) + return false; + + } + return true; + } + + private int[] result(int... i) { + return i; + } + + private VertexSequencePackedRtree createSPRtree(int... ords) { + int numCoord = ords.length / 2; + Coordinate[] pt = new Coordinate[numCoord]; + for (int i = 0 ; i < numCoord; i++) { + pt[i] = new Coordinate(ords[2*i], ords[2*i+1]); + } + return new VertexSequencePackedRtree(pt); + } +} diff --git a/modules/core/src/test/java/org/locationtech/jts/triangulate/tri/TriTest.java b/modules/core/src/test/java/org/locationtech/jts/triangulate/tri/TriTest.java new file mode 100644 index 0000000000..466f7efbd4 --- /dev/null +++ b/modules/core/src/test/java/org/locationtech/jts/triangulate/tri/TriTest.java @@ -0,0 +1,73 @@ +package org.locationtech.jts.triangulate.tri; + +import java.util.ArrayList; +import java.util.List; + +import org.locationtech.jts.algorithm.Orientation; +import org.locationtech.jts.geom.Coordinate; + +import junit.textui.TestRunner; +import test.jts.GeometryTestCase; + +public class TriTest extends GeometryTestCase { + + public static void main(String args[]) { + TestRunner.run(TriTest.class); + } + + private static Tri triCentre = createSimpleTriangulation(); + private static Tri tri0; + private static Tri tri1; + private static Tri tri2; + + public TriTest(String name) { + super(name); + } + + public void testAdjacent() { + assertTrue(tri0 == triCentre.getAdjacent(0)); + assertTrue(tri1 == triCentre.getAdjacent(1)); + assertTrue(tri2 == triCentre.getAdjacent(2)); + } + + public void testMidpoint() { + Tri tri = tri(0,0, 0,10, 10,0 ); + checkEqualXY(new Coordinate(0, 5), tri.midpoint(0)); + checkEqualXY(new Coordinate(5, 5), tri.midpoint(1)); + checkEqualXY(new Coordinate(5, 0), tri.midpoint(2)); + } + + public void testCoordinateIndex() { + Tri tri = tri(0,0, 0,10, 10,0 ); + assertEquals(0, tri.getIndex(new Coordinate(0,0))); + assertEquals(1, tri.getIndex(new Coordinate(0,10))); + assertEquals(2, tri.getIndex(new Coordinate(10,0))); + } + + private static Tri tri(double x0, double y0, double x1, double y1, double x2, double y2) { + Tri tri = Tri.create( + new Coordinate(x0, y0), + new Coordinate(x1, y1), + new Coordinate(x2, y2)); + assertTrue( Orientation.CLOCKWISE == Orientation.index( + tri.getCoordinate(0), tri.getCoordinate(1), tri.getCoordinate(2) )); + return tri; + } + + private static Tri createSimpleTriangulation() { + Tri tri = tri(10,10, 10,20, 20,10 ); + tri0 = tri(10,20, 10,10, 0,10 ); + tri1 = tri(20,10, 10,20, 20,20 ); + tri2 = tri(10,10, 20,10, 10,0 ); + build(tri, tri0, tri1, tri2); + return tri; + } + + private static void build(Tri... tri) { + List triList = new ArrayList(); + for (int i = 0; i < tri.length; i++) { + triList.add(tri[i]); + } + TriangulationBuilder.build(triList); + } +}