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 677d44e76396b2a48843649d825a81b3107eb5d8 Author: jsorel <[email protected]> AuthorDate: Thu Oct 15 09:44:58 2026 +0200 feat(Geometry): implement envelopes on Arc types --- .../sis/geometries/curve/ArcByCenterPoint.java | 7 + .../sis/geometries/internal/shared/Curves.java | 81 +++++++ .../internal/shared/DefaultArcByBulge.java | 36 ++- .../internal/shared/DefaultArcByCenterPoint.java | 70 +++++- .../internal/shared/DefaultCircularString.java | 18 +- .../geometries/internal/shared/DefaultConic.java | 249 ++++++++++++++++++++- .../sis/geometries/curve/ArcByBulgeTest.java | 57 ++++- .../sis/geometries/curve/ArcByCenterPointTest.java | 50 ++++- .../sis/geometries/curve/CircularStringTest.java | 45 +++- .../org/apache/sis/geometries/curve/ConicTest.java | 104 ++++++++- 10 files changed, 708 insertions(+), 9 deletions(-) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/ArcByCenterPoint.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/ArcByCenterPoint.java index 5781955ad5..4815ccbda0 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/ArcByCenterPoint.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/ArcByCenterPoint.java @@ -42,6 +42,13 @@ import org.apache.sis.maths.Array; * <p>Because the arc is defined by a bearing sweep around a centre, this parameterisation is * two-dimensional by nature; GML says as much.</p> * + * <p>Angles are measured from the direction of the first coordinate system axis, positive toward the + * direction of the second axis. With (<var>x</var>, <var>y</var>) axes this is the counter-clockwise + * angle from <var>x</var>; with (latitude, longitude) axes this is the clockwise azimuth from north. + * The arc goes from the start angle to the end angle, in the direction of increasing angles if the + * end angle is greater than the start angle and of decreasing angles otherwise. A difference of 360° + * or more between the two angles makes a full circle.</p> + * * @author Johann Sorel (Geomatys) * @see GML ArcByCenterPoint */ diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java index 47e6c59ea8..d68a9a3a61 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java @@ -21,6 +21,9 @@ import javax.measure.Quantity; import javax.measure.Unit; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.Point; +import org.apache.sis.geometry.GeneralEnvelope; +import org.apache.sis.maths.DataType; +import org.apache.sis.maths.ReadOnly; import org.apache.sis.maths.Vector; import org.apache.sis.measure.Units; @@ -101,6 +104,84 @@ public final class Curves { } } + /** + * Expands the given envelope for containing the circular arc going from {@code start} to {@code end} + * through {@code middle}. The envelope shall already contain the start and end points. This method + * adds each point where the circle reaches an extreme along an axis, if that point is on the arc. + * It works in any number of dimensions, the arc being in the plane of its three points. + * + * <p>If the three points are collinear, the arc is a straight segment and the envelope is unchanged. + * If the start and end points are identical, the arc is a full circle whose diameter goes from them + * to the middle point.</p> + * + * @param envelope the envelope to expand, containing the start and end points. + * @param start the start point. + * @param middle a point on the arc, distinct from its start and end points. + * @param end the end point. + * @throws IllegalStateException if the arc is a full circle in more than two dimensions, + * since its plane is then not determined. + */ + public static void addArc(final GeneralEnvelope envelope, final ReadOnly.Vector<?> start, + final ReadOnly.Vector<?> middle, final ReadOnly.Vector<?> end) + { + final int dimension = start.getDimension(); + final Vector<?> mid = start .copy(DataType.DOUBLE).lerp(end, 0.5); // Middle of the chord. + final Vector<?> u = end .copy(DataType.DOUBLE).subtract(start); // Unit vector along the chord. + final Vector<?> n = middle.copy(DataType.DOUBLE).subtract(mid); // Unit vector perpendicular to the chord, toward the arc. + final double chord = u.length(); + if (chord == 0) { + if (dimension != 2) { + throw new IllegalStateException("A circular arc with identical start and end points" + + " is a full circle, which plane is not determined in " + dimension + " dimensions."); + } + final double r = n.length() / 2; + final Vector<?> centre = start.copy(DataType.DOUBLE).lerp(middle, 0.5); + for (int i = 0; i < dimension; i++) { + extend(envelope, i, centre.get(i) - r); + extend(envelope, i, centre.get(i) + r); + } + return; + } + u.scale(1 / chord); + final double pu = n.dot(u); + n.perpendicular(u); + final double pn = n.length(); + if (!(pn > chord * 1E-12)) { + return; // Collinear points: a straight segment. + } + n.scale(1 / pn); + /* + * With the half chord c, the centre is at t·n from the middle of the chord, where t makes the + * centre equidistant from the start point and the middle point: c² + t² = pu² + (t − pn)². + * In the plane (u, n), the circle reaches its extremes along axis i in the directions + * ±(uᵢ·u + nᵢ·n), at distances ±r·hypot(uᵢ, nᵢ) from the centre along that axis. + * The arc is the part of the circle on the side of n from the chord. + */ + final double c = chord / 2; + final double t = (pu*pu + pn*pn - c*c) / (2*pn); + final double r = Math.hypot(c, t); + for (int i = 0; i < dimension; i++) { + final double h = Math.hypot(u.get(i), n.get(i)); + if (h == 0) { + continue; // The plane of the arc is perpendicular to this axis. + } + final double centre = mid.get(i) + t * n.get(i); + for (int sign = -1; sign <= 1; sign += 2) { + if (t + sign * r * n.get(i) / h >= 0) { + extend(envelope, i, centre + sign * r * h); + } + } + } + } + + /** + * Expands the range of the given envelope axis for containing the given value. + */ + private static void extend(final GeneralEnvelope envelope, final int axis, final double value) { + envelope.setRange(axis, Math.min(envelope.getMinimum(axis), value), + Math.max(envelope.getMaximum(axis), value)); + } + /** * Returns whether the given distance is expressed in percent of the curve length. * diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByBulge.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByBulge.java index 5f4cfb2040..53eb4a2a11 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByBulge.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByBulge.java @@ -19,6 +19,8 @@ package org.apache.sis.geometries.internal.shared; import org.opengis.geometry.Envelope; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.curve.ArcByBulge; +import org.apache.sis.geometry.GeneralEnvelope; +import org.apache.sis.maths.DataType; import org.apache.sis.maths.Vector; @@ -73,9 +75,41 @@ public non-sealed class DefaultArcByBulge extends AbstractGeometry implements Ar return normal; } + /** + * Returns the smallest envelope containing this arc. That envelope contains the two end points, + * plus each point where the circle reaches an extreme along an axis, if that point is on the side + * of the chord the arc bulges towards. This works in any number of dimensions, the arc being in + * the plane of the chord and the normal. + * + * @throws IllegalStateException if the two end points are identical, or if the normal is parallel to the chord. + */ @Override public Envelope getEnvelope() { - throw new UnsupportedOperationException("Not supported yet."); + final var envelope = new GeneralEnvelope(points.getEnvelope()); + if (bulge == 0) { + return envelope; // A straight segment. + } + final Vector<?> start = points.getPosition(0).copy(DataType.DOUBLE); + final Vector<?> end = points.getPosition(1).copy(DataType.DOUBLE); + final Vector<?> u = end.copy(DataType.DOUBLE).subtract(start); // Unit vector along the chord. + final Vector<?> n = normal.copy(DataType.DOUBLE); // Unit vector toward the side the arc bulges to. + final double chord = u.length(); + if (!(chord > 0)) { + throw new IllegalStateException("The two end points of an arc by bulge are identical, so they do not determine a circle."); + } + u.scale(1 / chord); + // Only the part of the normal perpendicular to the chord tells the side of the arc. + final double length = n.length(); + n.perpendicular(u); + final double perpendicular = n.length(); + if (!(perpendicular > length * 1E-12)) { + throw new IllegalStateException("The normal of an arc by bulge is parallel to its chord," + + " so it does not tell which side the arc bulges towards."); + } + // The point of the arc at the bulge distance from the middle of the chord. + final Vector<?> apex = start.copy(DataType.DOUBLE).lerp(end, 0.5).add(n.scale(bulge / perpendicular)); + Curves.addArc(envelope, start, apex, end); + return envelope; } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByCenterPoint.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByCenterPoint.java index a6c635be54..cc4007e538 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByCenterPoint.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArcByCenterPoint.java @@ -16,11 +16,14 @@ */ package org.apache.sis.geometries.internal.shared; +import javax.measure.IncommensurableException; import javax.measure.Unit; import org.opengis.geometry.Envelope; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.Point; import org.apache.sis.geometries.curve.ArcByCenterPoint; +import org.apache.sis.geometry.GeneralDirectPosition; +import org.apache.sis.geometry.GeneralEnvelope; /** @@ -93,9 +96,74 @@ public non-sealed class DefaultArcByCenterPoint extends AbstractGeometry impleme return endAngle; } + /** + * Returns the smallest envelope containing this arc. That envelope contains the start and end points, + * plus each point where the arc reaches an extreme along the first or second axis. Those extremes are + * at the multiples of 90° the arc sweeps over. + * + * @throws UnsupportedOperationException if the radius unit cannot be converted to the unit of an axis, + * for example a radius in metres around a centre in a geographic reference system. + */ @Override public Envelope getEnvelope() { - throw new UnsupportedOperationException("Not supported yet."); + final var envelope = new GeneralEnvelope(center.getEnvelope()); + final var position = new GeneralDirectPosition(envelope.getLowerCorner()); + final double cx = position.getCoordinate(0); + final double cy = position.getCoordinate(1); + final double rx = radius(0); + final double ry = radius(1); + // The sweep direction does not change the positions covered by the arc. + final double lo = Math.min(startAngle, endAngle); + final double hi = Math.max(startAngle, endAngle); + // The centre is not on the arc: the first two axes start from an end point, not from the centre. + double t = Math.toRadians(lo); + final double x = cx + rx * Math.cos(t); + final double y = cy + ry * Math.sin(t); + envelope.setRange(0, x, x); + envelope.setRange(1, y, y); + t = Math.toRadians(hi); + add(envelope, position, cx + rx * Math.cos(t), cy + ry * Math.sin(t)); + // Any sweep of 360° or more contains the four extremes, so the loop stops after one turn. + final double last = Math.min(hi, lo + 360); + for (double q = Math.ceil(lo / 90) * 90; q <= last; q += 90) { + switch (Math.floorMod((long) (q / 90), 4)) { + case 0: add(envelope, position, cx + rx, cy); break; + case 1: add(envelope, position, cx, cy + ry); break; + case 2: add(envelope, position, cx - rx, cy); break; + default: add(envelope, position, cx, cy - ry); break; + } + } + return envelope; + } + + /** + * Adds to the envelope the given position, with its first two coordinates replaced by the given ones. + * The other coordinates are those of the centre. + */ + private static void add(final GeneralEnvelope envelope, final GeneralDirectPosition position, + final double x, final double y) + { + position.setCoordinate(0, x); + position.setCoordinate(1, y); + envelope.add(position); + } + + /** + * Returns the radius expressed in the unit of the given axis. + * + * @throws UnsupportedOperationException if the radius unit cannot be converted to the axis unit. + */ + private double radius(final int axis) { + final Unit<?> axisUnit = getCoordinateReferenceSystem().getCoordinateSystem().getAxis(axis).getUnit(); + if (radiusUnit == null || axisUnit == null) { + return radius; + } + try { + return radiusUnit.getConverterToAny(axisUnit).convert(radius); + } catch (IncommensurableException e) { + throw new UnsupportedOperationException("The radius in " + radiusUnit + + " cannot be converted to the unit of axis " + axis + " (" + axisUnit + ").", e); + } } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCircularString.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCircularString.java index 064d9c6b41..7c5fb9dbf1 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCircularString.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCircularString.java @@ -19,6 +19,9 @@ package org.apache.sis.geometries.internal.shared; import org.opengis.geometry.Envelope; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.curve.CircularString; +import org.apache.sis.geometry.GeneralEnvelope; +import org.apache.sis.maths.DataType; +import org.apache.sis.maths.Vector; /** @@ -51,9 +54,22 @@ public non-sealed class DefaultCircularString extends AbstractGeometry implement return points; } + /** + * Returns the smallest envelope containing this circular string. That envelope contains the + * points of each arc, plus each point where an arc reaches an extreme along an axis. + * + * @throws IllegalStateException if an arc is a full circle in more than two dimensions. + */ @Override public Envelope getEnvelope() { - throw new UnsupportedOperationException("Not supported yet."); + final var envelope = new GeneralEnvelope(points.getEnvelope()); + Vector<?> start = (points.size() != 0) ? points.getPosition(0).copy(DataType.DOUBLE) : null; + for (int i = 2; i < points.size(); i += 2) { + final Vector<?> end = points.getPosition(i).copy(DataType.DOUBLE); + Curves.addArc(envelope, start, points.getPosition(i - 1).copy(DataType.DOUBLE), end); + start = end; + } + return envelope; } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConic.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConic.java index 7330f2aa54..b09ec81925 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConic.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConic.java @@ -16,13 +16,18 @@ */ package org.apache.sis.geometries.internal.shared; +import java.util.ArrayList; import java.util.Objects; import org.opengis.geometry.Envelope; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.curve.Conic; +import org.apache.sis.geometry.GeneralEnvelope; import org.apache.sis.maths.Array; +import org.apache.sis.maths.DataType; +import org.apache.sis.maths.Vector; +import org.apache.sis.maths.Vectors; /** @@ -99,9 +104,251 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { return points.isEmpty(); } + /** + * Returns the smallest envelope containing this conic. That envelope contains the data points, + * plus each point where an arc reaches an extreme along an axis. Each arc is computed in the plane + * of its five data points, with the coordinates handled as Cartesian ones: the exponential map at + * the control points is the identity. + * + * @throws IllegalStateException if the data points cannot be split in arcs of five points, + * or if the points of an arc are not coplanar, do not determine a single non-degenerate conic, + * are not in their order along the conic, or are joined through a point at infinity. + */ @Override public Envelope getEnvelope() { - throw new UnsupportedOperationException("Not supported yet."); + final var envelope = new GeneralEnvelope(points.getEnvelope()); + if (points.isEmpty()) { + return envelope; + } + final var positions = new ArrayList<Vector<?>>(); + for (int i = 0; i < points.size(); i++) { + positions.add(points.getPosition(i).copy(DataType.DOUBLE)); + } + if (cycle && positions.size() % 4 == 0) { + positions.add(positions.get(0)); // The first data point closes the curve. + } + final int n = positions.size(); + if (n < 5 || (n - 1) % 4 != 0) { + throw new IllegalStateException("A conic is made of arcs of five data points, consecutive arcs" + + " sharing a point, so it needs 4k+1 data points, but got " + n + '.'); + } + for (int i = 0; i + 4 < n; i += 4) { + addArc(envelope, positions.subList(i, i + 5).toArray(Vector<?>[]::new)); + } + return envelope; + } + + /** + * Expands the given envelope for containing the extremes along each axis of the conic arc + * going through the five given points. The envelope shall already contain those points. + */ + private static void addArc(final GeneralEnvelope envelope, final Vector<?>[] p) { + final int dimension = p[0].getDimension(); + /* + * Builds an orthonormal basis (e1, e2) of the plane of the points, with the centroid as origin. + * The coordinates in that plane are divided by the distance of the farthest point from the origin, + * for having values close to 1 in the computation of the conic coefficients. + */ + final Vector<?> origin = Vectors.createDouble(dimension); + for (final Vector<?> q : p) { + origin.add(q); + } + origin.scale(1.0 / p.length); + Vector<?> e1 = null, e2 = null; + double scale = 0, width = 0; + for (final Vector<?> q : p) { + final Vector<?> v = q.copy(DataType.DOUBLE).subtract(origin); + final double length = v.length(); + if (length > scale) { + scale = length; + e1 = v; + } + } + if (scale == 0) { + return; // All points are identical. + } + e1.scale(1 / scale); + for (final Vector<?> q : p) { + final Vector<?> v = q.copy(DataType.DOUBLE).subtract(origin).perpendicular(e1); + final double length = v.length(); + if (length > width) { + width = length; + e2 = v; + } + } + if (!(width > scale * 1E-10)) { + return; // Collinear points: a straight segment. + } + e2.scale(1 / width); + final double[] a = new double[p.length]; // Coordinates of the points in the plane. + final double[] b = new double[p.length]; + for (int k = 0; k < p.length; k++) { + final Vector<?> v = p[k].copy(DataType.DOUBLE).subtract(origin); + a[k] = v.dot(e1) / scale; + b[k] = v.dot(e2) / scale; + if (v.perpendicular(e1).perpendicular(e2).length() > scale * 1E-10) { + throw new IllegalStateException("The five data points of a conic arc are not coplanar."); + } + } + /* + * Coefficients of the conic A·a² + B·ab + C·b² + D·a + E·b + F = 0 through the five points, + * as the cofactors of the 5×6 matrix whose rows are (a², ab, b², a, b, 1) for each point. + */ + final Vector<?> coef = Vectors.createDouble(6); + for (int j = 0; j < 6; j++) { + final double[][] minor = new double[5][]; + for (int k = 0; k < 5; k++) { + final double[] row = {a[k]*a[k], a[k]*b[k], b[k]*b[k], a[k], b[k], 1}; + final double[] m = new double[5]; + for (int c = 0, d = 0; c < 6; c++) { + if (c != j) m[d++] = row[c]; + } + minor[k] = m; + } + coef.set(j, ((j & 1) == 0 ? 1 : -1) * determinant(minor)); + } + if (!(coef.length() > 1E-12)) { + throw new IllegalStateException("The five data points of a conic arc do not determine a single conic."); + } + coef.normalize(); + final double A = coef.get(0), B = coef.get(1), C = coef.get(2), D = coef.get(3), E = coef.get(4), F = coef.get(5); + if (Math.abs(determinant(new double[][] {{A, B/2, D/2}, {B/2, C, E/2}, {D/2, E/2, F}})) < 1E-12) { + throw new IllegalStateException("The five data points of a conic arc are on a degenerated conic," + + " made of straight lines, because three of them are collinear."); + } + /* + * The lines through the start point p₀ map the conic to the directions modulo π, continuously and + * one-to-one, p₀ itself being mapped to the direction of the tangent at p₀. The arc is therefore the + * interval of directions going from that tangent to the direction of p₄ through the other points. + */ + final double a0 = a[0], b0 = b[0]; + final double tangent = direction(-(B*a0 + 2*C*b0 + E), 2*A*a0 + B*b0 + D); + final double[] r = new double[5]; + for (int k = 1; k < 5; k++) { + r[k] = relative(direction(a[k] - a0, b[k] - b0), tangent); + } + final boolean increasing = r[2] < r[4]; + for (int k = 1; k < 4; k++) { + if (increasing ? !(r[k] < r[k+1]) : !(r[k] > r[k+1])) { + throw new IllegalStateException("The five data points of a conic arc are not in their order along the conic."); + } + } + final double end = r[4]; + /* + * Points at infinity, in the directions where A·cos²θ + B·cosθ·sinθ + C·sin²θ = 0, + * which is (A+C)/2 + R·cos(2θ − δ) = 0 with R and δ computed below. + */ + final double R = Math.hypot((A - C) / 2, B / 2); + if (R > 0) { + final double cos = -(A + C) / (2 * R); + if (Math.abs(cos) <= 1 + 1E-9) { + final double delta = Math.atan2(B / 2, (A - C) / 2); + final double acos = Math.acos(Math.max(-1, Math.min(1, cos))); + for (final double twice : new double[] {delta + acos, delta - acos}) { + final double rx = relative(direction(Math.cos(twice / 2), Math.sin(twice / 2)), tangent); + if (increasing ? rx <= end : rx >= end) { + throw new IllegalStateException("A conic arc goes through a point at infinity," + + " from one branch of a hyperbola to the other."); + } + } + } + } + /* + * Along axis i, the coordinate is a linear function α·a + β·b of the coordinates in the plane. + * Its extremes on the conic are where the gradient (Qa, Qb) is parallel to (α, β), which is the + * line β·Qa − α·Qb = 0, written g₁·a + g₂·b + g₃ = 0. Intersecting that line with the conic gives + * up to two points, which are kept if they are on the arc. + */ + for (int i = 0; i < dimension; i++) { + final double alpha = e1.get(i), beta = e2.get(i); + final double g1 = 2*A*beta - B*alpha; + final double g2 = B*beta - 2*C*alpha; + final double g3 = D*beta - E*alpha; + final double gg = g1*g1 + g2*g2; + if (!(gg > 0)) { + continue; + } + final double pa = -g3 * g1 / gg, pb = -g3 * g2 / gg; // A point on the line. + final double da = -g2, db = g1; // Direction of the line. + final double q2 = A*da*da + B*da*db + C*db*db; + final double q1 = (2*A*pa + B*pb + D) * da + (B*pa + 2*C*pb + E) * db; + final double q0 = A*pa*pa + B*pa*pb + C*pb*pb + D*pa + E*pb + F; + for (final double s : roots(q2, q1, q0)) { + final double xa = pa + s * da, xb = pb + s * db; + if (!Double.isFinite(xa) || !Double.isFinite(xb) || Math.hypot(xa - a0, xb - b0) < 1E-12) { + continue; + } + final double rx = relative(direction(xa - a0, xb - b0), tangent); + if (increasing ? rx <= end : rx >= end) { + final Vector<?> position = origin.copy(DataType.DOUBLE) + .add(e1.copy(DataType.DOUBLE).scale(scale * xa)) + .add(e2.copy(DataType.DOUBLE).scale(scale * xb)); + envelope.add(Vectors.asDirectPostion(position)); + } + } + } + } + + /** + * Returns the real roots of q₂·s² + q₁·s + q₀ = 0. + */ + private static double[] roots(final double q2, final double q1, final double q0) { + if (q2 == 0) { + return (q1 != 0) ? new double[] {-q0 / q1} : new double[0]; + } + final double discriminant = q1*q1 - 4*q2*q0; + if (discriminant < 0) { + return new double[0]; + } + final double q = -(q1 + Math.copySign(Math.sqrt(discriminant), q1)) / 2; + return (q != 0) ? new double[] {q / q2, q0 / q} : new double[] {0}; + } + + /** + * Returns the angle of the given direction, modulo π, in the [0 … π) range. + */ + private static double direction(final double x, final double y) { + double angle = Math.atan2(y, x); + if (angle < 0) angle += Math.PI; + if (angle >= Math.PI) angle -= Math.PI; + return angle; + } + + /** + * Returns the given direction relative to the reference one, in the [0 … π) range. + */ + private static double relative(final double angle, final double reference) { + final double r = angle - reference; + return (r < 0) ? r + Math.PI : r; + } + + /** + * Returns the determinant of the given square matrix, which is overwritten. + */ + private static double determinant(final double[][] m) { + final int n = m.length; + double det = 1; + for (int c = 0; c < n; c++) { + int pivot = c; + for (int r = c + 1; r < n; r++) { + if (Math.abs(m[r][c]) > Math.abs(m[pivot][c])) pivot = r; + } + if (m[pivot][c] == 0) { + return 0; + } + if (pivot != c) { + final double[] t = m[pivot]; m[pivot] = m[c]; m[c] = t; + det = -det; + } + det *= m[c][c]; + for (int r = c + 1; r < n; r++) { + final double f = m[r][c] / m[c][c]; + for (int k = c; k < n; k++) { + m[r][k] -= f * m[c][k]; + } + } + } + return det; } } diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByBulgeTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByBulgeTest.java index 64e6b83bb7..afc8878edd 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByBulgeTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByBulgeTest.java @@ -16,12 +16,14 @@ */ package org.apache.sis.geometries.curve; +import org.opengis.geometry.Envelope; import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.maths.Vector; import org.apache.sis.maths.Vectors; // Test dependencies import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertThrows; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -277,9 +279,62 @@ public class ArcByBulgeTest extends GeometryTest { public void testGetGeometryType2() { } + /** + * Creates a two-dimensional arc and returns its envelope. + */ + private static Envelope envelope(final double x0, final double y0, final double x1, final double y1, + final double bulge, final double nx, final double ny) + { + final Vector<?> normal = Vectors.createDouble(2); + normal.set(new double[] {nx, ny}); + return GeometryFactory.DEFAULT.createArcByBulge( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, x0, y0, x1, y1), bulge, normal).getEnvelope(); + } + + private static void assertEnvelope(final Envelope envelope, final double... bounds) { + final int dimension = bounds.length / 2; + assertEquals(dimension, envelope.getDimension()); + for (int i = 0; i < dimension; i++) { + assertEquals(bounds[i], envelope.getMinimum(i), TOLERANCE, "minimum"); + assertEquals(bounds[i + dimension], envelope.getMaximum(i), TOLERANCE, "maximum"); + } + } + @Test - @Disabled("Not implemented yet.") public void testGetEnvelope() { + // Minor arc: the box of the end points, extended to the top of the arc. + assertEnvelope(envelope(0, 0, 10, 0, 2, 0, 1), 0, 0, 10, 2); + assertEquals(CRS_2D, envelope(0, 0, 10, 0, 2, 0, 1).getCoordinateReferenceSystem()); + // A negative bulge bends the arc away from the normal. + assertEnvelope(envelope(0, 0, 10, 0, -2, 0, 1), 0, -2, 10, 0); + // The normal is neither normalized nor perpendicular to the chord: only its perpendicular part counts. + assertEnvelope(envelope(0, 0, 10, 0, 2, 3, 6), 0, 0, 10, 2); + // Half circle. + assertEnvelope(envelope(0, 0, 10, 0, 5, 0, 1), 0, 0, 10, 5); + // Major arc of radius 5.5625 centred on (5, 2.4375): it goes beyond the end points on both sides. + assertEnvelope(envelope(0, 0, 10, 0, 8, 0, 1), -0.5625, 0, 10.5625, 8); + // Null bulge: a straight segment. + assertEnvelope(envelope(0, 0, 10, 0, 0, 0, 1), 0, 0, 10, 0); + // Half circle on a diagonal chord, centred on (5, 5). + final double r = 5 * Math.sqrt(2); + assertEnvelope(envelope(0, 0, 10, 10, r, -1, 1), 5 - r, 0, 10, 5 + r); + } + + @Test + public void testGetEnvelope3D() { + final Vector<?> normal = Vectors.createDouble(3); + normal.set(new double[] {0, 0, 1}); + final ArcByBulge arc = GeometryFactory.DEFAULT.createArcByBulge( + GeometryFactory.DEFAULT.createDataPoints(CRS_3D, 0,0,0, 10,0,0), 2.0, normal); + assertEnvelope(arc.getEnvelope(), 0, 0, 0, 10, 0, 2); + } + + @Test + public void testGetEnvelopeUndefinedArc() { + // A normal parallel to the chord does not tell which side the arc bulges towards. + assertThrows(IllegalStateException.class, () -> envelope(0, 0, 10, 0, 2, 1, 0)); + // Two identical end points do not determine a circle. + assertThrows(IllegalStateException.class, () -> envelope(3, 3, 3, 3, 2, 0, 1)); } @Test diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByCenterPointTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByCenterPointTest.java index 732e9c6334..4c826c0b5a 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByCenterPointTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ArcByCenterPointTest.java @@ -16,11 +16,15 @@ */ package org.apache.sis.geometries.curve; +import javax.measure.Unit; +import org.opengis.geometry.Envelope; import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.Point; +import org.apache.sis.measure.Units; // Test dependencies import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertThrows; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -289,9 +293,53 @@ public class ArcByCenterPointTest extends GeometryTest { public void testGetGeometryType2() { } + /** + * Creates an arc around (10, 20) and returns its envelope. + */ + private static Envelope envelope(final double radius, final Unit<?> unit, final double startAngle, final double endAngle) { + final Point center = GeometryFactory.DEFAULT.createPoint(CRS_2D, 10, 20); + return GeometryFactory.DEFAULT.createArcByCenterPoint(center, radius, unit, startAngle, endAngle).getEnvelope(); + } + + private static void assertEnvelope(final Envelope envelope, final double xmin, final double ymin, + final double xmax, final double ymax) + { + assertEquals(CRS_2D, envelope.getCoordinateReferenceSystem()); + assertEquals(xmin, envelope.getMinimum(0), TOLERANCE, "xmin"); + assertEquals(ymin, envelope.getMinimum(1), TOLERANCE, "ymin"); + assertEquals(xmax, envelope.getMaximum(0), TOLERANCE, "xmax"); + assertEquals(ymax, envelope.getMaximum(1), TOLERANCE, "ymax"); + } + @Test - @Disabled("Not implemented yet.") public void testGetEnvelope() { + final double d = 5 * Math.sqrt(0.5); + // Quarter arc between two extremes: the envelope is the box of its end points. + assertEnvelope(envelope(5, null, 0, 90), 10, 20, 15, 25); + // Arc crossing the extreme in the direction of the second axis. + assertEnvelope(envelope(5, null, 45, 135), 10 - d, 20 + d, 10 + d, 25); + // Arc crossing the extreme in the direction of the first axis, at angle 0. + assertEnvelope(envelope(5, null, -45, 45), 10 + d, 20 - d, 15, 20 + d); + // Same arc, with angles beyond 360°. + assertEnvelope(envelope(5, null, 315, 405), 10 + d, 20 - d, 15, 20 + d); + // Clockwise sweep over the same quarter as the first case. + assertEnvelope(envelope(5, null, 90, 0), 10, 20, 15, 25); + // Clockwise sweep from 90° to -90° through 0°, not through 180°. + assertEnvelope(envelope(5, null, 90, -90), 10, 15, 15, 25); + // Sweep of more than half a turn, crossing three extremes. + assertEnvelope(envelope(5, null, 45, 315), 5, 15, 10 + d, 25); + // Full circle. + assertEnvelope(envelope(5, null, 0, 360), 5, 15, 15, 25); + // Degenerated arc: a single position. + assertEnvelope(envelope(5, null, 90, 90), 10, 25, 10, 25); + } + + @Test + public void testGetEnvelopeWithRadiusUnit() { + // 300 arc-minutes are 5 degrees, the unit of the axes of `CRS_2D`. + assertEnvelope(envelope(300, Units.ARC_MINUTE, 0, 360), 5, 15, 15, 25); + // A radius in metres cannot be converted to degrees without a geodesic computation. + assertThrows(UnsupportedOperationException.class, () -> envelope(5, Units.METRE, 0, 90)); } @Test diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/CircularStringTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/CircularStringTest.java index 70b5e97317..62407527cd 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/CircularStringTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/CircularStringTest.java @@ -16,6 +16,7 @@ */ package org.apache.sis.geometries.curve; +import org.opengis.geometry.Envelope; import org.apache.sis.geometries.GeometryFactory; // Test dependencies @@ -268,9 +269,51 @@ public class CircularStringTest extends GeometryTest { public void testGetGeometryType2() { } + /** + * Creates a two-dimensional circular string and returns its envelope. + */ + private static Envelope envelope(final double... coordinates) { + return GeometryFactory.DEFAULT.createCircularString( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, coordinates)).getEnvelope(); + } + + private static void assertEnvelope(final Envelope envelope, final double... bounds) { + final int dimension = bounds.length / 2; + assertEquals(dimension, envelope.getDimension()); + for (int i = 0; i < dimension; i++) { + assertEquals(bounds[i], envelope.getMinimum(i), TOLERANCE, "minimum"); + assertEquals(bounds[i + dimension], envelope.getMaximum(i), TOLERANCE, "maximum"); + } + } + @Test - @Disabled("Not implemented yet.") public void testGetEnvelope() { + // Half circle centred on (5, 0). + assertEnvelope(envelope(0,0, 5,5, 10,0), 0, 0, 10, 5); + assertEquals(CRS_2D, envelope(0,0, 5,5, 10,0).getCoordinateReferenceSystem()); + // Minor arc: the box of the points, the middle point being the top of the arc. + assertEnvelope(envelope(0,0, 5,2, 10,0), 0, 0, 10, 2); + // Major arc of radius 5.5625 centred on (5, 2.4375): it goes beyond the end points on both sides. + assertEnvelope(envelope(0,0, 5,8, 10,0), -0.5625, 0, 10.5625, 8); + // Arc of the circle of radius 5 centred on (0, 0), from 0° to about 233° counter-clockwise. + // The middle point is not where the arc is the furthest from the chord. + assertEnvelope(envelope(5,0, 0,5, -3,-4), -5, -4, 5, 5); + // Same arc in the other direction. + assertEnvelope(envelope(-3,-4, 0,5, 5,0), -5, -4, 5, 5); + // Two arcs making a full circle. + assertEnvelope(envelope(0,0, 5,5, 10,0, 5,-5, 0,0), 0, -5, 10, 5); + // Collinear points: a straight segment. + assertEnvelope(envelope(0,0, 5,0, 10,0), 0, 0, 10, 0); + // Identical start and end points: a full circle whose diameter goes from them to the middle point. + assertEnvelope(envelope(0,0, 10,0, 0,0), 0, -5, 10, 5); + } + + @Test + public void testGetEnvelope3D() { + // Half circle in the plane of the first and third axes. + final CircularString arc = GeometryFactory.DEFAULT.createCircularString( + GeometryFactory.DEFAULT.createDataPoints(CRS_3D, 0,0,0, 5,0,5, 10,0,0)); + assertEnvelope(arc.getEnvelope(), 0, 0, 0, 10, 0, 5); } @Test diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ConicTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ConicTest.java index 451a2311a2..b11b48a897 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ConicTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/ConicTest.java @@ -16,10 +16,12 @@ */ package org.apache.sis.geometries.curve; +import org.opengis.geometry.Envelope; import org.apache.sis.geometries.GeometryFactory; // Test dependencies import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertThrows; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -259,9 +261,107 @@ public class ConicTest extends GeometryTest { public void testGetGeometryType2() { } + /** + * Creates a two-dimensional conic and returns its envelope. + */ + private static Envelope envelope(final boolean cycle, final double... coordinates) { + return GeometryFactory.DEFAULT.createConic( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, coordinates), null, cycle).getEnvelope(); + } + + /** + * Returns the coordinates of the points of the ellipse x²/4 + y² = 1 at the given parameters. + */ + private static double[] ellipse(final double... t) { + final double[] coordinates = new double[t.length * 2]; + for (int i = 0; i < t.length; i++) { + coordinates[2*i] = 2 * Math.cos(t[i]); + coordinates[2*i + 1] = Math.sin(t[i]); + } + return coordinates; + } + + private static void assertEnvelope(final Envelope envelope, final double... bounds) { + final int dimension = bounds.length / 2; + assertEquals(dimension, envelope.getDimension()); + for (int i = 0; i < dimension; i++) { + assertEquals(bounds[i], envelope.getMinimum(i), TOLERANCE, "minimum"); + assertEquals(bounds[i + dimension], envelope.getMaximum(i), TOLERANCE, "maximum"); + } + } + @Test - @Disabled("Not implemented yet.") - public void testGetEnvelope() { + public void testGetEnvelopeEllipse() { + // Arc of the ellipse crossing its extremes along x at t = 0 and along y at t = π/2. + final Envelope envelope = envelope(false, ellipse(-0.3, 0.4, 1.2, 2.0, 2.6)); + assertEquals(CRS_2D, envelope.getCoordinateReferenceSystem()); + assertEnvelope(envelope, 2 * Math.cos(2.6), Math.sin(-0.3), 2, 1); + // Same arc in the other direction. + assertEnvelope(envelope(false, ellipse(2.6, 2.0, 1.2, 0.4, -0.3)), 2 * Math.cos(2.6), Math.sin(-0.3), 2, 1); + // Arc of the other side of the ellipse, from the same end points: it crosses the extremes at π and 3π/2. + assertEnvelope(envelope(false, ellipse(2.6, 3.4, 4.2, 5.0, 2*Math.PI - 0.3)), -2, -1, 2 * Math.cos(0.3), Math.sin(2.6)); + // Full ellipse made of two arcs, the first point closing the second arc. + assertEnvelope(envelope(true, ellipse(0.1, 0.9, 1.7, 2.5, 3.3, 4.1, 4.9, 5.7)), -2, -1, 2, 1); + } + + @Test + public void testGetEnvelopeHyperbola() { + // Right branch of x² − y² = 1, crossing its vertex at (1, 0). + final double[] u = {-1, -0.4, 0.2, 0.7, 1.2}; + final double[] coordinates = new double[10]; + for (int i = 0; i < u.length; i++) { + coordinates[2*i] = Math.cosh(u[i]); + coordinates[2*i + 1] = Math.sinh(u[i]); + } + assertEnvelope(envelope(false, coordinates), 1, Math.sinh(-1), Math.cosh(1.2), Math.sinh(1.2)); + /* + * Arc going from the right branch to the left branch: it passes through the point at infinity + * in the direction of the asymptote y = x. + */ + assertThrows(IllegalStateException.class, () -> envelope(false, + Math.cosh(0.5), Math.sinh(0.5), Math.cosh(1.5), Math.sinh(1.5), + -Math.cosh(1.5), -Math.sinh(1.5), -Math.cosh(0.5), -Math.sinh(0.5), -1, 0)); + } + + @Test + public void testGetEnvelopeParabola() { + // Parabola y = x², crossing its vertex at (0, 0). + assertEnvelope(envelope(false, -1,1, -0.2,0.04, 0.5,0.25, 1,1, 2,4), -1, 0, 2, 4); + } + + @Test + public void testGetEnvelopeCollinear() { + assertEnvelope(envelope(false, 0,0, 1,1, 2,2, 3,3, 4,4), 0, 0, 4, 4); + } + + @Test + public void testGetEnvelope3D() { + // The ellipse of `testGetEnvelopeEllipse()`, translated in the horizontal plane at height 5. + final double[] t = {-0.3, 0.4, 1.2, 2.0, 2.6}; + final double[] coordinates = new double[15]; + for (int i = 0; i < t.length; i++) { + coordinates[3*i] = 10 + 2 * Math.cos(t[i]); + coordinates[3*i + 1] = 20 + Math.sin(t[i]); + coordinates[3*i + 2] = 5; + } + final Conic conic = GeometryFactory.DEFAULT.createConic( + GeometryFactory.DEFAULT.createDataPoints(CRS_3D, coordinates), null, false); + assertEnvelope(conic.getEnvelope(), 10 + 2 * Math.cos(2.6), 20 + Math.sin(-0.3), 5, 12, 21, 5); + // The same points with one of them out of the plane. + coordinates[14] = 6; + final Conic skew = GeometryFactory.DEFAULT.createConic( + GeometryFactory.DEFAULT.createDataPoints(CRS_3D, coordinates), null, false); + assertThrows(IllegalStateException.class, () -> skew.getEnvelope()); + } + + @Test + public void testGetEnvelopeInvalid() { + // Six points are not a whole number of arcs. + assertThrows(IllegalStateException.class, () -> envelope(false, ellipse(0, 0.5, 1, 1.5, 2, 2.5))); + // Points not in their order along the ellipse. + assertThrows(IllegalStateException.class, () -> envelope(false, ellipse(0, 1.2, 0.4, 2.0, 2.6))); + // Three collinear points: the conic degenerates into a pair of straight lines. + assertThrows(IllegalStateException.class, () -> envelope(false, 0,0, 1,0, 2,0, 1,1, 0,2)); } @Test
