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 <[email protected]> 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.
