You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by js...@apache.org on 2026/10/16 07:05:08 UTC
(sis) 02/08: feat(Geometry): implement envelopes on Arc types
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 <jo...@geomatys.com>
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