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 6de6d0536ab7cd0d0debeaf07de279b490c3c324 Author: jsorel <[email protected]> AuthorDate: Thu Oct 15 12:14:48 2026 +0200 feat(Geometry): implement sphere and ellipsoid horizontal and vertical curves --- .../main/org/apache/sis/geometries/Curve.java | 27 ++-- .../internal/shared/DefaultBSplineCurve.java | 32 +++++ .../geometries/internal/shared/DefaultConic.java | 33 ++++- .../internal/shared/DefaultConicArc.java | 112 ++++++++++----- .../internal/shared/DefaultEllipsoid.java | 19 ++- .../internal/shared/DefaultEllipticArc.java | 18 +-- .../internal/shared/DefaultNurbCurve.java | 16 +++ .../internal/shared/DefaultReversedCurve.java | 24 ++++ .../geometries/internal/shared/DefaultSphere.java | 19 ++- .../internal/shared/EllipsoidalSurfaces.java | 159 +++++++++++++++++++++ .../apache/sis/geometries/surface/Ellipsoid.java | 3 +- .../org/apache/sis/geometries/surface/Sphere.java | 5 +- .../sis/geometries/curve/BSplineCurveTest.java | 25 +++- .../sis/geometries/curve/EllipticArcTest.java | 68 +++++++++ .../sis/geometries/curve/LineStringTest.java | 9 +- .../apache/sis/geometries/curve/NurbCurveTest.java | 13 ++ .../internal/shared/DefaultReversedCurveTest.java | 16 ++- .../sis/geometries/surface/EllipsoidTest.java | 122 ++++++++++++++++ .../apache/sis/geometries/surface/SphereTest.java | 142 ++++++++++++++++++ 19 files changed, 783 insertions(+), 79 deletions(-) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Curve.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Curve.java index b5ac88ed88..d85ee70f5a 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Curve.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Curve.java @@ -272,14 +272,21 @@ public sealed interface Curve extends Orientable * <li>In a composite curve, the last knot of a segment is the first knot of the next segment.</li> * </ul> * + * <p>By default, the knot space is the sequence of integers starting at 0: the knot of + * each data point is its index. Curves with another construction parameterization, + * such as splines, override this method.</p> + * * @return knot values in the construction (knot) space. * * @see ISO 19107:2019 - 6.4.18.4 */ @UML(identifier="knot", specification=ISO_19107) default double[] getKnots() { - //TODO - throw new UnsupportedOperationException(); + final double[] knots = new double[getDataPoints().size()]; + for (int i = 0; i < knots.length; i++) { + knots[i] = i; + } + return knots; } /** @@ -310,14 +317,16 @@ public sealed interface Curve extends Orientable * necessarily smaller than the {@linkplain #getEndConstrParam() end parameter}.</li> * </ul> * - * @return construction parameter of the start point. + * <p>By default, this is the first {@linkplain #getKnots() knot}.</p> + * + * @return construction parameter of the start point, or NaN if this curve has no knot. * * @see ISO 19107:2019 - 6.4.18.15 */ @UML(identifier="startConstrParam", specification=ISO_19107) default double getStartConstrParam() { - //TODO - throw new UnsupportedOperationException(); + final double[] knots = getKnots(); + return (knots.length != 0) ? knots[0] : Double.NaN; } /** @@ -330,14 +339,16 @@ public sealed interface Curve extends Orientable * necessarily greater than the {@linkplain #getStartConstrParam() start parameter}.</li> * </ul> * - * @return construction parameter of the end point. + * <p>By default, this is the last {@linkplain #getKnots() knot}.</p> + * + * @return construction parameter of the end point, or NaN if this curve has no knot. * * @see ISO 19107:2019 - 6.4.18.15 */ @UML(identifier="endConstrParam", specification=ISO_19107) default double getEndConstrParam() { - //TODO - throw new UnsupportedOperationException(); + final double[] knots = getKnots(); + return (knots.length != 0) ? knots[knots.length - 1] : Double.NaN; } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBSplineCurve.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBSplineCurve.java index 099c50ba8d..005174a4b0 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBSplineCurve.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBSplineCurve.java @@ -132,6 +132,38 @@ public non-sealed class DefaultBSplineCurve extends AbstractGeometry implements return knots.clone(); } + /** + * Returns the start of the curve domain: the first knot, or the knot at index {@code degree} of a + * complete knot vector, or 0 if no knot is given since the knots are then uniform between 0 and 1. + */ + @Override + public double getStartConstrParam() { + return domain()[0]; + } + + /** + * Returns the end of the curve domain: the last knot, or the knot at the index of the number of + * control points in a complete knot vector, or 1 if no knot is given. + */ + @Override + public double getEndConstrParam() { + return domain()[1]; + } + + /** + * Returns the start and end of the curve domain, as used by the evaluator. + */ + private double[] domain() { + if (knots.length == 0) { + return new double[] {0, 1}; + } + final int count = (controlPoints != null) ? Math.toIntExact(controlPoints.getLength()) : points.size(); + if (knots.length == count + degree + 1) { + return new double[] {knots[degree], knots[count]}; + } + return new double[] {knots[0], knots[knots.length - 1]}; + } + /** * Returns the number of intervals in the knot array, which is the number of arcs of this curve. */ 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 8b7b5953f2..3e10240c0b 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 @@ -149,7 +149,25 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { } /** - * Returns the arcs of this conic, made of five data points each, consecutive arcs sharing a point. + * Returns the number of data points of each arc: 5 for a general conic. + */ + protected int pointsPerArc() { + return 5; + } + + /** + * Creates the arc going through the given data points. + * + * @param p the {@link #pointsPerArc()} data points of the arc. + * @return the arc going through the given points. + * @throws IllegalStateException if the points do not define a valid arc. + */ + protected DefaultConicArc arc(final Vector<?>[] p) { + return new DefaultConicArc(p, false); + } + + /** + * Returns the arcs of this conic, made of {@link #pointsPerArc()} data points each, consecutive arcs sharing a point. * * @throws IllegalStateException if the data points cannot be split in arcs, or do not define valid arcs. */ @@ -161,17 +179,18 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { for (int i = 0; i < points.size(); i++) { positions.add(points.getPosition(i).copy(DataType.DOUBLE)); } - if (cycle && positions.size() % 4 == 0) { + final int size = pointsPerArc(), step = size - 1; + if (cycle && positions.size() % step == 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 + '.'); + if (n < size || (n - 1) % step != 0) { + throw new IllegalStateException("This curve is made of arcs of " + size + " data points, consecutive arcs" + + " sharing a point, so it needs " + step + "k+1 data points, but got " + n + '.'); } final var arcs = new ArrayList<DefaultConicArc>(); - for (int i = 0; i + 4 < n; i += 4) { - arcs.add(new DefaultConicArc(positions.subList(i, i + 5).toArray(Vector<?>[]::new))); + for (int i = 0; i + step < n; i += step) { + arcs.add(arc(positions.subList(i, i + size).toArray(Vector<?>[]::new))); } return arcs; } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConicArc.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConicArc.java index 0bc3273ece..1fdaed1aa6 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConicArc.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConicArc.java @@ -23,14 +23,18 @@ import org.apache.sis.maths.Vectors; /** - * A conic arc going through five points, computed in the plane of those points. This is one arc of the chain + * A conic arc going through its data points, computed in the plane of those points. This is one arc of the chain * of arcs of a {@link DefaultConic}, the coordinates being handled as Cartesian ones. * + * <p>A general conic arc goes through five points. An elliptic arc goes through four points, since its equation + * has no cross term in the basis of its plane whose first axis is the projection of the first coordinate system + * axis which is not perpendicular to that plane.</p> + * * @author Johann Sorel (Geomatys) */ public final class DefaultConicArc implements ConicArc { /** - * The five data points of the arc. + * The data points of the arc: five for a general conic, four for an elliptic arc. */ private final Vector<?>[] p; @@ -61,7 +65,7 @@ public final class DefaultConicArc implements ConicArc { * Direction of the tangent at p₀, and directions from p₀ to the data points relative to the tangent. */ private double tangent; - private final double[] r = new double[5]; + private final double[] r; /** * Whether the directions from p₀ to the next data points increase. @@ -71,12 +75,19 @@ public final class DefaultConicArc implements ConicArc { /** * Computes the conic through the given points. * - * @param p the five data points of the arc, in their order along the arc. + * @param p the data points of the arc, in their order along the arc. + * @param elliptic {@code true} for an elliptic arc of four points, {@code false} for a conic arc of five points. * @throws IllegalStateException if the points 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. + * conic (or an ellipse), are not in their order along the conic, or are joined through a point at infinity. */ - DefaultConicArc(final Vector<?>[] p) { + DefaultConicArc(final Vector<?>[] p, final boolean elliptic) { + final int[] terms = elliptic ? new int[] {0, 2, 3, 4, 5} : new int[] {0, 1, 2, 3, 4, 5}; + if (p.length != terms.length - 1) { + throw new IllegalArgumentException("Expected " + (terms.length - 1) + " points, but got " + p.length + '.'); + } + final String what = elliptic ? "The four data points of an elliptic arc" : "The five data points of a conic arc"; this.p = p; + this.r = new double[p.length]; final int dimension = p[0].getDimension(); origin = Vectors.createDouble(dimension); for (final Vector<?> q : p) { @@ -111,6 +122,23 @@ public final class DefaultConicArc implements ConicArc { } linear = false; e2.scale(1 / width); + if (elliptic) { + /* + * The axes of the ellipse are along the projection on the plane of the first coordinate system axis + * which is not perpendicular to the plane, and along the perpendicular direction in the plane. + */ + for (int i = 0; i < dimension; i++) { + final double c1 = e1.get(i), c2 = e2.get(i); + final double h = Math.hypot(c1, c2); + if (h > 1E-10) { + final Vector<?> u1 = e1.copy(DataType.DOUBLE).scale( c1 / h).add(e2.copy(DataType.DOUBLE).scale(c2 / h)); + final Vector<?> u2 = e1.copy(DataType.DOUBLE).scale(-c2 / h).add(e2.copy(DataType.DOUBLE).scale(c1 / h)); + e1 = u1; + e2 = u2; + break; + } + } + } 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++) { @@ -118,50 +146,60 @@ public final class DefaultConicArc implements ConicArc { 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."); + throw new IllegalStateException(what + " are not coplanar."); } } /* - * Coefficients of the conic through the five points, - * as the cofactors of the 5×6 matrix whose rows are (a², ab, b², a, b, 1) for each point. + * Coefficients of the conic through the points, as the cofactors of the matrix whose rows are the terms + * (a², ab, b², a, b, 1) for each point. The ab term is omitted for an elliptic arc. */ - 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 int m = terms.length; + final Vector<?> solution = Vectors.createDouble(m); + for (int j = 0; j < m; j++) { + final double[][] minor = new double[m - 1][]; + for (int k = 0; k < m - 1; 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]; + final double[] values = new double[m - 1]; + for (int c = 0, d = 0; c < m; c++) { + if (c != j) values[d++] = row[terms[c]]; } - minor[k] = m; + minor[k] = values; } - coef.set(j, ((j & 1) == 0 ? 1 : -1) * determinant(minor)); + solution.set(j, ((j & 1) == 0 ? 1 : -1) * determinant(minor)); + } + if (!(solution.length() > 1E-12)) { + throw new IllegalStateException(what + " do not determine a single conic." + (elliptic + ? " Two pairs of points having the same coordinate along an axis of the ellipse may be the cause." : "")); } - if (!(coef.length() > 1E-12)) { - throw new IllegalStateException("The five data points of a conic arc do not determine a single conic."); + solution.normalize(); + final double[] coef = new double[6]; + for (int c = 0; c < m; c++) { + coef[terms[c]] = solution.get(c); } - coef.normalize(); - A = coef.get(0); B = coef.get(1); C = coef.get(2); D = coef.get(3); E = coef.get(4); F = coef.get(5); + A = coef[0]; B = coef[1]; C = coef[2]; D = coef[3]; E = coef[4]; F = coef[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," + throw new IllegalStateException(what + " are on a degenerated conic," + " made of straight lines, because three of them are collinear."); } + if (elliptic && !(A * C > 0)) { + throw new IllegalStateException(what + " are not on an ellipse whose axes are aligned with its plane basis."); + } /* * 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. + * interval of directions going from that tangent to the direction of the last point through the other points. */ + final int last = p.length - 1; a0 = a[0]; b0 = b[0]; tangent = direction(-(B*a0 + 2*C*b0 + E), 2*A*a0 + B*b0 + D); - for (int k = 1; k < 5; k++) { + for (int k = 1; k <= last; k++) { r[k] = relative(direction(a[k] - a0, b[k] - b0), tangent); } - increasing = r[2] < r[4]; - for (int k = 1; k < 4; k++) { + increasing = r[1] < r[last]; + for (int k = 1; k < last; 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."); + throw new IllegalStateException(what + " are not in their order along the conic."); } } /* @@ -188,7 +226,8 @@ public final class DefaultConicArc implements ConicArc { * Returns whether the given direction from p₀, relative to the tangent, is the direction of a point of the arc. */ private boolean onArc(final double rx) { - return increasing ? rx <= r[4] : rx >= r[4]; + final double end = r[r.length - 1]; + return increasing ? rx <= end : rx >= end; } /** @@ -246,14 +285,17 @@ public final class DefaultConicArc implements ConicArc { * in that direction. If the points are collinear, the parameter is the index of the data points. */ CurveLinearizer.Piece piece() { + final double[] t = new double[p.length]; if (linear) { - return new CurveLinearizer.Piece((t) -> { - final int i = Math.max(0, Math.min(3, (int) Math.floor(t))); - return p[i].copy(DataType.DOUBLE).lerp(p[i + 1], t - i); - }, new double[] {0, 1, 2, 3, 4}, p); + for (int k = 0; k < t.length; k++) { + t[k] = k; + } + return new CurveLinearizer.Piece((index) -> { + final int i = Math.max(0, Math.min(p.length - 2, (int) Math.floor(index))); + return p[i].copy(DataType.DOUBLE).lerp(p[i + 1], index - i); + }, t, p); } - final double[] t = new double[5]; - for (int k = 1; k < 5; k++) { + for (int k = 1; k < t.length; k++) { t[k] = increasing ? r[k] : Math.PI - r[k]; } return new CurveLinearizer.Piece((angle) -> { diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipsoid.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipsoid.java index 8a4992f9a4..6231bdb79e 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipsoid.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipsoid.java @@ -215,19 +215,32 @@ public non-sealed class DefaultEllipsoid extends AbstractGeometry implements Ell throw new UnsupportedOperationException("Not supported yet."); } + /** + * Returns the line of constant latitude, {@code v} being the latitude in radians. + * At the poles, the returned line string is reduced to the pole. + */ @Override public Curve getHorizontalCurve(double v) { - throw new UnsupportedOperationException("Not supported yet."); + return EllipsoidalSurfaces.parallel(getCoordinateReferenceSystem(), getCenter(), semiAxes, v, false); } + /** + * Returns the line of constant longitude from the south pole to the north pole, + * {@code u} being the longitude in radians. + */ @Override public Curve getVerticalCurve(double u) { - throw new UnsupportedOperationException("Not supported yet."); + return EllipsoidalSurfaces.meridian(getCoordinateReferenceSystem(), getCenter(), semiAxes, u, false); } + /** + * Returns the position of this ellipsoid at longitude {@code u} and latitude {@code v}, in radians. + * + * @throws IllegalStateException if the coordinate system has less than three dimensions. + */ @Override public Vector<?> getSurface(double u, double v) { - throw new UnsupportedOperationException("Not supported yet."); + return EllipsoidalSurfaces.position(getCenter(), semiAxes, u, v); } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipticArc.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipticArc.java index 6b216cacf7..5ea670c3ca 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipticArc.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultEllipticArc.java @@ -16,12 +16,10 @@ */ package org.apache.sis.geometries.internal.shared; -import javax.measure.Quantity; -import org.opengis.geometry.Envelope; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.curve.EllipticArc; -import org.apache.sis.geometries.curve.LineString; import org.apache.sis.maths.Array; +import org.apache.sis.maths.Vector; /** @@ -44,20 +42,18 @@ public non-sealed class DefaultEllipticArc extends DefaultConic implements Ellip } /** - * Not yet supported: the arcs of an ellipse are determined by four data points, not by the five - * data points of the arcs of a general conic. + * Returns 4: an elliptic arc is determined by four data points. */ @Override - public Envelope getEnvelope() { - throw new UnsupportedOperationException("Not supported yet."); + protected int pointsPerArc() { + return 4; } /** - * Not yet supported, for the same reason as {@link #getEnvelope()}. + * Creates the elliptic arc going through the given four data points. */ @Override - public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { - throw new UnsupportedOperationException("Not supported yet."); + protected DefaultConicArc arc(final Vector<?>[] p) { + return new DefaultConicArc(p, true); } - } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbCurve.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbCurve.java index b74cb7ab43..c97d14f19e 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbCurve.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbCurve.java @@ -73,6 +73,22 @@ public non-sealed class DefaultNurbCurve extends AbstractGeometry implements Nur return knots.clone(); } + /** + * Returns the start of the curve domain, which is the knot at index {@code degree}. + */ + @Override + public double getStartConstrParam() { + return domainStart(); + } + + /** + * Returns the end of the curve domain, which is the knot at the index of the number of control points. + */ + @Override + public double getEndConstrParam() { + return domainEnd(); + } + @Override public DataPoints getDataPoints() { return points; diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultReversedCurve.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultReversedCurve.java index bc1d1388c2..f63f980739 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultReversedCurve.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultReversedCurve.java @@ -156,6 +156,30 @@ public non-sealed class DefaultReversedCurve extends AbstractGeometry implements return base.getDataPoints(); } + /** + * Returns the knots of the base curve: reversing a curve keeps its construction parameterization. + */ + @Override + public double[] getKnots() { + return base.getKnots(); + } + + /** + * Returns the end construction parameter of the base curve, which is where this curve starts. + */ + @Override + public double getStartConstrParam() { + return base.getEndConstrParam(); + } + + /** + * Returns the start construction parameter of the base curve, which is where this curve ends. + */ + @Override + public double getEndConstrParam() { + return base.getStartConstrParam(); + } + /** * Returns the polyline of the base curve with its points in reverse order. */ diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultSphere.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultSphere.java index 118f9dc4ec..65fc6416d8 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultSphere.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultSphere.java @@ -147,19 +147,32 @@ public non-sealed class DefaultSphere extends AbstractGeometry implements Sphere throw new UnsupportedOperationException("Not supported yet."); } + /** + * Returns the line of constant latitude, {@code v} being the latitude in radians. + * At the poles, the returned line string is reduced to the pole. + */ @Override public Curve getHorizontalCurve(double v) { - throw new UnsupportedOperationException("Not supported yet."); + return EllipsoidalSurfaces.parallel(getCoordinateReferenceSystem(), getCenter(), new double[] {radius, radius, radius}, v, true); } + /** + * Returns the line of constant longitude from the south pole to the north pole, + * {@code u} being the longitude in radians. + */ @Override public Curve getVerticalCurve(double u) { - throw new UnsupportedOperationException("Not supported yet."); + return EllipsoidalSurfaces.meridian(getCoordinateReferenceSystem(), getCenter(), new double[] {radius, radius, radius}, u, true); } + /** + * Returns the position of this sphere at longitude {@code u} and latitude {@code v}, in radians. + * + * @throws IllegalStateException if the coordinate system has less than three dimensions. + */ @Override public Vector<?> getSurface(double u, double v) { - throw new UnsupportedOperationException("Not supported yet."); + return EllipsoidalSurfaces.position(getCenter(), new double[] {radius, radius, radius}, u, v); } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/EllipsoidalSurfaces.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/EllipsoidalSurfaces.java new file mode 100644 index 0000000000..3f18400a95 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/EllipsoidalSurfaces.java @@ -0,0 +1,159 @@ +/* + * 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.internal.shared; + +import java.util.ArrayList; +import java.util.List; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.apache.sis.geometries.Curve; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.maths.Array; +import org.apache.sis.maths.DataType; +import org.apache.sis.maths.NDArrays; +import org.apache.sis.maths.ReadOnly; +import org.apache.sis.maths.SampleSystem; +import org.apache.sis.maths.Vector; + + +/** + * Parameterization shared by the spheres and the ellipsoids, whose axes are along the first three axes + * of the coordinate system. The <var>u</var> parameter is the longitude and the <var>v</var> parameter + * the latitude, both in radians: + * + * <blockquote><var>S</var>(<var>u</var>, <var>v</var>) = <var>c</var> + + * (<var>a</var>⋅cos <var>u</var>⋅cos <var>v</var>, + * <var>b</var>⋅sin <var>u</var>⋅cos <var>v</var>, + * <var>c</var>⋅sin <var>v</var>)</blockquote> + * + * where <var>a</var>, <var>b</var> and <var>c</var> are the semi-axis lengths. The coordinates on the + * axes after the third one are those of the center. The horizontal curves are the lines of constant + * latitude and the vertical curves the lines of constant longitude, from the south pole to the north pole. + * + * @author Johann Sorel (Geomatys) + */ +final class EllipsoidalSurfaces { + /** + * Do not allow instantiation of this class. + */ + private EllipsoidalSurfaces() { + } + + /** + * Returns the position at the given longitude and latitude. + * + * @param center center of the surface. + * @param semiAxes semi-axis lengths along the first three axes. + * @param u longitude in radians. + * @param v latitude in radians. + * @return position on the surface. + * @throws IllegalStateException if the coordinate system has less than three dimensions. + */ + static Vector<?> position(final ReadOnly.Vector<?> center, final double[] semiAxes, final double u, final double v) { + final int dimension = center.getDimension(); + if (dimension < 3) { + throw new IllegalStateException("A sphere or an ellipsoid is a parametric surface in at least three" + + " dimensions, but the coordinate system has " + dimension + '.'); + } + final Vector<?> position = center.copy(DataType.DOUBLE); + position.set(0, center.get(0) + semiAxes[0] * Math.cos(u) * Math.cos(v)); + position.set(1, center.get(1) + semiAxes[1] * Math.sin(u) * Math.cos(v)); + position.set(2, center.get(2) + semiAxes[2] * Math.sin(v)); + return position; + } + + /** + * Returns the line of constant latitude: a circle if the first two semi-axes are equal, + * an ellipse otherwise. At the poles, the line is reduced to a point, returned as a line string + * of two identical positions. + * + * @param crs coordinate reference system of the surface. + * @param center center of the surface. + * @param semiAxes semi-axis lengths along the first three axes. + * @param v latitude in radians. + * @param circle whether to return a {@link org.apache.sis.geometries.curve.Circle} instead of an + * {@link org.apache.sis.geometries.curve.EllipticArc}. + * @return the line of constant latitude. + */ + static Curve parallel(final CoordinateReferenceSystem crs, final ReadOnly.Vector<?> center, + final double[] semiAxes, final double v, final boolean circle) + { + final GeometryFactory factory = GeometryFactory.DEFAULT; + if (Math.abs(Math.cos(v)) < 1E-12) { + final Vector<?> pole = position(center, semiAxes, 0, v); + return factory.createLineString(factory.createDataPoints(array(crs, List.of(pole, pole)))); + } + final Vector<?> middle = position(center, semiAxes, 0, v); + middle.set(0, center.get(0)); + middle.set(1, center.get(1)); + final var points = new ArrayList<Vector<?>>(); + final var centers = new ArrayList<Vector<?>>(); + if (circle) { + // Four quarters of circle, since an arc of a circle shall be less than a half turn. + for (int i = 0; i <= 4; i++) { + points.add(position(center, semiAxes, (i % 4) * Math.PI / 2, v)); + if (i < 4) centers.add(middle); + } + return factory.createCircle(factory.createDataPoints(array(crs, points)), array(crs, centers), null); + } + /* + * Two halves of ellipse of four points each, the first point closing the curve. The longitudes are shifted + * by π/12 because four points made of two pairs of points having the same coordinate along an axis of the + * ellipse do not determine a single ellipse, which happens with shifts that are multiples of π/6. + */ + for (int i = 0; i < 6; i++) { + points.add(position(center, semiAxes, (2*i + 0.5) * Math.PI / 6, v)); + } + return factory.createEllipticArc(factory.createDataPoints(array(crs, points)), array(crs, List.of(middle, middle)), true); + } + + /** + * Returns the line of constant longitude, from the south pole to the north pole: half a great circle + * of a sphere, or half an ellipse of an ellipsoid. + * + * @param crs coordinate reference system of the surface. + * @param center center of the surface. + * @param semiAxes semi-axis lengths along the first three axes. + * @param u longitude in radians. + * @param circle whether to return an {@link org.apache.sis.geometries.curve.Arc} instead of an + * {@link org.apache.sis.geometries.curve.EllipticArc}. + * @return the line of constant longitude. + */ + static Curve meridian(final CoordinateReferenceSystem crs, final ReadOnly.Vector<?> center, + final double[] semiAxes, final double u, final boolean circle) + { + final GeometryFactory factory = GeometryFactory.DEFAULT; + final Vector<?> c = center.copy(DataType.DOUBLE); + if (circle) { + // Two quarters of great circle, from the south pole to the equator then to the north pole. + final List<Vector<?>> points = List.of(position(center, semiAxes, u, -Math.PI / 2), + position(center, semiAxes, u, 0), position(center, semiAxes, u, Math.PI / 2)); + return factory.createArc(factory.createDataPoints(array(crs, points)), array(crs, List.of(c, c)), null, false); + } + // Not symmetric around the equator, for the reason explained in `parallel(…)`. + final List<Vector<?>> points = List.of(position(center, semiAxes, u, -Math.PI / 2), + position(center, semiAxes, u, -Math.PI / 6), position(center, semiAxes, u, Math.PI / 4), + position(center, semiAxes, u, Math.PI / 2)); + return factory.createEllipticArc(factory.createDataPoints(array(crs, points)), array(crs, List.of(c)), false); + } + + /** + * Returns the given positions as an array in the given reference system. + */ + private static Array array(final CoordinateReferenceSystem crs, final List<? extends ReadOnly.Vector<?>> positions) { + return NDArrays.of(positions, SampleSystem.of(crs), DataType.DOUBLE); + } +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Ellipsoid.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Ellipsoid.java index 1c8639c6f7..bdebeedc99 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Ellipsoid.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Ellipsoid.java @@ -94,7 +94,8 @@ public sealed interface Ellipsoid extends ParametricCurveSurface } /** - * Returns {@link GeometryType#ELLIPTICARC}: the vertical curves are the ellipses of constant longitude. + * Returns {@link GeometryType#ELLIPTICARC}: the vertical curves are the lines of constant longitude, + * which are halves of ellipses going from the south pole to the north pole. */ @Override default GeometryType getVerticalCurveType() { diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Sphere.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Sphere.java index 7fcc4743c8..e070990a28 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Sphere.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/Sphere.java @@ -103,13 +103,14 @@ public sealed interface Sphere extends ParametricCurveSurface } /** - * Returns {@link GeometryType#CIRCLE}: the vertical curves are the great circles of constant longitude. + * Returns {@link GeometryType#ARC}: the vertical curves are the lines of constant longitude, + * which are halves of great circles going from the south pole to the north pole. * * @see ISO 19107:2019 - 8.5.2 */ @Override default GeometryType getVerticalCurveType() { - return GeometryType.CIRCLE; + return GeometryType.ARC; } /** diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BSplineCurveTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BSplineCurveTest.java index b950a83da9..802d743121 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BSplineCurveTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BSplineCurveTest.java @@ -21,6 +21,7 @@ import org.apache.sis.measure.Quantities; import org.apache.sis.measure.Units; // Test dependencies +import static org.junit.jupiter.api.Assertions.assertArrayEquals; import static org.junit.jupiter.api.Assertions.assertEquals; import static org.junit.jupiter.api.Assertions.assertThrows; import org.junit.jupiter.api.Disabled; @@ -45,8 +46,8 @@ public class BSplineCurveTest extends GeometryTest { } @Test - @Disabled("Not implemented yet.") public void testGetKnots() { + assertArrayEquals(new double[] {0, 1, 2}, createPolyline().getKnots()); } @Test @@ -130,14 +131,32 @@ public class BSplineCurveTest extends GeometryTest { public void testGetControlPoints() { } + /** + * The construction parameters span the domain of the curve: the first and last distinct knots, + * [0 … 1] if no knot is given, or the knots at index degree and count of a complete knot vector. + */ @Test - @Disabled("Not implemented yet.") public void testGetStartConstrParam() { + assertEquals(0, createArch().getStartConstrParam()); + assertEquals(0, createPolyline().getStartConstrParam()); + assertEquals(2, createUnclamped().getStartConstrParam()); } @Test - @Disabled("Not implemented yet.") public void testGetEndConstrParam() { + assertEquals(1, createArch().getEndConstrParam()); + assertEquals(2, createPolyline().getEndConstrParam()); + assertEquals(5, createUnclamped().getEndConstrParam()); + } + + /** + * Creates a quadratic b-spline of 5 control points with a complete, uniform and unclamped knot vector. + * Its domain is from the knot at index 2 to the knot at index 5. + */ + private static BSplineCurve createUnclamped() { + return GeometryFactory.DEFAULT.createBSplineCurve( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 1,2, 2,0, 3,2, 4,0), null, + new double[] {0, 1, 2, 3, 4, 5, 6, 7}, 2, null, null, false); } @Test diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/EllipticArcTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/EllipticArcTest.java index afbe307900..b6b1ced3fb 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/EllipticArcTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/EllipticArcTest.java @@ -16,7 +16,12 @@ */ package org.apache.sis.geometries.curve; +import org.opengis.geometry.Envelope; +import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.GeometryFactory; + // Test dependencies +import static org.junit.jupiter.api.Assertions.*; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -32,4 +37,67 @@ public class EllipticArcTest extends GeometryTest { @Disabled("Not implemented yet.") public void testGetInterpolation() { } + + /** + * 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 EllipticArc arc(final boolean cycle, final double... coordinates) { + return GeometryFactory.DEFAULT.createEllipticArc(GeometryFactory.DEFAULT.createDataPoints(CRS_2D, coordinates), null, cycle); + } + + private static void assertEnvelope(final Envelope e, final double xmin, final double ymin, final double xmax, final double ymax) { + assertEquals(xmin, e.getMinimum(0), TOLERANCE, "xmin"); + assertEquals(ymin, e.getMinimum(1), TOLERANCE, "ymin"); + assertEquals(xmax, e.getMaximum(0), TOLERANCE, "xmax"); + assertEquals(ymax, e.getMaximum(1), TOLERANCE, "ymax"); + } + + /** + * An arc of four points crossing the extremes of the ellipse along x at t = 0 and along y at t = π/2. + */ + @Test + public void testGetEnvelope() { + assertEnvelope(arc(false, ellipse(-0.3, 0.6, 1.5, 2.6)).getEnvelope(), 2 * Math.cos(2.6), Math.sin(-0.3), 2, 1); + // Complete ellipse made of two arcs, the first data point closing the second arc. + assertEnvelope(arc(true, ellipse(0.1, 1.1, 2.1, 3.1, 4.1, 5.1)).getEnvelope(), -2, -1, 2, 1); + } + + @Test + public void testAsLine() { + final DataPoints points = arc(false, ellipse(-0.3, 0.6, 1.5, 2.6)).asLine(null, null).getDataPoints(); + assertTrue(points.size() > 10); + double previous = Double.NEGATIVE_INFINITY; + for (int i = 0; i < points.size(); i++) { + final double[] p = points.getPosition(i).toArrayDouble(); + assertEquals(1, p[0]*p[0]/4 + p[1]*p[1], 1E-9, "Point not on the ellipse."); + final double angle = Math.atan2(p[1], p[0] / 2); + assertTrue(angle > previous, "Points not in their order along the ellipse."); + previous = angle; + } + } + + @Test + public void testInvalid() { + // Five points are not a whole number of arcs of four points. + assertThrows(IllegalStateException.class, () -> arc(false, ellipse(0, 0.5, 1, 1.5, 2)).getEnvelope()); + // Four points of the hyperbola x² − y² = 1, which is not an ellipse. + final double[] u = {-1, -0.3, 0.4, 1.1}; + final double[] h = new double[8]; + for (int i = 0; i < u.length; i++) { + h[2*i] = Math.cosh(u[i]); + h[2*i + 1] = Math.sinh(u[i]); + } + assertThrows(IllegalStateException.class, () -> arc(false, h).getEnvelope()); + // Two pairs of points having the same y coordinate: the pair of lines y = 0 and y = sin(π/3) goes through them too. + assertThrows(IllegalStateException.class, () -> arc(false, ellipse(0, Math.PI / 3, 2 * Math.PI / 3, Math.PI)).getEnvelope()); + } } diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/LineStringTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/LineStringTest.java index 533239332c..609ae66704 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/LineStringTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/LineStringTest.java @@ -187,19 +187,22 @@ public abstract class LineStringTest extends GeometryTest { assertEquals(2, ends.getNumGeometries()); } + /** + * The knots of a line string are the indices of its data points. + */ @Test - @Disabled("Not implemented yet.") public void testGetKnots() { + assertArrayEquals(new double[] {0, 1, 2}, createLineString(CRS_2D, 0,0, 10,0, 10,10).getKnots()); } @Test - @Disabled("Not implemented yet.") public void testGetStartConstrParam() { + assertEquals(0, createLineString(CRS_2D, 0,0, 10,0, 10,10).getStartConstrParam()); } @Test - @Disabled("Not implemented yet.") public void testGetEndConstrParam() { + assertEquals(2, createLineString(CRS_2D, 0,0, 10,0, 10,10).getEndConstrParam()); } @Test diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/NurbCurveTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/NurbCurveTest.java index 1cca871fb5..01cdd54c7c 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/NurbCurveTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/NurbCurveTest.java @@ -41,6 +41,19 @@ import org.apache.sis.geometries.GeometryTest; * @author Johann Sorel (Geomatys) */ public class NurbCurveTest extends GeometryTest { + /** + * The construction parameters span the domain of the curve, from the knot at index degree + * to the knot at index count, which differ from the first and last knots of an unclamped curve. + */ + @Test + public void testConstrParamDomain() { + final NurbCurve curve = GeometryFactory.DEFAULT.createNurbCurve( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 1,2, 2,0), + new double[] {1, 1, 1}, new double[] {0, 1, 2, 3, 4, 5}, 2); + assertEquals(2, curve.getStartConstrParam()); + assertEquals(3, curve.getEndConstrParam()); + } + @Test @Disabled("Not implemented yet.") public void testGetGeometryType() { diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/internal/shared/DefaultReversedCurveTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/internal/shared/DefaultReversedCurveTest.java index 6260a35a83..f728e2758c 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/internal/shared/DefaultReversedCurveTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/internal/shared/DefaultReversedCurveTest.java @@ -24,6 +24,7 @@ import org.apache.sis.measure.Units; // Test dependencies import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertArrayEquals; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -88,9 +89,14 @@ public class DefaultReversedCurveTest extends GeometryTest { public void testGetControlPoints() { } + /** + * A reversed curve keeps the construction parameterization of its base curve. + */ @Test - @Disabled("Not implemented yet.") public void testGetKnots() { + final LineString base = GeometryFactory.DEFAULT.createLineString( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 10,0, 10,10)); + assertArrayEquals(base.getKnots(), GeometryFactory.DEFAULT.createReversed(base).getKnots()); } @Test @@ -99,13 +105,17 @@ public class DefaultReversedCurveTest extends GeometryTest { } @Test - @Disabled("Not implemented yet.") public void testGetStartConstrParam() { + final LineString base = GeometryFactory.DEFAULT.createLineString( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 10,0, 10,10)); + assertEquals(2, GeometryFactory.DEFAULT.createReversed(base).getStartConstrParam()); } @Test - @Disabled("Not implemented yet.") public void testGetEndConstrParam() { + final LineString base = GeometryFactory.DEFAULT.createLineString( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 10,0, 10,10)); + assertEquals(0, GeometryFactory.DEFAULT.createReversed(base).getEndConstrParam()); } @Test diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/surface/EllipsoidTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/surface/EllipsoidTest.java new file mode 100644 index 0000000000..bcdc4ab045 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/surface/EllipsoidTest.java @@ -0,0 +1,122 @@ +/* + * 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.surface; + +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.apache.sis.geometries.Curve; +import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.GeometryType; +import org.apache.sis.geometries.curve.EllipticArc; +import org.apache.sis.maths.Vector; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; +import org.apache.sis.geometries.GeometryTest; +import static org.apache.sis.geometries.surface.SphereTest.assertEnvelope; +import static org.apache.sis.geometries.surface.SphereTest.assertPosition; +import static org.apache.sis.geometries.surface.SphereTest.line; + + +/** + * Tests the parameterization of {@link Ellipsoid}: u is the longitude and v the latitude, in radians. + * + * @author Johann Sorel (Geomatys) + */ +public class EllipsoidTest extends GeometryTest { + + private static final CoordinateReferenceSystem CRS = Geometries.getUndefinedCRS(3); + + /** + * Creates the ellipsoid of semi-axes 3, 2 and 1 centered on the origin. + */ + private static Ellipsoid ellipsoid() { + return GeometryFactory.DEFAULT.createEllipsoid(GeometryFactory.DEFAULT.createDataPoints(CRS, 0, 0, 0), 3, 2, 1); + } + + /** + * Returns the value of the ellipsoid equation at the given position, which is 1 on the ellipsoid. + */ + private static double equation(final double[] p) { + return p[0]*p[0]/9 + p[1]*p[1]/4 + p[2]*p[2]; + } + + @Test + public void testGetSurface() { + final Ellipsoid ellipsoid = ellipsoid(); + assertPosition(3, 0, 0, ellipsoid.getSurface(0, 0)); + assertPosition(0, 2, 0, ellipsoid.getSurface(Math.PI / 2, 0)); + assertPosition(0, 0, 1, ellipsoid.getSurface(0.7, Math.PI / 2)); + final Vector<?> p = ellipsoid.getSurface(0.7, 0.3); + assertEquals(1, equation(p.toArrayDouble()), TOLERANCE); + } + + /** + * The horizontal curves are the ellipses of constant latitude. + */ + @Test + public void testGetHorizontalCurve() { + final Curve curve = ellipsoid().getHorizontalCurve(Math.PI / 6); + assertInstanceOf(EllipticArc.class, curve); + assertEquals(GeometryType.ELLIPTICARC, ellipsoid().getHorizontalCurveType()); + final double c = Math.cos(Math.PI / 6); + for (final double[] p : line(curve)) { + assertEquals(0.5, p[2], TOLERANCE, "Point not at the latitude of the curve."); + assertEquals(1, equation(p), 1E-9, "Point not on the ellipsoid."); + } + assertEnvelope(curve.getEnvelope(), -3*c, -2*c, 0.5, 3*c, 2*c, 0.5); + } + + @Test + public void testGetHorizontalCurveAtPole() { + for (final double[] p : line(ellipsoid().getHorizontalCurve(-Math.PI / 2))) { + assertArrayEquals(new double[] {0, 0, -1}, p, TOLERANCE); + } + } + + /** + * The vertical curves are the half ellipses of constant longitude, from the south pole to the north pole. + */ + @Test + public void testGetVerticalCurve() { + final Curve curve = ellipsoid().getVerticalCurve(Math.PI / 4); + assertInstanceOf(EllipticArc.class, curve); + assertPosition(0, 0, -1, curve.getStartPoint().getPosition()); + assertPosition(0, 0, 1, curve.getEndPoint().getPosition()); + final double s = Math.sqrt(0.5); + for (final double[] p : line(curve)) { + assertEquals(1, equation(p), 1E-9, "Point not on the ellipsoid."); + assertEquals(2 * p[0] / 3, p[1], 1E-9, "Point not at the longitude of the curve."); + } + assertEnvelope(curve.getEnvelope(), 0, 0, -1, 3*s, 2*s, 1); + } + + /** + * At longitude 90°, the first axis is perpendicular to the plane of the curve, + * so the axes of the ellipse are the second and the third coordinate system axes. + */ + @Test + public void testGetVerticalCurvePerpendicularToFirstAxis() { + final Curve curve = ellipsoid().getVerticalCurve(Math.PI / 2); + for (final double[] p : line(curve)) { + assertEquals(0, p[0], TOLERANCE); + assertEquals(1, equation(p), 1E-9, "Point not on the ellipsoid."); + } + assertEnvelope(curve.getEnvelope(), 0, 0, -1, 0, 2, 1); + } +} diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/surface/SphereTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/surface/SphereTest.java new file mode 100644 index 0000000000..9c49cd6439 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/surface/SphereTest.java @@ -0,0 +1,142 @@ +/* + * 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.surface; + +import org.opengis.geometry.Envelope; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.apache.sis.geometries.Curve; +import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.GeometryType; +import org.apache.sis.geometries.curve.Arc; +import org.apache.sis.geometries.curve.Circle; +import org.apache.sis.maths.Vector; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; +import org.apache.sis.geometries.GeometryTest; + + +/** + * Tests the parameterization of {@link Sphere}: u is the longitude and v the latitude, in radians. + * + * @author Johann Sorel (Geomatys) + */ +public class SphereTest extends GeometryTest { + + private static final CoordinateReferenceSystem CRS = Geometries.getUndefinedCRS(3); + + /** + * Creates the sphere of radius 2 centered on (1, 2, 3). + */ + private static Sphere sphere() { + return GeometryFactory.DEFAULT.createSphere(GeometryFactory.DEFAULT.createDataPoints(CRS, 1, 2, 3), 2); + } + + static void assertPosition(final double x, final double y, final double z, final Vector<?> p) { + assertEquals(x, p.get(0), TOLERANCE, "x"); + assertEquals(y, p.get(1), TOLERANCE, "y"); + assertEquals(z, p.get(2), TOLERANCE, "z"); + } + + static void assertEnvelope(final Envelope e, final double... bounds) { + for (int i = 0; i < 3; i++) { + assertEquals(bounds[i], e.getMinimum(i), TOLERANCE, "minimum"); + assertEquals(bounds[i + 3], e.getMaximum(i), TOLERANCE, "maximum"); + } + } + + /** + * Returns the positions of the linear approximation of the given curve. + */ + static double[][] line(final Curve curve) { + final DataPoints points = curve.asLine(null, null).getDataPoints(); + final double[][] positions = new double[points.size()][]; + for (int i = 0; i < positions.length; i++) { + positions[i] = points.getPosition(i).toArrayDouble(); + } + return positions; + } + + @Test + public void testGetSurface() { + final Sphere sphere = sphere(); + assertPosition(3, 2, 3, sphere.getSurface(0, 0)); + assertPosition(1, 4, 3, sphere.getSurface(Math.PI / 2, 0)); + assertPosition(-1, 2, 3, sphere.getSurface(Math.PI, 0)); + assertPosition(1, 2, 5, sphere.getSurface(0.7, Math.PI / 2)); + assertPosition(1, 2, 1, sphere.getSurface(0.7, -Math.PI / 2)); + final Vector<?> p = sphere.getSurface(0.7, 0.3); + assertEquals(2, Math.sqrt(Math.pow(p.get(0) - 1, 2) + Math.pow(p.get(1) - 2, 2) + Math.pow(p.get(2) - 3, 2)), TOLERANCE); + } + + /** + * A sphere in two dimensions is a circle, which has no position for two parameters. + */ + @Test + public void testGetSurface2D() { + final Sphere circle = GeometryFactory.DEFAULT.createSphere( + GeometryFactory.DEFAULT.createDataPoints(Geometries.getUndefinedCRS(2), 0, 0), 1); + assertThrows(IllegalStateException.class, () -> circle.getSurface(0, 0)); + } + + /** + * The horizontal curves are the circles of constant latitude. + */ + @Test + public void testGetHorizontalCurve() { + final Curve curve = sphere().getHorizontalCurve(Math.PI / 6); + assertInstanceOf(Circle.class, curve); + assertEquals(GeometryType.CIRCLE, sphere().getHorizontalCurveType()); + final double r = Math.sqrt(3); // 2⋅cos(30°) + for (final double[] p : line(curve)) { + assertEquals(4, p[2], TOLERANCE, "Point not at the latitude of the curve."); + assertEquals(r, Math.hypot(p[0] - 1, p[1] - 2), TOLERANCE, "Point not on the circle."); + } + assertEnvelope(curve.getEnvelope(), 1 - r, 2 - r, 4, 1 + r, 2 + r, 4); + } + + /** + * At a pole, the circle of constant latitude is reduced to the pole. + */ + @Test + public void testGetHorizontalCurveAtPole() { + for (final double[] p : line(sphere().getHorizontalCurve(Math.PI / 2))) { + assertArrayEquals(new double[] {1, 2, 5}, p, TOLERANCE); + } + } + + /** + * The vertical curves are the half great circles of constant longitude, from the south pole to the north pole. + */ + @Test + public void testGetVerticalCurve() { + final Curve curve = sphere().getVerticalCurve(Math.PI / 2); + assertInstanceOf(Arc.class, curve); + assertEquals(GeometryType.ARC, sphere().getVerticalCurveType()); + assertPosition(1, 2, 1, curve.getStartPoint().getPosition()); + assertPosition(1, 2, 5, curve.getEndPoint().getPosition()); + for (final double[] p : line(curve)) { + assertEquals(1, p[0], TOLERANCE, "Point not at the longitude of the curve."); + assertTrue(p[1] >= 2 - TOLERANCE, "Point on the other half of the great circle."); + assertEquals(2, Math.hypot(p[1] - 2, p[2] - 3), TOLERANCE, "Point not on the sphere."); + } + assertEnvelope(curve.getEnvelope(), 1, 2, 1, 1, 4, 5); + } +}
