You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by js...@apache.org on 2026/10/21 07:02:11 UTC
(sis) 09/13: feat(Geometry): improve earclipping, add contains and intersections tests
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 a81b6f8cf270b81e947a9b5a5067a5b755a05487
Author: jsorel <jo...@geomatys.com>
AuthorDate: Mon Oct 19 10:11:02 2026 +0200
feat(Geometry): improve earclipping, add contains and intersections tests
---
.../apache/sis/geometries/operation/Contains.java | 141 +++++-
.../sis/geometries/operation/Intersection.java | 549 +++++++++++++++++++++
.../operation/triangulate/EarClipping.java | 65 ++-
.../main/org/apache/sis/maths/Maths.java | 28 ++
.../sis/geometries/operation/Contains3DTest.java | 182 +++++++
.../geometries/operation/Intersection3DTest.java | 264 ++++++++++
.../operation/triangulate/EarClippingTest.java | 42 ++
.../test/org/apache/sis/maths/MathsTest.java | 28 ++
8 files changed, 1258 insertions(+), 41 deletions(-)
diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Contains.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Contains.java
index 293338f1bf..ec0c60ddfc 100644
--- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Contains.java
+++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Contains.java
@@ -18,56 +18,72 @@ package org.apache.sis.geometries.operation;
import org.apache.sis.geometries.Point;
import org.apache.sis.geometries.DataPointsType;
+import org.apache.sis.geometries.GeometryCollection;
import org.apache.sis.geometries.curve.LineString;
+import org.apache.sis.geometries.operation.triangulate.EarClipping;
+import org.apache.sis.geometries.solid.Polyhedron;
import org.apache.sis.geometries.surface.Polygon;
import org.apache.sis.maths.Array;
import org.apache.sis.maths.Cursor;
import org.apache.sis.maths.Maths;
-import org.apache.sis.maths.Vectors;
+import org.apache.sis.maths.ReadOnly;
import org.apache.sis.maths.Vector;
/**
- * Constains 2D processors.
+ * Contains processors.
*
* @author Johann Sorel (Geomatys)
*/
public final class Contains {
+ /**
+ * The axes of the (x,y) plane, for the tests done in two dimensions.
+ */
+ private static final int[] XY = {0, 1};
+
private Contains(){}
/**
- * Test if point is within polygon using Winding Number algorithm.
+ * Test if point is within a ring using Winding Number algorithm, in the plane of the two given axes.
+ * The ring and the point can have more dimensions: for example a ring in a vertical plane can be tested
+ * in the plane given by {@link Maths#projectionAxes(double[])}. A point on the ring gets an arbitrary answer.
* http://geomalgorithms.com/a03-_inclusion.html
* http://en.wikipedia.org/wiki/Point_in_polygon
*
- * @param ring not null
- * @param point not null
+ * @param ring positions of the closed ring, the last one at the location of the first one
+ * @param point the point to test
+ * @param axes indices of the two coordinates of the plane in which the test is done
+ * @return whether the point is inside the ring
*/
- private static boolean contains(Array ring, Vector<?> point) {
+ public static boolean contains(Array ring, ReadOnly.Vector<?> point, int[] axes) {
+ final long n = ring.getLength();
+ if (n == 0) {
+ return false;
+ }
+ final int u = axes[0], v = axes[1];
+ final double px = point.get(u), py = point.get(v);
final Cursor cursor = ring.cursor();
-
- int windingNumber = 0;
- Vector<?> current;
- Vector<?> previous;
cursor.moveTo(0);
- current = cursor.samples();
- previous = Vectors.create(current.getSampleSystem(), current.getDataType());
- final double pointY = point.get(1);
- for (long i = 1, n = ring.getLength(); i < n; i++){
- previous.set(current);
+ Vector<?> current = cursor.samples();
+ double x = current.get(u), y = current.get(v);
+ int windingNumber = 0;
+ for (long i = 1; i < n; i++){
+ final double previousX = x, previousY = y;
cursor.moveTo(i);
current = cursor.samples();
+ x = current.get(u);
+ y = current.get(v);
- if (previous.get(1) <= pointY){
- if (current.get(1) > pointY){
- if (Maths.lineSide(previous, current, point) > 0){
+ if (previousY <= py){
+ if (y > py){
+ if (Maths.lineSide(previousX, previousY, x, y, px, py) > 0){
windingNumber++;
}
}
} else {
- if (current.get(1) <= pointY){
- if (Maths.lineSide(previous, current, point) < 0){
+ if (y <= py){
+ if (Maths.lineSide(previousX, previousY, x, y, px, py) < 0){
windingNumber--;
}
}
@@ -78,6 +94,87 @@ public final class Contains {
return windingNumber != 0;
}
+ /**
+ * Returns the generalized winding number of a surface made of triangles around a point: the sum of the solid
+ * angles of the triangles seen from the point, divided by 4π. The number is 1 inside a closed surface whose
+ * triangles are oriented outward (counterclockwise seen from outside), −1 if they are oriented inward, and 0
+ * outside. The numbers of several surfaces add: a point inside the void of a solid is inside the outward
+ * surface of the solid and inside the inward surface of the void, which gives 0. The number of a surface
+ * which is not closed is between those values, and the number of a point on a surface is a fraction,
+ * for example ½ on a face.
+ *
+ * <p>The solid angle of each triangle is computed with the formula of Van Oosterom and Strackee.
+ * A triangle whose plane contains the point is ignored: seen from its plane, a triangle has no solid
+ * angle, unless the point is on it, which is the case of a point on the surface.</p>
+ *
+ * @param xyz the coordinates of the vertices, three per vertex.
+ * @param triangles indices of the three vertices of each triangle, one triangle after the other.
+ * @param point the point, of three coordinates.
+ * @return the winding number, close to an integer for a point which is not on a closed surface.
+ *
+ * @see <a href="https://doi.org/10.1145/2461912.2461916">Jacobson et al., Robust inside-outside segmentation using generalized winding numbers</a>
+ */
+ public static double windingNumber(final double[] xyz, final int[] triangles, final double[] point) {
+ double sum = 0;
+ for (int i = 0; i + 2 < triangles.length; i += 3) {
+ final int a = triangles[i] * 3, b = triangles[i+1] * 3, c = triangles[i+2] * 3;
+ final double ax = xyz[a] - point[0], ay = xyz[a+1] - point[1], az = xyz[a+2] - point[2];
+ final double bx = xyz[b] - point[0], by = xyz[b+1] - point[1], bz = xyz[b+2] - point[2];
+ final double cx = xyz[c] - point[0], cy = xyz[c+1] - point[1], cz = xyz[c+2] - point[2];
+ final double la = Math.sqrt(ax*ax + ay*ay + az*az);
+ final double lb = Math.sqrt(bx*bx + by*by + bz*bz);
+ final double lc = Math.sqrt(cx*cx + cy*cy + cz*cz);
+ final double det = ax * (by*cz - bz*cy) + ay * (bz*cx - bx*cz) + az * (bx*cy - by*cx);
+ if (det != 0) {
+ final double den = la*lb*lc + (ax*bx + ay*by + az*bz) * lc + (ax*cx + ay*cy + az*cz) * lb + (bx*cx + by*cy + bz*cz) * la;
+ sum += Math.atan2(det, den); // half of the solid angle.
+ }
+ }
+ return sum / (2 * Math.PI);
+ }
+
+ /**
+ * Polyhedron contains Point test. The point is inside the polyhedron if it is inside the exterior shell and
+ * outside all the interior shells, as told by the {@linkplain #windingNumber winding number} of each shell.
+ * Each shell is tested alone, so the orientation of its faces does not matter as long as it is consistent.
+ * A point on the boundary of the polyhedron gets an arbitrary answer.
+ *
+ * @throws OperationException if the geometries do not have the same CRS, or are not three-dimensional,
+ * or if a face cannot be triangulated.
+ */
+ public static boolean contains(Polyhedron polyhedron, Point candidate) throws OperationException {
+ ProcessorUtils.ensureSameCRS(polyhedron, candidate);
+ final double[] point = candidate.getPosition().toArrayDouble();
+ if (point.length != 3) {
+ throw new OperationException("A polyhedron contains a point in 3 dimensions only.");
+ }
+ if (!isInside(polyhedron.getExteriorShell(), point)) {
+ return false;
+ }
+ for (final GeometryCollection<Polygon> shell : polyhedron.getInteriorShells()) {
+ if (isInside(shell, point)) {
+ //point is within a void
+ return false;
+ }
+ }
+ return true;
+ }
+
+ /**
+ * Returns whether the given point is inside the given closed shell, whatever the orientation of its faces.
+ * The winding numbers of the faces are added, each face being split in triangles.
+ */
+ private static boolean isInside(final GeometryCollection<Polygon> shell, final double[] point) {
+ final EarClipping triangulator = new EarClipping();
+ double winding = 0;
+ for (int i = 0, n = shell.getNumGeometries(); i < n; i++) {
+ final EarClipping.Triangulation face = triangulator.triangulateIndexed(shell.getGeometryN(i));
+ final double[] xyz = face.vertices().getAttributeArray(DataPointsType.ATT_POSITION).toArrayDouble();
+ winding += windingNumber(xyz, face.triangles(), point);
+ }
+ return Math.abs(winding) > 0.5;
+ }
+
/**
* Polygon contains Point test.
*/
@@ -86,7 +183,7 @@ public final class Contains {
{ //check exterior
final Array coords = polygon.getExteriorRing().getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION);
- if (!contains(coords, candidate.getPosition())) {
+ if (!contains(coords, candidate.getPosition(), XY)) {
//point is outside the exterior ring
return false;
}
@@ -96,7 +193,7 @@ public final class Contains {
for (int i = 0, n = polygon.getNumInteriorRing(); i < n; i++) {
final LineString hole = polygon.getInteriorRingN(i);
final Array coords = hole.getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION);
- if (contains(coords, candidate.getPosition())) {
+ if (contains(coords, candidate.getPosition(), XY)) {
//point is within a hole
return false;
}
diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java
index db21ba0719..44526319ae 100644
--- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java
+++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java
@@ -18,6 +18,7 @@ package org.apache.sis.geometries.operation;
import java.util.ArrayList;
import java.util.Arrays;
+import java.util.Comparator;
import java.util.Iterator;
import java.util.List;
import java.util.Optional;
@@ -39,10 +40,12 @@ import org.apache.sis.geometries.surface.PreparedTIN;
import org.apache.sis.geometries.surface.Triangle;
import org.apache.sis.maths.Array;
import org.apache.sis.maths.Cursor;
+import org.apache.sis.maths.Maths;
import org.apache.sis.maths.Matrices;
import org.apache.sis.maths.NDArrays;
import org.apache.sis.maths.ReadOnly;
import org.apache.sis.maths.Vector;
+import org.apache.sis.maths.Vectors;
import static org.apache.sis.maths.Maths.clamp;
@@ -208,6 +211,552 @@ public final class Intersection {
return intersection;
}
+ /**
+ * Intersection of two segments in a plane. Each point of the intersection is given by its position on both
+ * segments, as a ratio from 0 at the start of the segment to 1 at its end. A ratio is exactly 0 or 1 when the
+ * point is an end of the segment.
+ *
+ * @param type whether the segments are disjoint, meet at one point, or overlap along their common line.
+ * @param startOnFirst ratio on the first segment of the single point, or of the start of the overlap, NaN if none.
+ * @param endOnFirst ratio on the first segment of the end of the overlap, greater than {@code startOnFirst},
+ * or {@code startOnFirst} for a single point.
+ * @param startOnSecond ratio on the second segment of the point at {@code startOnFirst}.
+ * @param endOnSecond ratio on the second segment of the point at {@code endOnFirst}.
+ *
+ * @see #segmentSegment2D(double[], double[], double[], double[], int[])
+ */
+ public record SegmentSegment(Type type, double startOnFirst, double endOnFirst, double startOnSecond, double endOnSecond) {
+ /**
+ * Kinds of intersection of two segments.
+ */
+ public enum Type {
+ /** The segments do not intersect. */
+ NONE,
+ /** The segments cross or touch at a single point. */
+ POINT,
+ /** The segments are on the same line and share a part of it. */
+ OVERLAP
+ }
+
+ /**
+ * The intersection of disjoint segments.
+ */
+ public static final SegmentSegment NONE = new SegmentSegment(Type.NONE, Double.NaN, Double.NaN, Double.NaN, Double.NaN);
+
+ /**
+ * Returns the intersection made of a single point, at the given ratios on the two segments.
+ */
+ static SegmentSegment point(final double onFirst, final double onSecond) {
+ return new SegmentSegment(Type.POINT, onFirst, onFirst, onSecond, onSecond);
+ }
+ }
+
+ /**
+ * Point where a segment touches a triangle, located on both. Edge <var>i</var> of the triangle goes from
+ * its corner <var>i</var> to its corner (<var>i</var>+1) modulo 3.
+ *
+ * @param type where the point is on the triangle, or why there is no point.
+ * @param feature index of the edge or of the corner of the triangle on which the point is, or -1.
+ * @param ratio position of the point on the segment, from 0 at its start to 1 at its end,
+ * exactly 0 or 1 if the point is an end of the segment. NaN if there is no point.
+ *
+ * @see #segmentTriangle(double[], double[], double[], double[], double[], double)
+ * @see #segmentInTriangle2D(double[], double[], double[], double[], double[], int[])
+ */
+ public record SegmentTriangle(Type type, int feature, double ratio) {
+ /**
+ * Where a segment touches a triangle.
+ */
+ public enum Type {
+ /** The segment does not touch the triangle. */
+ NONE,
+ /** The segment lies in the plane of the triangle, and shall be intersected with it in that plane. */
+ COPLANAR,
+ /** The segment touches the triangle inside it. */
+ INTERIOR,
+ /** The segment touches the triangle on the edge given by the feature index. */
+ EDGE,
+ /** The segment touches the triangle on the corner given by the feature index. */
+ VERTEX
+ }
+
+ /**
+ * A segment which does not touch the triangle.
+ */
+ public static final SegmentTriangle NONE = new SegmentTriangle(Type.NONE, -1, Double.NaN);
+
+ /**
+ * A segment in the plane of the triangle.
+ */
+ public static final SegmentTriangle COPLANAR = new SegmentTriangle(Type.COPLANAR, -1, Double.NaN);
+
+ /**
+ * Returns the point inside the triangle, on an edge or on a corner, from the number of edges on which
+ * the point is. This is used after testing the side of the point relative to the three edges.
+ *
+ * @param zeros number of edges on which the point is.
+ * @param zeroEdge index of an edge on which the point is.
+ * @param otherEdge index of an edge on which the point is not.
+ * @param ratio position of the point on the segment.
+ * @return the point, or {@code null} if it is on the three edges, which happens with a degenerate triangle.
+ */
+ static SegmentTriangle locate(final int zeros, final int zeroEdge, final int otherEdge, final double ratio) {
+ switch (zeros) {
+ case 0: return new SegmentTriangle(Type.INTERIOR, -1, ratio);
+ case 1: return new SegmentTriangle(Type.EDGE, zeroEdge, ratio);
+ case 2: return new SegmentTriangle(Type.VERTEX, (otherEdge + 2) % 3, ratio); // corner shared by the two other edges.
+ default: return null;
+ }
+ }
+ }
+
+ /**
+ * End of the intersection of two triangles: a point where an edge of one triangle touches the other triangle.
+ *
+ * @param edgeOfFirst whether the edge belongs to the first triangle, the point being then on the second triangle.
+ * @param edge index of the edge, which goes from the corner of same index to the next corner.
+ * @param contact where the edge touches the other triangle, with the position of the point on the edge.
+ * @param position coordinates of the point, interpolated on the edge.
+ *
+ * @see #triangleTriangle(double[], double[], double[], double[], double[], double[], double)
+ */
+ public record TriangleEdgeContact(boolean edgeOfFirst, int edge, SegmentTriangle contact, double[] position) {
+ }
+
+ /**
+ * Intersection of two triangles: nothing, a point or a segment, unless the triangles are in the same plane.
+ *
+ * @param type kind of intersection.
+ * @param start first end of the intersection, or {@code null} if the triangles are disjoint or coplanar.
+ * @param end second end of the intersection, the same as {@code start} for a single point.
+ *
+ * @see #triangleTriangle(double[], double[], double[], double[], double[], double[], double)
+ */
+ public record TriangleTriangle(Type type, TriangleEdgeContact start, TriangleEdgeContact end) {
+ /**
+ * Kinds of intersection of two triangles.
+ */
+ public enum Type {
+ /** The triangles do not intersect. */
+ DISJOINT,
+ /** The triangles touch at a single point. */
+ POINT,
+ /** The triangles cross or touch along a segment. */
+ SEGMENT,
+ /** The triangles are in the same plane, and shall be intersected in that plane. */
+ COPLANAR
+ }
+
+ /** Intersection of disjoint triangles. */
+ public static final TriangleTriangle DISJOINT = new TriangleTriangle(Type.DISJOINT, null, null);
+
+ /** Intersection of triangles in the same plane. */
+ public static final TriangleTriangle COPLANAR = new TriangleTriangle(Type.COPLANAR, null, null);
+ }
+
+ /**
+ * Computes the intersection of two segments in the plane of the two given axes. The side tests are robust,
+ * so whether the segments intersect, and whether a point of the intersection is an end of a segment, are exact.
+ * Only the ratios of the points inside the segments are approximations.
+ *
+ * @param a0 start of the first segment.
+ * @param a1 end of the first segment.
+ * @param b0 start of the second segment.
+ * @param b1 end of the second segment.
+ * @param axes indices of the two coordinates of the plane, for example from {@link Maths#projectionAxes(double[])}.
+ * @return the intersection of the two segments.
+ */
+ public static SegmentSegment segmentSegment2D(double[] a0, double[] a1, double[] b0, double[] b1, int[] axes) {
+ final int u = axes[0], v = axes[1];
+ //sides of the ends of each segment relative to the line of the other segment.
+ final double d0 = side(b0, b1, a0, u, v), d1 = side(b0, b1, a1, u, v);
+ final double e0 = side(a0, a1, b0, u, v), e1 = side(a0, a1, b1, u, v);
+ if (sameSide(d0, d1) || sameSide(e0, e1)) {
+ return SegmentSegment.NONE;
+ }
+ if ((d0 == 0 && d1 == 0) || (e0 == 0 && e1 == 0)) {
+ return collinear(a0, a1, b0, b1, u, v);
+ }
+ return SegmentSegment.point(ratio(d0, d1), ratio(e0, e1));
+ }
+
+ /**
+ * Computes the intersection of two segments on the same line. The positions are compared along the axis
+ * where the segments are the longest, so that a line parallel to one of the axes is not degenerate.
+ */
+ private static SegmentSegment collinear(double[] a0, double[] a1, double[] b0, double[] b1, int u, int v) {
+ final int k = (Math.abs(a1[u] - a0[u]) + Math.abs(b1[u] - b0[u]) >= Math.abs(a1[v] - a0[v]) + Math.abs(b1[v] - b0[v])) ? u : v;
+ final double lower = Math.max(Math.min(a0[k], a1[k]), Math.min(b0[k], b1[k]));
+ final double upper = Math.min(Math.max(a0[k], a1[k]), Math.max(b0[k], b1[k]));
+ if (!(lower <= upper)) {
+ return SegmentSegment.NONE;
+ }
+ double startOnFirst = along(a0[k], a1[k], lower), startOnSecond = along(b0[k], b1[k], lower);
+ if (lower == upper) {
+ return SegmentSegment.point(startOnFirst, startOnSecond);
+ }
+ double endOnFirst = along(a0[k], a1[k], upper), endOnSecond = along(b0[k], b1[k], upper);
+ if (startOnFirst > endOnFirst) {
+ //the first segment goes toward decreasing values along the axis.
+ double t = startOnFirst; startOnFirst = endOnFirst; endOnFirst = t;
+ t = startOnSecond; startOnSecond = endOnSecond; endOnSecond = t;
+ }
+ return new SegmentSegment(SegmentSegment.Type.OVERLAP, startOnFirst, endOnFirst, startOnSecond, endOnSecond);
+ }
+
+ /**
+ * Returns the position of a value between the values at the start and at the end of a segment along an axis,
+ * exactly 0 or 1 at the ends.
+ */
+ private static double along(final double start, final double end, final double value) {
+ if (value == start) return 0;
+ if (value == end) return 1;
+ return (value - start) / (end - start);
+ }
+
+ /**
+ * Returns the position of the point where a segment crosses a line or a plane, from the signed distances of
+ * the segment ends to that line or plane, or any values proportional to those distances. The position is
+ * exactly 0 or 1 if an end is on the line or plane.
+ */
+ private static double ratio(final double start, final double end) {
+ if (start == 0) return 0;
+ if (end == 0) return 1;
+ return clamp(start / (start - end), 0, 1);
+ }
+
+ /**
+ * Returns whether two values are both strictly positive or both strictly negative.
+ */
+ private static boolean sameSide(final double s1, final double s2) {
+ return (s1 > 0 && s2 > 0) || (s1 < 0 && s2 < 0);
+ }
+
+ /**
+ * Returns whether one value is strictly positive and the other one strictly negative.
+ */
+ private static boolean opposite(final double s1, final double s2) {
+ return (s1 > 0 && s2 < 0) || (s1 < 0 && s2 > 0);
+ }
+
+ /**
+ * Returns the robust side of point r relative to the line from p to q, in the plane of the axes u and v.
+ * Positive on the left, negative on the right, zero on the line.
+ */
+ private static double side(final double[] p, final double[] q, final double[] r, final int u, final int v) {
+ return Maths.lineSide(p[u], p[v], q[u], q[v], r[u], r[v]);
+ }
+
+ /**
+ * Returns the normal (b−a)×(c−a) of a triangle, whose length is twice the area of the triangle.
+ */
+ private static double[] normal(final double[] a, final double[] b, final double[] c) {
+ final double ux = b[0] - a[0], uy = b[1] - a[1], uz = b[2] - a[2];
+ final double vx = c[0] - a[0], vy = c[1] - a[1], vz = c[2] - a[2];
+ return new double[] {uy*vz - uz*vy, uz*vx - ux*vz, ux*vy - uy*vx};
+ }
+
+ /**
+ * Computes where a segment touches a triangle, in three dimensions. An end of the segment within the tolerance
+ * of the plane of the triangle is on that plane, as decided by {@link Maths#planeSide Maths.planeSide(…)}.
+ * If both ends are on the plane, the segment is {@linkplain SegmentTriangle.Type#COPLANAR coplanar} and its
+ * intersection with the triangle is given by {@link #segmentInTriangle2D segmentInTriangle2D(…)}.
+ *
+ * <p>When the segment crosses the plane, where it crosses the triangle is decided by the robust sides of
+ * its line relative to the three edges, so that a segment crossing an edge or a corner is found exactly on it.
+ * Only the ratio of a point inside the segment is an approximation. When an end of the segment is on the plane,
+ * that end is located in the triangle with robust side tests in the plane.</p>
+ *
+ * @param s0 start of the segment.
+ * @param s1 end of the segment.
+ * @param a first corner of the triangle.
+ * @param b second corner of the triangle.
+ * @param c third corner of the triangle.
+ * @param tolerance largest distance to the plane of the triangle for considering an end of the segment on it.
+ * @return where the segment touches the triangle.
+ */
+ public static SegmentTriangle segmentTriangle(double[] s0, double[] s1, double[] a, double[] b, double[] c, double tolerance) {
+ final int side0 = Maths.planeSide(a, b, c, s0, tolerance);
+ final int side1 = Maths.planeSide(a, b, c, s1, tolerance);
+ if (side0 == 0 && side1 == 0) {
+ return SegmentTriangle.COPLANAR;
+ }
+ if (side0 == side1) {
+ return SegmentTriangle.NONE; // both ends strictly on the same side.
+ }
+ final double[][] corners = {a, b, c};
+ final double[] normal = normal(a, b, c);
+ if (side0 == 0 || side1 == 0) {
+ final int[] axes = Maths.projectionAxes(normal);
+ return (side0 == 0) ? locate2D(s0, corners, axes[0], axes[1], 0)
+ : locate2D(s1, corners, axes[0], axes[1], 1);
+ }
+ /*
+ * The segment crosses the plane between its ends. The sign of the tetrahedron made by the segment and
+ * an edge tells on which side of that edge the line of the segment passes. The line passes inside the
+ * triangle if the three signs are the same, and on an edge or a corner if one or two of them are zero.
+ */
+ int positive = 0, negative = 0, zeros = 0, zeroEdge = -1, otherEdge = -1;
+ for (int i = 0; i < 3; i++) {
+ final double[] p = corners[i], q = corners[(i + 1) % 3];
+ final double sign = ShewchukPredicates.orient3d(s0[0], s0[1], s0[2], s1[0], s1[1], s1[2],
+ p[0], p[1], p[2], q[0], q[1], q[2]);
+ if (sign > 0) {positive++; otherEdge = i;}
+ else if (sign < 0) {negative++; otherEdge = i;}
+ else {zeros++; zeroEdge = i;}
+ }
+ if (positive != 0 && negative != 0) {
+ return SegmentTriangle.NONE;
+ }
+ final double d0 = normal[0] * (s0[0] - a[0]) + normal[1] * (s0[1] - a[1]) + normal[2] * (s0[2] - a[2]);
+ final double d1 = normal[0] * (s1[0] - a[0]) + normal[1] * (s1[1] - a[1]) + normal[2] * (s1[2] - a[2]);
+ final SegmentTriangle point = SegmentTriangle.locate(zeros, zeroEdge, otherEdge, ratio(d0, d1));
+ return (point != null) ? point : SegmentTriangle.COPLANAR;
+ }
+
+ /**
+ * Locates a point of the plane of a triangle relative to that triangle, with robust side tests in the plane
+ * of the two given axes.
+ *
+ * @param p the point to locate.
+ * @param corners the three corners of the triangle.
+ * @param ratio position of the point on the segment for which the point is located.
+ * @return where the point is on the triangle, or {@link SegmentTriangle#NONE} if outside or if the triangle is degenerate.
+ */
+ private static SegmentTriangle locate2D(final double[] p, final double[][] corners, final int u, final int v, final double ratio) {
+ final double orientation = Math.signum(side(corners[0], corners[1], corners[2], u, v));
+ if (orientation == 0) {
+ return SegmentTriangle.NONE;
+ }
+ int zeros = 0, zeroEdge = -1, otherEdge = -1;
+ for (int i = 0; i < 3; i++) {
+ final double s = orientation * side(corners[i], corners[(i + 1) % 3], p, u, v);
+ if (s < 0) {
+ return SegmentTriangle.NONE;
+ }
+ if (s == 0) {zeros++; zeroEdge = i;}
+ else {otherEdge = i;}
+ }
+ final SegmentTriangle point = SegmentTriangle.locate(zeros, zeroEdge, otherEdge, ratio);
+ return (point != null) ? point : SegmentTriangle.NONE;
+ }
+
+ /**
+ * Computes the part of a segment inside a triangle, both in the same plane, with robust side tests in the plane
+ * of the two given axes. The result gives the two ends of that part, each with its position on the segment and
+ * its location on the triangle: an end is an end of the segment inside the triangle, or a point where the
+ * segment crosses an edge or passes through a corner. Whether a point is on an edge or a corner, or is an end
+ * of the segment, is exact. Only the ratios of the points inside the segment are approximations.
+ *
+ * @param s0 start of the segment.
+ * @param s1 end of the segment.
+ * @param a first corner of the triangle.
+ * @param b second corner of the triangle.
+ * @param c third corner of the triangle.
+ * @param axes indices of the two coordinates of the plane, for example from {@link Maths#projectionAxes(double[])}.
+ * @return the ends of the part of the segment inside the triangle by increasing ratio: none if the segment
+ * does not touch the triangle, one if it touches it at a single point, two otherwise.
+ * Their types are {@code INTERIOR}, {@code EDGE} or {@code VERTEX}.
+ */
+ public static List<SegmentTriangle> segmentInTriangle2D(double[] s0, double[] s1, double[] a, double[] b, double[] c, int[] axes) {
+ final int u = axes[0], v = axes[1];
+ final double[][] corners = {a, b, c};
+ if (side(a, b, c, u, v) == 0) {
+ return List.of(); // a degenerate triangle has no interior.
+ }
+ final List<SegmentTriangle> points = new ArrayList<>(4);
+ final SegmentTriangle start = locate2D(s0, corners, u, v, 0);
+ if (start != SegmentTriangle.NONE) {
+ points.add(start);
+ }
+ final boolean isPoint = (s0[u] == s1[u] && s0[v] == s1[v]);
+ if (!isPoint) {
+ final SegmentTriangle end = locate2D(s1, corners, u, v, 1);
+ if (end != SegmentTriangle.NONE) {
+ points.add(end);
+ }
+ /*
+ * The corners on the line of the segment, between its ends.
+ * The ends of the segment have been located above.
+ */
+ final double dx = s1[u] - s0[u], dy = s1[v] - s0[v];
+ final double[] sides = new double[3];
+ for (int j = 0; j < 3; j++) {
+ final double[] p = corners[j];
+ sides[j] = side(s0, s1, p, u, v);
+ if (sides[j] == 0 && !(p[u] == s0[u] && p[v] == s0[v]) && !(p[u] == s1[u] && p[v] == s1[v])) {
+ final double t = ((p[u] - s0[u]) * dx + (p[v] - s0[v]) * dy) / (dx*dx + dy*dy);
+ if (t > 0 && t < 1) {
+ points.add(new SegmentTriangle(SegmentTriangle.Type.VERTEX, j, t));
+ }
+ }
+ }
+ /*
+ * The edges crossed between their corners, by the segment between its ends.
+ */
+ for (int i = 0; i < 3; i++) {
+ final int j = (i + 1) % 3;
+ if (opposite(sides[i], sides[j])) {
+ final double d0 = side(corners[i], corners[j], s0, u, v);
+ final double d1 = side(corners[i], corners[j], s1, u, v);
+ if (opposite(d0, d1)) {
+ points.add(new SegmentTriangle(SegmentTriangle.Type.EDGE, i, ratio(d0, d1)));
+ }
+ }
+ }
+ }
+ switch (points.size()) {
+ case 0: return List.of();
+ case 1: return List.of(points.get(0));
+ default: {
+ //the triangle is convex, so the part inside it goes from the first point to the last one.
+ points.sort(Comparator.comparingDouble(SegmentTriangle::ratio));
+ return List.of(points.get(0), points.get(points.size() - 1));
+ }
+ }
+ }
+
+ /**
+ * Computes the intersection of two triangles in three dimensions. A corner of a triangle within the tolerance
+ * of the plane of the other triangle is on that plane, as decided by {@link Maths#planeSide Maths.planeSide(…)}.
+ * If the three corners of a triangle are on the plane of the other one, the triangles are
+ * {@linkplain TriangleTriangle.Type#COPLANAR coplanar} and shall be intersected in that plane,
+ * for example with {@link #segmentInTriangle2D segmentInTriangle2D(…)}.
+ *
+ * <p>Otherwise the intersection is on the line where the two planes meet, and each of its ends is a point
+ * where an edge of a triangle touches the other triangle, as given by {@link #segmentTriangle segmentTriangle(…)}.
+ * An edge in the plane of the other triangle is intersected with it in that plane.</p>
+ *
+ * @param a0 first corner of the first triangle.
+ * @param a1 second corner of the first triangle.
+ * @param a2 third corner of the first triangle.
+ * @param b0 first corner of the second triangle.
+ * @param b1 second corner of the second triangle.
+ * @param b2 third corner of the second triangle.
+ * @param tolerance largest distance to the plane of a triangle for considering a point on it.
+ * @return the intersection of the two triangles.
+ */
+ public static TriangleTriangle triangleTriangle(double[] a0, double[] a1, double[] a2,
+ double[] b0, double[] b1, double[] b2, double tolerance)
+ {
+ final double[][] first = {a0, a1, a2}, second = {b0, b1, b2};
+ final int[] sidesOfSecond = sides(first, second, tolerance);
+ final int[] sidesOfFirst = sides(second, first, tolerance);
+ if (allOnSameSide(sidesOfSecond) || allOnSameSide(sidesOfFirst)) {
+ return TriangleTriangle.DISJOINT;
+ }
+ if (Arrays.equals(sidesOfSecond, new int[3]) || Arrays.equals(sidesOfFirst, new int[3])) {
+ return TriangleTriangle.COPLANAR;
+ }
+ final List<TriangleEdgeContact> contacts = new ArrayList<>();
+ addEdgeContacts(first, second, true, tolerance, contacts);
+ addEdgeContacts(second, first, false, tolerance, contacts);
+ if (contacts.isEmpty()) {
+ return TriangleTriangle.DISJOINT;
+ }
+ /*
+ * All points are on the line where the two planes meet.
+ * The ends of the intersection are the two extreme points along that line.
+ */
+ final double[] direction = Vectors.cross(normal(a0, a1, a2), normal(b0, b1, b2));
+ TriangleEdgeContact start = null, end = null;
+ double min = Double.POSITIVE_INFINITY, max = Double.NEGATIVE_INFINITY;
+ for (final TriangleEdgeContact contact : contacts) {
+ final double t = Vectors.dot(contact.position(), direction);
+ if (t < min) {min = t; start = contact;}
+ if (t > max) {max = t; end = contact;}
+ }
+ if (start == end || locations(start).equals(locations(end)) || Arrays.equals(start.position(), end.position())) {
+ return new TriangleTriangle(TriangleTriangle.Type.POINT, start, start);
+ }
+ return new TriangleTriangle(TriangleTriangle.Type.SEGMENT, start, end);
+ }
+
+ /**
+ * Returns the sides of the corners of a triangle relative to the plane of another triangle.
+ */
+ private static int[] sides(final double[][] plane, final double[][] triangle, final double tolerance) {
+ final int[] sides = new int[3];
+ for (int i = 0; i < 3; i++) {
+ sides[i] = Maths.planeSide(plane[0], plane[1], plane[2], triangle[i], tolerance);
+ }
+ return sides;
+ }
+
+ /**
+ * Returns whether the three sides are all strictly positive or all strictly negative.
+ */
+ private static boolean allOnSameSide(final int[] sides) {
+ return sameSide(sides[0], sides[1]) && sameSide(sides[1], sides[2]);
+ }
+
+ /**
+ * Adds the points where the edges of a triangle touch another triangle.
+ *
+ * @param edges the triangle whose edges are tested.
+ * @param other the triangle touched by the edges.
+ * @param edgeOfFirst whether the triangle of the edges is the first one given to {@code triangleTriangle(…)}.
+ * @param tolerance largest distance to the plane of a triangle for considering a point on it.
+ * @param contacts where to add the points.
+ */
+ private static void addEdgeContacts(final double[][] edges, final double[][] other, final boolean edgeOfFirst,
+ final double tolerance, final List<TriangleEdgeContact> contacts)
+ {
+ for (int i = 0; i < 3; i++) {
+ final double[] p = edges[i], q = edges[(i + 1) % 3];
+ final SegmentTriangle contact = segmentTriangle(p, q, other[0], other[1], other[2], tolerance);
+ switch (contact.type()) {
+ case NONE: break;
+ case COPLANAR: {
+ final int[] axes = Maths.projectionAxes(normal(other[0], other[1], other[2]));
+ for (final SegmentTriangle point : segmentInTriangle2D(p, q, other[0], other[1], other[2], axes)) {
+ contacts.add(new TriangleEdgeContact(edgeOfFirst, i, point, interpolate(p, q, point.ratio())));
+ }
+ break;
+ }
+ default: {
+ contacts.add(new TriangleEdgeContact(edgeOfFirst, i, contact, interpolate(p, q, contact.ratio())));
+ break;
+ }
+ }
+ }
+ }
+
+ /**
+ * Returns the three first coordinates of the point at the given ratio between two points,
+ * exactly one of them at ratio 0 or 1.
+ */
+ private static double[] interpolate(final double[] p, final double[] q, final double ratio) {
+ if (ratio == 0) return Arrays.copyOf(p, 3);
+ if (ratio == 1) return Arrays.copyOf(q, 3);
+ final double[] point = new double[3];
+ for (int k = 0; k < 3; k++) {
+ point[k] = p[k] + ratio * (q[k] - p[k]);
+ }
+ return point;
+ }
+
+ /**
+ * Where a point is on a triangle: inside, or on the edge or the corner of the given index.
+ */
+ private record Location(SegmentTriangle.Type type, int index) {
+ }
+
+ /**
+ * Returns where the point of the given contact is on the first triangle and on the second triangle.
+ * Two contacts found from different edges are the same point if they are at the same locations.
+ */
+ private static List<Location> locations(final TriangleEdgeContact contact) {
+ final double ratio = contact.contact().ratio();
+ final int edge = contact.edge();
+ final Location onEdge;
+ if (ratio == 0) onEdge = new Location(SegmentTriangle.Type.VERTEX, edge);
+ else if (ratio == 1) onEdge = new Location(SegmentTriangle.Type.VERTEX, (edge + 1) % 3);
+ else onEdge = new Location(SegmentTriangle.Type.EDGE, edge);
+ final Location onOther = new Location(contact.contact().type(), contact.contact().feature());
+ return contact.edgeOfFirst() ? List.of(onEdge, onOther) : List.of(onOther, onEdge);
+ }
+
/**
* Intersection operations on BSpline curves and surfaces.
*/
diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java
index 00ed3a038d..d5ace07d12 100644
--- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java
+++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/triangulate/EarClipping.java
@@ -30,7 +30,7 @@ import org.apache.sis.geometries.surface.Polygon;
import org.apache.sis.maths.Array;
import org.apache.sis.maths.Maths;
import org.apache.sis.maths.NDArrays;
-import org.apache.sis.maths.Vector;
+import org.apache.sis.maths.ReadOnly;
import org.apache.sis.util.ArraysExt;
@@ -82,15 +82,16 @@ public class EarClipping {
private DataPoints source;
/**
- * Positions of the polygon, with their coordinates in the order used by the triangulation.
- * The triangulation works on the two first coordinates, in the plane where the polygon has its largest area.
+ * Positions of the vertices, with their coordinates in the order used by the triangulation.
+ * The triangulation works on the two first coordinates. For a polygon, those positions are
+ * the ones of {@link #source} in the plane where the polygon has its largest area.
*/
- private Vector<?>[] working;
+ private ReadOnly.Vector<?>[] working;
/**
- * Indices in {@link #source} of the vertices of the ring joining the holes to the exterior ring, without
- * repeated consecutive positions and without the closing position. The ring has the orientation of the
- * exterior ring.
+ * Indices in {@link #working} of the vertices of the ring to triangulate, without repeated consecutive
+ * positions and without the closing position. For a polygon, this is the ring joining the holes to the
+ * exterior ring, which has the orientation of the exterior ring.
*/
private int[] ring;
@@ -112,7 +113,7 @@ public class EarClipping {
private Angle[] angles;
/**
- * Indices in {@link #source} of the three vertices of each triangle found so far.
+ * Indices in {@link #working} of the three vertices of each triangle found so far.
*/
private final List<int[]> triangles = new ArrayList<>();
@@ -127,13 +128,39 @@ public class EarClipping {
public Triangulation triangulateIndexed(Polygon geometry) throws OperationException {
source = geometry.getDataPoints();
final int[] axes = workingAxes(geometry.getExteriorRing().getDataPoints().getAttributeArray(DataPointsType.ATT_POSITION));
- working = new Vector<?>[source.size()];
+ working = new ReadOnly.Vector<?>[source.size()];
Arrays.setAll(working, (i) -> source.getPosition(i).swizzle(axes));
triangles.clear();
run(HoleToExterior.ring(geometry));
return result();
}
+ /**
+ * Triangulates a simple ring given by the indices of its vertices in an array of positions.
+ * The triangulation works on the two first coordinates of the positions: a ring in three dimensions
+ * shall be projected first, for example with {@link ReadOnly.Vector#swizzle(int...)} and the axes given by
+ * {@link Maths#projectionAxes(double[])}. Consecutive vertices at the same position are a single vertex.
+ *
+ * @param positions positions of the vertices, of which the ring may use only some.
+ * @param indices indices in {@code positions} of the vertices of the ring, in order. The ring can be closed,
+ * with its last vertex at the location of the first one, or not.
+ * @return indices in {@code positions} of the three corners of each triangle, one triangle after the other.
+ * The triangles have the orientation of the ring.
+ * @throws OperationException if the ring cannot be triangulated, for example because it intersects itself.
+ */
+ public int[] triangulate(final ReadOnly.Vector<?>[] positions, int[] indices) throws OperationException {
+ working = positions;
+ triangles.clear();
+ final int n = indices.length;
+ if (n != 0 && !positions[indices[0]].equals(positions[indices[n - 1]])) {
+ //the triangulation expects a closed ring.
+ indices = Arrays.copyOf(indices, n + 1);
+ indices[n] = indices[0];
+ }
+ run(indices);
+ return triangles.stream().flatMapToInt(Arrays::stream).toArray();
+ }
+
/**
* Triangulates the given polygon into an indexed mesh holding all the attributes of the polygon.
*
@@ -152,7 +179,7 @@ public class EarClipping {
/**
* Returns the order of the coordinates to use for triangulating a polygon with the given positions,
- * to be given to {@link Vector#swizzle(int...)}. The triangulation works on the two first coordinates
+ * to be given to {@link ReadOnly.Vector#swizzle(int...)}. The triangulation works on the two first coordinates
* of that order. With three or more dimensions, the three first axes are in the {@linkplain Maths#projectionAxes
* projection order} of the {@linkplain Maths#calculatePolygonNormal(Array) normal of the polygon}: the axis along
* which the normal is the largest goes after the two axes of the projection plane, so that a vertical polygon is
@@ -189,7 +216,7 @@ public class EarClipping {
* triangle with its two neighbors contains no other vertex. Cutting it removes the vertex from the
* ring, which stays a valid polygon, until only one triangle is left.
*
- * @param joined indices in {@link #source} of the ring positions, the last one at the location of the first one.
+ * @param joined indices in {@link #working} of the ring positions, the last one at the location of the first one.
* @throws OperationException if the vertices left without ear still enclose an area.
*/
private void run(final int[] joined) throws OperationException {
@@ -239,19 +266,19 @@ public class EarClipping {
* @param start a vertex which is not cut yet.
*/
private boolean isDegenerate(final int start) {
- final Vector<?> origin = position(start);
+ final ReadOnly.Vector<?> origin = position(start);
final double ox = origin.get(0), oy = origin.get(1);
double area = 0; // twice the signed area.
int vertex = start;
do {
- final Vector<?> a = position(vertex), b = position(next[vertex]);
+ final ReadOnly.Vector<?> a = position(vertex), b = position(next[vertex]);
area += (a.get(0) - ox) * (b.get(1) - oy) - (b.get(0) - ox) * (a.get(1) - oy);
vertex = next[vertex];
} while (vertex != start);
double xmin = Double.POSITIVE_INFINITY, xmax = Double.NEGATIVE_INFINITY;
double ymin = Double.POSITIVE_INFINITY, ymax = Double.NEGATIVE_INFINITY;
for (final int index : ring) {
- final Vector<?> p = working[index];
+ final ReadOnly.Vector<?> p = working[index];
xmin = Math.min(xmin, p.get(0)); xmax = Math.max(xmax, p.get(0));
ymin = Math.min(ymin, p.get(1)); ymax = Math.max(ymax, p.get(1));
}
@@ -262,7 +289,7 @@ public class EarClipping {
/**
* Returns the working position of the given vertex of the ring.
*/
- private Vector<?> position(final int vertex) {
+ private ReadOnly.Vector<?> position(final int vertex) {
return working[ring[vertex]];
}
@@ -283,11 +310,11 @@ public class EarClipping {
* so a copy located on a corner of the triangle does not prevent the ear.
*/
private boolean isEar(final int vertex) {
- final Vector<?> t1 = position(previous[vertex]);
- final Vector<?> t2 = position(vertex);
- final Vector<?> t3 = position(next[vertex]);
+ final ReadOnly.Vector<?> t1 = position(previous[vertex]);
+ final ReadOnly.Vector<?> t2 = position(vertex);
+ final ReadOnly.Vector<?> t3 = position(next[vertex]);
for (int i = next[next[vertex]]; i != previous[vertex]; i = next[i]) {
- final Vector<?> p = position(i);
+ final ReadOnly.Vector<?> p = position(i);
if (!isConvex(i) && !p.equals(t1) && !p.equals(t2) && !p.equals(t3) && Maths.inTriangle(t1, t2, t3, p)) {
return false;
}
diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Maths.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Maths.java
index 185505bc55..799a9bed0e 100644
--- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Maths.java
+++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Maths.java
@@ -167,6 +167,34 @@ public final class Maths {
return (normal[drop] < 0) ? new int[] {v, u, drop} : new int[] {u, v, drop};
}
+ /**
+ * Returns the side of a point relative to the plane of a triangle, with a tolerance. The point is on the plane
+ * if its distance to the plane is not greater than the tolerance. Otherwise its side is given by the robust
+ * {@link ShewchukPredicates#orient3d orient3d} predicate, so that the side is right even very close to the plane
+ * when the tolerance is zero. Only the three first coordinates are used.
+ *
+ * @param a first corner of the triangle
+ * @param b second corner of the triangle
+ * @param c third corner of the triangle
+ * @param p point to test
+ * @param tolerance largest distance to the plane for considering the point on it, zero for an exact test
+ * @return +1 if the point is on the side toward which the normal (b−a)×(c−a) points, −1 if it is on the
+ * other side, 0 if it is on the plane. A degenerate triangle has no plane: every point is on it.
+ */
+ public static int planeSide(double[] a, double[] b, double[] c, double[] p, double tolerance) {
+ final double ux = b[0] - a[0], uy = b[1] - a[1], uz = b[2] - a[2];
+ final double vx = c[0] - a[0], vy = c[1] - a[1], vz = c[2] - a[2];
+ final double nx = uy*vz - uz*vy, ny = uz*vx - ux*vz, nz = ux*vy - uy*vx;
+ //the distance to the plane multiplied by the length of the normal.
+ final double distance = nx * (p[0] - a[0]) + ny * (p[1] - a[1]) + nz * (p[2] - a[2]);
+ if (Math.abs(distance) <= tolerance * Math.sqrt(nx*nx + ny*ny + nz*nz)) {
+ return 0;
+ }
+ //orient3d is positive below the plane, on the side opposite to the normal.
+ return (int) -Math.signum(ShewchukPredicates.orient3d(a[0], a[1], a[2], b[0], b[1], b[2],
+ c[0], c[1], c[2], p[0], p[1], p[2]));
+ }
+
public static double[] subtract(float[] A, float[] B){
return new double[] {
A[0] - (double) B[0],
diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/Contains3DTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/Contains3DTest.java
new file mode 100644
index 0000000000..9e703703ab
--- /dev/null
+++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/Contains3DTest.java
@@ -0,0 +1,182 @@
+/*
+ * 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;
+
+import java.util.Arrays;
+import java.util.List;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.apache.sis.geometries.Geometries;
+import org.apache.sis.geometries.GeometryCollection;
+import org.apache.sis.geometries.GeometryFactory;
+import org.apache.sis.geometries.curve.LinearRing;
+import org.apache.sis.geometries.solid.Polyhedron;
+import org.apache.sis.geometries.solid.polyhedron.Icosahedron;
+import org.apache.sis.geometries.surface.Polygon;
+import org.apache.sis.maths.Array;
+import org.apache.sis.maths.Maths;
+import org.apache.sis.maths.NDArrays;
+import org.apache.sis.maths.Vectors;
+import org.apache.sis.referencing.CommonCRS;
+
+// Test dependencies
+import static org.junit.jupiter.api.Assertions.*;
+import org.junit.jupiter.api.Test;
+
+
+/**
+ * Tests the three-dimensional and planar tests of {@link Contains}: winding number of a surface,
+ * polyhedron containing a point, and ring containing a point in any plane.
+ *
+ * @author Johann Sorel (Geomatys)
+ */
+public class Contains3DTest {
+
+ private static final double EPS = 1E-9;
+
+ private static final CoordinateReferenceSystem CRS3D = Geometries.getUndefinedCRS(3);
+
+ /**
+ * Indices of the corners of the twelve triangles of a cube, oriented outward. The corner of index
+ * <var>x</var> + 2<var>y</var> + 4<var>z</var> is at the minimum (0) or maximum (1) of each axis.
+ */
+ private static final int[] CUBE = {
+ 0,2,3, 0,3,1, // z = 0
+ 4,5,7, 4,7,6, // z = 1
+ 0,1,5, 0,5,4, // y = 0
+ 2,6,7, 2,7,3, // y = 1
+ 0,4,6, 0,6,2, // x = 0
+ 1,3,7, 1,7,5 // x = 1
+ };
+
+ /**
+ * Returns the coordinates of the eight corners of the cube between the given minimum and maximum on all axes.
+ */
+ private static double[] corners(final double min, final double max) {
+ final double[] xyz = new double[8 * 3];
+ for (int i = 0; i < 8; i++) {
+ xyz[i*3 ] = ((i & 1) == 0) ? min : max;
+ xyz[i*3+1] = ((i & 2) == 0) ? min : max;
+ xyz[i*3+2] = ((i & 4) == 0) ? min : max;
+ }
+ return xyz;
+ }
+
+ /**
+ * Returns the triangles of the cube, reversed so that they are oriented inward, with indices shifted by the given offset.
+ */
+ private static int[] inward(final int offset) {
+ final int[] triangles = new int[CUBE.length];
+ for (int i = 0; i < CUBE.length; i += 3) {
+ triangles[i ] = CUBE[i ] + offset;
+ triangles[i+1] = CUBE[i+2] + offset;
+ triangles[i+2] = CUBE[i+1] + offset;
+ }
+ return triangles;
+ }
+
+ /**
+ * The winding number is 1 inside a closed surface oriented outward, −1 if oriented inward, 0 outside,
+ * and ½ on a face.
+ */
+ @Test
+ public void testWindingNumber() {
+ final double[] xyz = corners(0, 1);
+ assertEquals( 1, Contains.windingNumber(xyz, CUBE, new double[] {0.5, 0.5, 0.5}), EPS);
+ assertEquals( 1, Contains.windingNumber(xyz, CUBE, new double[] {0.9, 0.1, 0.3}), EPS);
+ assertEquals( 0, Contains.windingNumber(xyz, CUBE, new double[] {2.0, 0.5, 0.5}), EPS);
+ assertEquals( 0, Contains.windingNumber(xyz, CUBE, new double[] {-3, 7, 0.2}), EPS);
+ assertEquals(-1, Contains.windingNumber(xyz, inward(0), new double[] {0.5, 0.5, 0.5}), EPS);
+ assertEquals( 0.5, Contains.windingNumber(xyz, CUBE, new double[] {0.3, 0.6, 0.0}), EPS, "On a face.");
+ }
+
+ /**
+ * The winding numbers of an outer surface and of the inward surface of its void add,
+ * giving 0 inside the void.
+ */
+ @Test
+ public void testWindingNumberWithVoid() {
+ final double[] outer = corners(0, 4), inner = corners(1, 3);
+ final double[] xyz = Arrays.copyOf(outer, outer.length + inner.length);
+ System.arraycopy(inner, 0, xyz, outer.length, inner.length);
+ final int[] triangles = Arrays.copyOf(CUBE, CUBE.length * 2);
+ System.arraycopy(inward(8), 0, triangles, CUBE.length, CUBE.length);
+ assertEquals(1, Contains.windingNumber(xyz, triangles, new double[] {0.5, 0.5, 0.5}), EPS, "In the solid.");
+ assertEquals(0, Contains.windingNumber(xyz, triangles, new double[] {2, 2, 2 }), EPS, "In the void.");
+ assertEquals(0, Contains.windingNumber(xyz, triangles, new double[] {5, 2, 2 }), EPS, "Outside.");
+ }
+
+ /**
+ * A polyhedron contains the points inside its exterior shell and outside its interior shells.
+ * The interior shell is oriented as the exterior one, which shall not matter.
+ */
+ @Test
+ public void testPolyhedron() {
+ final GeometryFactory f = GeometryFactory.DEFAULT;
+ final Polyhedron cube = f.createPolyhedron(cube(0, 4), List.of());
+ assertTrue (Contains.contains(cube, f.createPoint(CRS3D, 2, 2, 2)));
+ assertTrue (Contains.contains(cube, f.createPoint(CRS3D, 0.1, 3.9, 0.2)));
+ assertFalse(Contains.contains(cube, f.createPoint(CRS3D, 5, 2, 2)));
+
+ final Polyhedron hollow = f.createPolyhedron(cube(0, 4), List.of(cube(1, 3)));
+ assertTrue (Contains.contains(hollow, f.createPoint(CRS3D, 0.5, 2, 2)), "In the solid.");
+ assertFalse(Contains.contains(hollow, f.createPoint(CRS3D, 2, 2, 2)), "In the void.");
+ assertFalse(Contains.contains(hollow, f.createPoint(CRS3D, 5, 2, 2)), "Outside.");
+
+ final Polyhedron icosahedron = new Icosahedron();
+ assertTrue (Contains.contains(icosahedron, f.createPoint(icosahedron.getCoordinateReferenceSystem(), 0.1, -0.2, 0.3)));
+ assertFalse(Contains.contains(icosahedron, f.createPoint(icosahedron.getCoordinateReferenceSystem(), 0.9, 0.9, 0)));
+
+ assertThrows(OperationException.class, () -> Contains.contains(cube, f.createPoint(CommonCRS.WGS84.geographic3D(), 2, 2, 2)));
+ }
+
+ /**
+ * A ring in a vertical plane contains a point when tested in the plane of its normal.
+ */
+ @Test
+ public void testRingInPlane() {
+ final Array ring = NDArrays.of(3, new double[] {0,2,0, 4,2,0, 4,2,4, 0,2,4, 0,2,0});
+ final int[] axes = Maths.projectionAxes(Maths.calculatePolygonNormal(ring));
+ assertTrue (Contains.contains(ring, Vectors.createDouble(3).set(new double[] {1, 2, 3}), axes));
+ assertFalse(Contains.contains(ring, Vectors.createDouble(3).set(new double[] {5, 2, 3}), axes));
+ assertFalse(Contains.contains(ring, Vectors.createDouble(3).set(new double[] {1, 2, 5}), axes));
+ }
+
+ /**
+ * Returns the six faces of the cube between the given minimum and maximum on all axes, oriented outward.
+ */
+ private static GeometryCollection<Polygon> cube(final double min, final double max) {
+ final double[] xyz = corners(min, max);
+ final GeometryFactory f = GeometryFactory.DEFAULT;
+ final Polygon[] faces = new Polygon[6];
+ for (int i = 0; i < 6; i++) {
+ //the two triangles of a face are (a,b,c) and (a,c,d), so the ring of the face is a, b, c, d.
+ final int a = CUBE[i*6], b = CUBE[i*6+1], c = CUBE[i*6+2], d = CUBE[i*6+5];
+ final double[] coordinates = new double[5 * 3];
+ int k = 0;
+ for (final int corner : new int[] {a, b, c, d, a}) {
+ System.arraycopy(xyz, corner * 3, coordinates, k, 3);
+ k += 3;
+ }
+ faces[i] = f.createPolygon(ring(coordinates), List.of());
+ }
+ return f.createMultiPolygon(CRS3D, faces);
+ }
+
+ private static LinearRing ring(final double... coordinates) {
+ return GeometryFactory.DEFAULT.createLinearRing(GeometryFactory.DEFAULT.createDataPoints(NDArrays.of(CRS3D, coordinates)));
+ }
+}
diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/Intersection3DTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/Intersection3DTest.java
new file mode 100644
index 0000000000..3c241dfac0
--- /dev/null
+++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/Intersection3DTest.java
@@ -0,0 +1,264 @@
+/*
+ * 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;
+
+import java.util.List;
+import org.apache.sis.geometries.operation.Intersection.SegmentSegment;
+import org.apache.sis.geometries.operation.Intersection.SegmentTriangle;
+import org.apache.sis.geometries.operation.Intersection.TriangleTriangle;
+import org.apache.sis.maths.Maths;
+import org.apache.sis.maths.Vectors;
+
+// Test dependencies
+import static org.junit.jupiter.api.Assertions.*;
+import org.junit.jupiter.api.Test;
+
+
+/**
+ * Tests the intersections of segments and triangles of {@link Intersection}.
+ *
+ * @author Johann Sorel (Geomatys)
+ */
+public class Intersection3DTest {
+
+ private static final double EPS = 1E-12;
+
+ /**
+ * The axes of the (x,y) plane.
+ */
+ private static final int[] XY = {0, 1};
+
+ /**
+ * The triangle used by most tests, in the (x,y) plane, counterclockwise seen from above.
+ */
+ private static final double[] A = {0, 0, 0}, B = {4, 0, 0}, C = {0, 4, 0};
+
+ private static double[] p(final double... coordinates) {
+ return coordinates;
+ }
+
+ /**
+ * Segments crossing, touching at an end, parallel, or apart.
+ */
+ @Test
+ public void testSegmentSegment() {
+ SegmentSegment r = Intersection.segmentSegment2D(p(0,0), p(2,2), p(0,2), p(2,0), XY);
+ assertEquals(SegmentSegment.Type.POINT, r.type());
+ assertEquals(0.5, r.startOnFirst(), EPS);
+ assertEquals(0.5, r.startOnSecond(), EPS);
+
+ r = Intersection.segmentSegment2D(p(0,0), p(1,0), p(1,0), p(1,1), XY);
+ assertEquals(SegmentSegment.Type.POINT, r.type(), "Touching at their ends.");
+ assertEquals(1, r.startOnFirst(), "The end of a segment shall be exact.");
+ assertEquals(0, r.startOnSecond(), "The end of a segment shall be exact.");
+
+ r = Intersection.segmentSegment2D(p(0,0), p(2,0), p(1,0), p(1,1), XY);
+ assertEquals(SegmentSegment.Type.POINT, r.type(), "An end on the other segment.");
+ assertEquals(0.5, r.startOnFirst(), EPS);
+ assertEquals(0, r.startOnSecond());
+
+ assertSame(SegmentSegment.NONE, Intersection.segmentSegment2D(p(0,0), p(1,0), p(0,1), p(1,1), XY), "Parallel.");
+ assertSame(SegmentSegment.NONE, Intersection.segmentSegment2D(p(0,0), p(1,0), p(2,-1), p(2,1), XY), "Lines crossing apart.");
+
+ //segments in the vertical plane y = 5, in its (x,z) axes.
+ r = Intersection.segmentSegment2D(p(0,5,0), p(2,5,2), p(0,5,2), p(2,5,0), new int[] {0, 2});
+ assertEquals(SegmentSegment.Type.POINT, r.type());
+ assertEquals(0.5, r.startOnFirst(), EPS);
+ }
+
+ /**
+ * Segments on the same line, overlapping, touching at an end, or apart.
+ */
+ @Test
+ public void testCollinearSegments() {
+ SegmentSegment r = Intersection.segmentSegment2D(p(0,0), p(4,0), p(3,0), p(1,0), XY);
+ assertEquals(SegmentSegment.Type.OVERLAP, r.type());
+ assertEquals(0.25, r.startOnFirst(), EPS);
+ assertEquals(0.75, r.endOnFirst(), EPS);
+ assertEquals(1, r.startOnSecond(), "The point at (1,0) is the end of the second segment.");
+ assertEquals(0, r.endOnSecond(), "The point at (3,0) is the start of the second segment.");
+
+ r = Intersection.segmentSegment2D(p(4,0), p(0,0), p(1,0), p(3,0), XY);
+ assertEquals(SegmentSegment.Type.OVERLAP, r.type(), "First segment in the opposite direction.");
+ assertEquals(0.25, r.startOnFirst(), EPS);
+ assertEquals(0.75, r.endOnFirst(), EPS);
+ assertEquals(1, r.startOnSecond());
+ assertEquals(0, r.endOnSecond());
+
+ r = Intersection.segmentSegment2D(p(0,0), p(1,1), p(1,1), p(2,2), XY);
+ assertEquals(SegmentSegment.Type.POINT, r.type(), "Touching at their ends.");
+ assertEquals(1, r.startOnFirst());
+ assertEquals(0, r.startOnSecond());
+
+ assertSame(SegmentSegment.NONE, Intersection.segmentSegment2D(p(0,0), p(1,0), p(2,0), p(3,0), XY));
+ }
+
+ /**
+ * A segment crossing the plane of the triangle inside it, on an edge, on a corner, or outside.
+ */
+ @Test
+ public void testSegmentTriangle() {
+ assertContact(SegmentTriangle.Type.INTERIOR, -1, 0.5, p(1,1,-1), p(1,1,1));
+ assertContact(SegmentTriangle.Type.EDGE, 0, 0.5, p(2,0,-1), p(2,0,1));
+ assertContact(SegmentTriangle.Type.EDGE, 1, 0.5, p(2,2,-1), p(2,2,1));
+ assertContact(SegmentTriangle.Type.EDGE, 2, 0.5, p(0,2,-1), p(0,2,1));
+ assertContact(SegmentTriangle.Type.VERTEX, 1, 0.5, p(4,0,-1), p(4,0,1));
+ assertContact(SegmentTriangle.Type.VERTEX, 0, 0.25, p(0,0,-1), p(0,0,3));
+ assertSame(SegmentTriangle.NONE, Intersection.segmentTriangle(p(5,5,-1), p(5,5,1), A, B, C, 0), "Outside.");
+ assertSame(SegmentTriangle.NONE, Intersection.segmentTriangle(p(1,1, 1), p(1,1,2), A, B, C, 0), "Above.");
+ }
+
+ /**
+ * A segment with an end on the plane of the triangle, or in that plane.
+ */
+ @Test
+ public void testSegmentEndOnTriangle() {
+ assertContact(SegmentTriangle.Type.INTERIOR, -1, 0, p(1,1,0), p(1,1,5));
+ assertContact(SegmentTriangle.Type.EDGE, 0, 1, p(2,0,3), p(2,0,0));
+ assertContact(SegmentTriangle.Type.VERTEX, 2, 0, p(0,4,0), p(1,1,1));
+ assertSame(SegmentTriangle.NONE, Intersection.segmentTriangle(p(5,5,0), p(5,5,1), A, B, C, 0), "End on the plane, outside.");
+ assertSame(SegmentTriangle.COPLANAR, Intersection.segmentTriangle(p(1,1,0), p(3,0,0), A, B, C, 0));
+ //an end near the plane is on it only within the tolerance.
+ assertSame(SegmentTriangle.NONE, Intersection.segmentTriangle(p(1,1,1E-12), p(1,1,5), A, B, C, 0));
+ final SegmentTriangle r = Intersection.segmentTriangle(p(1,1,1E-12), p(1,1,5), A, B, C, 1E-9);
+ assertEquals(SegmentTriangle.Type.INTERIOR, r.type());
+ assertEquals(0, r.ratio());
+ }
+
+ /**
+ * A segment crossing a triangle in an oblique plane, with coordinates which are not integers.
+ */
+ @Test
+ public void testSegmentObliqueTriangle() {
+ final double[] o = {0.1, 0.2, 0.3}, u = {1/3., 2/3., 2/3.}, v = {2/3., 1/3., -2/3.};
+ final double[] n = Vectors.cross(u, v);
+ final double[] a = o, b = Vectors.add(o, Vectors.scale(u, 4)), c = Vectors.add(o, Vectors.scale(v, 4));
+ final double[] inside = Vectors.add(o, Vectors.add(u, v));
+ final double[] outside = Vectors.add(o, Vectors.scale(Vectors.add(u, v), 3));
+ SegmentTriangle r = Intersection.segmentTriangle(Vectors.add(inside, Vectors.scale(n, -1)), Vectors.add(inside, Vectors.scale(n, 3)), a, b, c, 0);
+ assertEquals(SegmentTriangle.Type.INTERIOR, r.type());
+ assertEquals(0.25, r.ratio(), EPS);
+ r = Intersection.segmentTriangle(Vectors.add(outside, Vectors.scale(n, -1)), Vectors.add(outside, Vectors.scale(n, 1)), a, b, c, 0);
+ assertEquals(SegmentTriangle.Type.NONE, r.type());
+ }
+
+ /**
+ * Segments in the plane of the triangle: inside, crossing, through a corner, along an edge, or outside.
+ */
+ @Test
+ public void testSegmentInTriangle() {
+ List<SegmentTriangle> r = Intersection.segmentInTriangle2D(p(1,1,0), p(2,1,0), A, B, C, XY);
+ assertPoints(r, SegmentTriangle.Type.INTERIOR, -1, 0, SegmentTriangle.Type.INTERIOR, -1, 1);
+
+ r = Intersection.segmentInTriangle2D(p(-1,1,0), p(5,1,0), A, B, C, XY);
+ assertPoints(r, SegmentTriangle.Type.EDGE, 2, 1/6., SegmentTriangle.Type.EDGE, 1, 2/3.);
+
+ r = Intersection.segmentInTriangle2D(p(1,1,0), p(5,1,0), A, B, C, XY);
+ assertPoints(r, SegmentTriangle.Type.INTERIOR, -1, 0, SegmentTriangle.Type.EDGE, 1, 0.5);
+
+ r = Intersection.segmentInTriangle2D(p(-1,-1,0), p(2,2,0), A, B, C, XY);
+ assertPoints(r, SegmentTriangle.Type.VERTEX, 0, 1/3., SegmentTriangle.Type.EDGE, 1, 1);
+
+ r = Intersection.segmentInTriangle2D(p(-1,0,0), p(2,0,0), A, B, C, XY);
+ assertPoints(r, SegmentTriangle.Type.VERTEX, 0, 1/3., SegmentTriangle.Type.EDGE, 0, 1);
+
+ r = Intersection.segmentInTriangle2D(p(4,0,0), p(5,-1,0), A, B, C, XY);
+ assertEquals(1, r.size(), "Touching at a corner only.");
+ assertEquals(SegmentTriangle.Type.VERTEX, r.get(0).type());
+ assertEquals(1, r.get(0).feature());
+ assertEquals(0, r.get(0).ratio());
+
+ assertTrue(Intersection.segmentInTriangle2D(p(5,5,0), p(6,6,0), A, B, C, XY).isEmpty());
+
+ //a triangle in the vertical plane y = 2, projected on the plane of its normal.
+ final double[] a = {0,2,0}, b = {4,2,0}, c = {0,2,4};
+ final int[] axes = Maths.projectionAxes(Maths.calculateNormal(a, b, c));
+ r = Intersection.segmentInTriangle2D(p(1,2,1), p(1,2,5), a, b, c, axes);
+ assertPoints(r, SegmentTriangle.Type.INTERIOR, -1, 0, SegmentTriangle.Type.EDGE, 1, 0.5);
+ }
+
+ /**
+ * Triangles crossing, touching at a point or along an edge, coplanar or apart.
+ */
+ @Test
+ public void testTriangleTriangle() {
+ //a vertical triangle crossing the interior of the horizontal one.
+ TriangleTriangle r = Intersection.triangleTriangle(A, B, C, p(1,1,-1), p(1,1,1), p(2.5,1,0), 0);
+ assertEquals(TriangleTriangle.Type.SEGMENT, r.type());
+ assertEnds(r, p(1,1,0), p(2.5,1,0));
+
+ //a vertical triangle crossing an edge of the horizontal one.
+ r = Intersection.triangleTriangle(A, B, C, p(1,-1,-1), p(1,-1,1), p(1,2,0), 0);
+ assertEquals(TriangleTriangle.Type.SEGMENT, r.type());
+ assertEnds(r, p(1,0,0), p(1,2,0));
+
+ //a vertical triangle standing on an edge of the horizontal one.
+ r = Intersection.triangleTriangle(A, B, C, p(1,0,0), p(3,0,0), p(2,0,2), 0);
+ assertEquals(TriangleTriangle.Type.SEGMENT, r.type());
+ assertEnds(r, p(1,0,0), p(3,0,0));
+
+ //a triangle touching the horizontal one by a corner.
+ r = Intersection.triangleTriangle(A, B, C, p(1,1,0), p(2,1,1), p(1,2,1), 0);
+ assertEquals(TriangleTriangle.Type.POINT, r.type());
+ assertSame(r.start(), r.end());
+ assertArrayEquals(p(1,1,0), r.start().position());
+
+ assertSame(TriangleTriangle.COPLANAR, Intersection.triangleTriangle(A, B, C, p(1,1,0), p(5,1,0), p(1,5,0), 0));
+ assertSame(TriangleTriangle.DISJOINT, Intersection.triangleTriangle(A, B, C, p(10,10,1), p(11,10,1), p(10,11,2), 0), "Apart.");
+ assertSame(TriangleTriangle.DISJOINT, Intersection.triangleTriangle(A, B, C, p(5,5,-1), p(5,5,1), p(6,6,0), 0),
+ "Crossing the plane outside of the triangle.");
+ }
+
+ /**
+ * Verifies the result of the intersection of a segment with the triangle {@link #A}, {@link #B}, {@link #C}.
+ */
+ private static void assertContact(SegmentTriangle.Type type, int feature, double ratio, double[] s0, double[] s1) {
+ final SegmentTriangle r = Intersection.segmentTriangle(s0, s1, A, B, C, 0);
+ assertEquals(type, r.type());
+ assertEquals(feature, r.feature());
+ assertEquals(ratio, r.ratio(), EPS);
+ }
+
+ /**
+ * Verifies the two ends of the part of a segment inside a triangle.
+ */
+ private static void assertPoints(final List<SegmentTriangle> points,
+ final SegmentTriangle.Type type1, final int feature1, final double ratio1,
+ final SegmentTriangle.Type type2, final int feature2, final double ratio2)
+ {
+ assertEquals(2, points.size());
+ assertEquals(type1, points.get(0).type());
+ assertEquals(feature1, points.get(0).feature());
+ assertEquals(ratio1, points.get(0).ratio(), EPS);
+ assertEquals(type2, points.get(1).type());
+ assertEquals(feature2, points.get(1).feature());
+ assertEquals(ratio2, points.get(1).ratio(), EPS);
+ }
+
+ /**
+ * Verifies the ends of the intersection of two triangles, in any order.
+ */
+ private static void assertEnds(final TriangleTriangle r, final double[] p1, final double[] p2) {
+ final double[] start = r.start().position(), end = r.end().position();
+ if (Math.abs(start[0] - p1[0]) + Math.abs(start[1] - p1[1]) + Math.abs(start[2] - p1[2]) < EPS) {
+ assertArrayEquals(p2, end, EPS);
+ } else {
+ assertArrayEquals(p1, end, EPS);
+ assertArrayEquals(p2, start, EPS);
+ }
+ }
+}
diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/EarClippingTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/EarClippingTest.java
index a649a46445..0a1c535918 100644
--- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/EarClippingTest.java
+++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/operation/triangulate/EarClippingTest.java
@@ -28,6 +28,7 @@ import org.apache.sis.geometries.operation.OperationException;
import org.apache.sis.geometries.surface.Polygon;
import org.apache.sis.maths.Maths;
import org.apache.sis.maths.NDArrays;
+import org.apache.sis.maths.ReadOnly;
import org.apache.sis.maths.Vectors;
import org.apache.sis.referencing.CommonCRS;
@@ -127,6 +128,47 @@ public class EarClippingTest {
assertThrows(OperationException.class, () -> new EarClipping().triangulateIndexed(bowTie));
}
+ /**
+ * Triangulates a concave ring given by the indices of its vertices in an array of positions, some of them
+ * not in the ring. The ring gives n − 2 triangles covering its area, made of its vertices only and turning
+ * in its direction, whether it is closed or not.
+ */
+ @Test
+ public void testTriangulateRing() {
+ final ReadOnly.Vector<?>[] positions = {
+ vector(9, 9), // not in the ring.
+ vector(0, 0), vector(4, 0), vector(4, 4), vector(2, 1), vector(0, 4)
+ };
+ assertRingTriangles(positions, new int[] {1, 2, 3, 4, 5}, true);
+ assertRingTriangles(positions, new int[] {1, 2, 3, 4, 5, 1}, true);
+ assertRingTriangles(positions, new int[] {5, 4, 3, 2, 1}, false);
+
+ final ReadOnly.Vector<?>[] bowTie = {vector(0, 0), vector(2, 2), vector(2, 0), vector(0, 2)};
+ assertThrows(OperationException.class, () -> new EarClipping().triangulate(bowTie, new int[] {0, 1, 2, 3}));
+ }
+
+ /**
+ * Triangulates the given ring of five vertices, of area 10, and verifies the triangles.
+ */
+ private static void assertRingTriangles(final ReadOnly.Vector<?>[] positions, final int[] ring, final boolean counterclockwise) {
+ final int[] triangles = new EarClipping().triangulate(positions, ring);
+ assertEquals((5 - 2) * 3, triangles.length, "Expected n − 2 triangles.");
+ double area = 0;
+ for (int i = 0; i < triangles.length; i += 3) {
+ for (int k = 0; k < 3; k++) {
+ assertNotEquals(0, triangles[i+k], "The position is not a vertex of the ring.");
+ }
+ final double side = Maths.lineSide(positions[triangles[i]], positions[triangles[i+1]], positions[triangles[i+2]]);
+ assertTrue(counterclockwise ? side > 0 : side < 0, "The triangle shall turn in the direction of the ring.");
+ area += Math.abs(side) / 2;
+ }
+ assertEquals(10, area, 1E-12);
+ }
+
+ private static ReadOnly.Vector<?> vector(final double x, final double y) {
+ return Vectors.createDouble(2).set(new double[] {x, y});
+ }
+
/**
* Triangulates the given polygon and verifies that each triangle turns in the direction of the exterior ring,
* which is that its normal has the direction of the normal of the ring. Also verifies the area and the number
diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/MathsTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/MathsTest.java
index 140dde08a1..d997e4f5c0 100644
--- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/MathsTest.java
+++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/MathsTest.java
@@ -217,6 +217,34 @@ public class MathsTest {
assertThrows(IllegalArgumentException.class, () -> Maths.projectionAxes(new double[] {0, 1}));
}
+ /**
+ * Tests {@link Maths#planeSide(double[], double[], double[], double[], double)} on a horizontal triangle,
+ * with and without tolerance, and on a triangle in an oblique plane with coordinates which are not integers.
+ */
+ @Test
+ public void testPlaneSide() {
+ final double[] a = {0, 0, 0}, b = {1, 0, 0}, c = {0, 1, 0};
+ assertEquals( 1, Maths.planeSide(a, b, c, new double[] {0.3, 0.3, 1}, 0), "Above, toward the normal.");
+ assertEquals(-1, Maths.planeSide(a, b, c, new double[] {0.3, 0.3, -1}, 0), "Below.");
+ assertEquals(-1, Maths.planeSide(a, c, b, new double[] {0.3, 0.3, 1}, 0), "Reversed triangle.");
+ assertEquals( 0, Maths.planeSide(a, b, c, new double[] {5, 5, 0}, 0), "On the plane, outside of the triangle.");
+ assertEquals( 1, Maths.planeSide(a, b, c, new double[] {0.3, 0.3, 1E-12}, 0));
+ assertEquals( 0, Maths.planeSide(a, b, c, new double[] {0.3, 0.3, 1E-12}, 1E-9), "Within the tolerance.");
+ assertEquals( 1, Maths.planeSide(a, b, c, new double[] {0.3, 0.3, 2E-9}, 1E-9), "Beyond the tolerance.");
+ assertEquals( 0, Maths.planeSide(a, b, new double[] {2, 0, 0}, new double[] {0, 0, 1}, 0), "Degenerate triangle.");
+ /*
+ * A triangle made of two orthonormal vectors u and v of an oblique plane, whose normal is u × v.
+ */
+ final double[] o = {0.1, 0.2, 0.3}, u = {1/3., 2/3., 2/3.}, v = {2/3., 1/3., -2/3.};
+ final double[] n = Vectors.cross(u, v);
+ final double[] ou = Vectors.add(o, u), ov = Vectors.add(o, v);
+ final double[] inPlane = Vectors.add(o, Vectors.add(Vectors.scale(u, 0.4), Vectors.scale(v, 0.3)));
+ assertEquals( 1, Maths.planeSide(o, ou, ov, Vectors.add(inPlane, Vectors.scale(n, 1E-3)), 0));
+ assertEquals(-1, Maths.planeSide(o, ou, ov, Vectors.add(inPlane, Vectors.scale(n, -1E-3)), 0));
+ assertEquals( 0, Maths.planeSide(o, ou, ov, Vectors.add(inPlane, Vectors.scale(n, 1E-3)), 1E-2));
+ assertEquals( 0, Maths.planeSide(o, ou, ov, inPlane, 1E-12), "Rounding errors within the tolerance.");
+ }
+
/**
* Tests that {@link Maths#projectionAxes(double[])} keeps the orientation: a triangle whose normal has the
* direction given to the method is counterclockwise once projected, and the reversed triangle clockwise.