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 4f0239d852ac7f75dd8f3bc0c7b943962a915802 Author: jsorel <[email protected]> AuthorDate: Thu Oct 15 11:15:32 2026 +0200 feat(Geometry): implement Curve.asLine on multiple types --- .../main/org/apache/sis/geometries/Curve.java | 3 +- .../main/org/apache/sis/geometries/curve/Line.java | 7 +- .../apache/sis/geometries/curve/LineString.java | 12 +- .../internal/shared/BSplineEvaluator.java | 64 +++ .../internal/shared/CurveLinearizer.java | 493 +++++++++++++++++++++ .../sis/geometries/internal/shared/DefaultArc.java | 89 ++++ .../internal/shared/DefaultArcByBulge.java | 34 +- .../internal/shared/DefaultArcByCenterPoint.java | 28 ++ .../internal/shared/DefaultBSplineCurve.java | 13 + .../geometries/internal/shared/DefaultBezier.java | 12 + .../internal/shared/DefaultCircularString.java | 22 + .../internal/shared/DefaultCompoundCurve.java | 11 + .../geometries/internal/shared/DefaultConic.java | 259 ++--------- .../{DefaultConic.java => DefaultConicArc.java} | 228 +++++----- .../internal/shared/DefaultEllipticArc.java | 20 + .../internal/shared/DefaultNurbCurve.java | 13 + .../internal/shared/DefaultReversedCurve.java | 9 + .../internal/shared/CurveLinearizerTest.java | 266 +++++++++++ 18 files changed, 1229 insertions(+), 354 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 4303695aa9..b5ac88ed88 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 @@ -506,8 +506,7 @@ public sealed interface Curve extends Orientable */ @UML(identifier="asLine", specification=ISO_19107) default LineString asLine(Quantity<?> spacing, Quantity<?> offset) { - //TODO - throw new UnsupportedOperationException(); + throw new UnsupportedOperationException("Linear approximation of " + getGeometryType() + " curves is not supported yet."); } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/Line.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/Line.java index 1298ae7230..9bee4f89a2 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/Line.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/Line.java @@ -29,6 +29,7 @@ import org.apache.sis.geometries.GeometryType; import org.apache.sis.geometries.Point; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.internal.shared.CurveLinearizer; import org.apache.sis.geometries.internal.shared.Curves; import org.apache.sis.geometries.internal.shared.DefaultLine; import org.apache.sis.geometries.internal.shared.DefaultLineString; @@ -88,13 +89,15 @@ public sealed interface Line extends Curve } /** - * Returns {@code this}, since a line string is already a linear approximation of itself. + * Returns a line string with the points of this line, plus points added on its segments if needed + * for meeting the spacing. The attributes of the added points are interpolated between the ends of + * their segment. The offset is ignored, since the returned line has the shape of this line. * * @see ISO 19107:2019 - 6.4.18.17 */ @Override default LineString asLine(Quantity<?> spacing, Quantity<?> offset) { - return new DefaultLineString(getDataPoints()); + return new DefaultLineString(CurveLinearizer.densify(this, getDataPoints(), spacing)); } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/LineString.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/LineString.java index cffc5b9109..cc3f52d9c9 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/LineString.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/LineString.java @@ -28,6 +28,8 @@ import org.apache.sis.geometries.GeometryType; import org.apache.sis.geometries.Point; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.internal.shared.CurveLinearizer; import org.apache.sis.geometries.internal.shared.Curves; import org.apache.sis.geometries.internal.shared.DefaultLineString; import org.apache.sis.geometries.mesh.MeshPrimitive; @@ -92,13 +94,19 @@ public sealed interface LineString extends Curve } /** - * Returns {@code this}, since a line string is already a linear approximation of itself. + * Returns this line string, with points added on its segments if needed for meeting the spacing. + * The attributes of the added points are interpolated between the ends of their segment. + * The offset is ignored, since the returned line is this line string. + * + * @return {@code this} if no point needs to be added, otherwise a new line string. * * @see ISO 19107:2019 - 6.4.18.17 */ @Override default LineString asLine(Quantity<?> spacing, Quantity<?> offset) { - return this; + final DataPoints points = getDataPoints(); + final DataPoints densified = CurveLinearizer.densify(this, points, spacing); + return (densified == points) ? this : GeometryFactory.DEFAULT.createLineString(densified); } /** diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/BSplineEvaluator.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/BSplineEvaluator.java index 5750e57e7b..378313a124 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/BSplineEvaluator.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/BSplineEvaluator.java @@ -23,16 +23,19 @@ import javax.measure.Quantity; import javax.measure.Unit; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.AttributeInterpolation; +import org.apache.sis.geometries.Curve; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.Point; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometries.operation.GaussLegendreQuadrature; import org.apache.sis.maths.Array; import org.apache.sis.maths.DataType; import org.apache.sis.maths.NDArrays; import org.apache.sis.maths.SampleSystem; import org.apache.sis.maths.Vector; +import org.apache.sis.maths.Vectors; /** @@ -254,6 +257,67 @@ final class BSplineEvaluator { return pointAt(u, k); } + /** + * Returns the distinct knot values from the start to the end of the curve domain. + * The curve is a polynomial or rational function between two consecutive values. + */ + double[] breakpoints() { + final double[] values = new double[last - degree + 2]; + int n = 0; + for (int k = degree; k <= last + 1; k++) { + if (n == 0 || knots[k] > values[n - 1]) { + values[n++] = knots[k]; + } + } + return Arrays.copyOf(values, n); + } + + /** + * Returns a polyline whose points are on the curve, including the positions at the knots. + * + * @param curve the curve evaluated by this object. + * @param spacing maximal distance along the curve between two points, or zero or {@code null} for no limit. + * @param offset maximal distance between the curve and the polyline, or zero or {@code null} for no limit. + * @return the polyline approximating the curve. + */ + LineString asLine(final Curve curve, final Quantity<?> spacing, final Quantity<?> offset) { + return CurveLinearizer.linearize(curve, spacing, offset, + List.of(new CurveLinearizer.Piece(this::position, breakpoints(), null))); + } + + /** + * Returns the position of the curve at parameter {@code u}, without the other attributes. + * + * @param u the parameter, in the curve domain. + * @return the position at the given parameter. + */ + Vector<?> position(final double u) { + int k = degree; + for (int i = degree; i <= last; i++) { + if (knots[i] < knots[i+1] && (u >= knots[i] || i == degree)) { + k = i; + } + } + final double[] R = new double[degree + 1]; + basis(u, k, R, null); + final int dim = positions[0].length; + final double[] sum = new double[dim]; + double W = 0; + for (int r = 0; r <= degree; r++) { + final int i = k - degree + r; + final double w = R[r] * weights[i]; + W += w; + for (int c = 0; c < dim; c++) { + sum[c] += w * positions[i][c]; + } + } + final Vector<?> position = Vectors.create(SampleSystem.of(controls.getCoordinateReferenceSystem()), DataType.DOUBLE); + for (int c = 0; c < dim; c++) { + position.set(c, sum[c] / W); + } + return position; + } + /** * Integrates the speed over [a … b] in span {@code k}, which gives the arc length. */ diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/CurveLinearizer.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/CurveLinearizer.java new file mode 100644 index 0000000000..ec0a2dbff4 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/CurveLinearizer.java @@ -0,0 +1,493 @@ +/* + * 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.Collections; +import java.util.List; +import java.util.function.DoubleFunction; +import javax.measure.Quantity; +import javax.measure.Unit; +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.DataPointsType; +import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.Point; +import org.apache.sis.geometries.curve.LineString; +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; +import org.apache.sis.measure.Quantities; + + +/** + * Approximates curves by polylines whose points lie on the curves, for implementing + * {@link Curve#asLine(Quantity, Quantity)}. + * + * <p>A curve is described by one or more {@linkplain Piece pieces}, each one being a function from a parameter + * to a position, with the parameters of the positions which must be in the polyline (the start and end of the + * piece, the data points, the knots). Each interval between two such parameters is halved until both criteria + * of {@code asLine} are met:</p> + * <ul> + * <li>the arc length of each interval is at most the spacing. It is estimated from the lengths of the chords + * joining 2 and 4 sub-intervals, extrapolated to an infinite number of chords;</li> + * <li>the distance between the curve and the chord of each interval is at most the offset. + * It is estimated from the curve positions at one quarter, half and three quarters of the interval.</li> + * </ul> + * + * <p>If both the spacing and the offset are zero, the offset is a fraction of the extent of the curve.</p> + * + * @author Johann Sorel (Geomatys) + */ +public final class CurveLinearizer { + /** + * Offset used when no criterion is given, relative to the extent of the curve. + */ + private static final double DEFAULT_OFFSET = 1E-3; + + /** + * Offset used for estimating the length of the curve, relative to its extent. + * The length is needed when the spacing or the offset is given in percent of the length. + */ + private static final double LENGTH_OFFSET = 1E-6; + + /** + * Maximal number of interval halvings, which stops the refinement at the precision of the parameter. + */ + private static final int MAX_DEPTH = 50; + + /** + * A part of a curve, as a function from a parameter to a position. + * + * @param function the position on the curve for a given parameter. + * @param parameters the parameters of the positions which shall be in the polyline, in their order along + * the curve. The first and last ones are the start and end of the piece. The values may + * be increasing or decreasing. + * @param positions the exact positions at the given parameters, or {@code null} for computing all of them + * with the function. A {@code null} element is also computed with the function. + */ + public record Piece(DoubleFunction<? extends ReadOnly.Vector<?>> function, double[] parameters, + ReadOnly.Vector<?>[] positions) + { + /** + * Returns the position at the parameter of the given index. + */ + ReadOnly.Vector<?> position(final int i) { + if (positions != null && positions[i] != null) { + return positions[i]; + } + return function.apply(parameters[i]); + } + } + + /** + * Maximal arc length and offset of the intervals, in the unit of the coordinates. Zero means no limit. + */ + private final double spacing, offset; + + /** + * The points of the polyline, in order. + */ + private final List<ReadOnly.Vector<?>> points = new ArrayList<>(); + + private CurveLinearizer(final double spacing, final double offset) { + this.spacing = spacing; + this.offset = offset; + } + + /** + * Returns a polyline approximating a curve made of the given pieces, consecutive pieces sharing their end point. + * + * @param curve the curve to approximate, for its reference system and the unit of its coordinates. + * @param spacing maximal distance along the curve between two points of the polyline, or zero or {@code null} + * for no limit. May be in percent of the curve length. + * @param offset maximal distance between the curve and the polyline, or zero or {@code null} for no limit. + * May be in percent of the curve length. + * @param pieces the pieces of the curve. + * @return the polyline approximating the curve. + * @throws IllegalArgumentException if the spacing or the offset is negative, or cannot be converted to + * the unit of the coordinates. + */ + public static LineString linearize(final Curve curve, final Quantity<?> spacing, final Quantity<?> offset, + final List<Piece> pieces) + { + final Unit<?> unit = Geometries.getLinearUnit(curve); + double length = Double.NaN; + if (isPercent(spacing) || isPercent(offset)) { + length = length(run(0, LENGTH_OFFSET * extent(pieces), pieces)); + } + double s = distance(spacing, unit, length); + double o = distance(offset, unit, length); + if (s == 0 && o == 0) { + o = DEFAULT_OFFSET * extent(pieces); + } + return lineString(curve.getCoordinateReferenceSystem(), run(s, o, pieces)); + } + + /** + * Computes the points of the polyline approximating the given pieces. + */ + private static List<ReadOnly.Vector<?>> run(final double spacing, final double offset, final List<Piece> pieces) { + final var linearizer = new CurveLinearizer(spacing, offset); + for (final Piece piece : pieces) { + final double[] t = piece.parameters(); + ReadOnly.Vector<?> previous = piece.position(0); + if (linearizer.points.isEmpty()) { + linearizer.points.add(previous); + } + for (int i = 1; i < t.length; i++) { + final ReadOnly.Vector<?> next = piece.position(i); + linearizer.refine(piece.function(), t[i-1], previous, t[i], next, 0); + previous = next; + } + } + return linearizer.points; + } + + /** + * Adds the points of the polyline between the given positions, excluding the first one. + */ + private void refine(final DoubleFunction<? extends ReadOnly.Vector<?>> function, + final double ta, final ReadOnly.Vector<?> pa, + final double tb, final ReadOnly.Vector<?> pb, final int depth) + { + final double tm = (ta + tb) / 2; + if (depth < MAX_DEPTH && tm != ta && tm != tb) { + final ReadOnly.Vector<?> pm = function.apply(tm); + final ReadOnly.Vector<?> q1 = function.apply((ta + tm) / 2); + final ReadOnly.Vector<?> q3 = function.apply((tm + tb) / 2); + boolean split = false; + if (spacing > 0) { + final double half = distance(pa, pm) + distance(pm, pb); + final double quarter = distance(pa, q1) + distance(q1, pm) + distance(pm, q3) + distance(q3, pb); + split = quarter + (quarter - half) / 3 > spacing; + } + if (!split && offset > 0) { + split = Math.max(distance(pm, pa, pb), Math.max(distance(q1, pa, pb), distance(q3, pa, pb))) > offset; + } + if (split) { + refine(function, ta, pa, tm, pm, depth + 1); + refine(function, tm, pm, tb, pb, depth + 1); + return; + } + } + points.add(pb); + } + + /** + * Returns the circular arc going from {@code start} to {@code end} through {@code middle}, as a function of the + * angle in radians from the start point. The arc is in the plane of its three points, so it works in any number + * of dimensions. If the three points are collinear, the piece is the straight segment from start to end. + * If the start and end points are identical, the piece is a full circle whose diameter goes from them to + * the middle point, turning counterclockwise. + * + * @param start the start point. + * @param middle a point on the arc, distinct from its start and end points. + * @param end the end point. + * @param keep whether the middle point is a data point which shall be in the polyline. + * @return the arc as a piece of curve. + * @throws IllegalStateException if the arc is a full circle in more than two dimensions, + * since its plane is then not determined. + */ + public static Piece circularArc(final ReadOnly.Vector<?> start, final ReadOnly.Vector<?> middle, + final ReadOnly.Vector<?> end, final boolean keep) + { + final int dimension = start.getDimension(); + final Vector<?> u = end.copy(DataType.DOUBLE).subtract(start); + final double chord = u.length(); + final Vector<?> centre, e1, e2; + final double radius; + 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."); + } + centre = start.copy(DataType.DOUBLE).lerp(middle, 0.5); + radius = distance(start, centre); + e1 = start.copy(DataType.DOUBLE).subtract(centre).scale(1 / radius); + e2 = e1.copy(DataType.DOUBLE); + e2.set(0, -e1.get(1)); + e2.set(1, e1.get(0)); + } else { + u.scale(1 / chord); + final Vector<?> mid = start.copy(DataType.DOUBLE).lerp(end, 0.5); + final Vector<?> n = middle.copy(DataType.DOUBLE).subtract(mid); + final double pu = n.dot(u); + n.perpendicular(u); + final double pn = n.length(); + if (!(pn > chord * 1E-12)) { + // Collinear points: a straight segment. + final double tm = middle.copy(DataType.DOUBLE).subtract(start).dot(u) / chord; + return new Piece((t) -> start.copy(DataType.DOUBLE).lerp(end, t), + keep ? new double[] {0, tm, 1} : new double[] {0, 1}, + keep ? new ReadOnly.Vector<?>[] {start, middle, end} : new ReadOnly.Vector<?>[] {start, end}); + } + n.scale(1 / pn); + // Same circle as in `Curves.addArc(…)`. + final double c = chord / 2; + final double t = (pu*pu + pn*pn - c*c) / (2*pn); + radius = Math.hypot(c, t); + centre = mid.add(n.copy(DataType.DOUBLE).scale(t)); + e1 = start.copy(DataType.DOUBLE).subtract(centre).scale(1 / radius); + e2 = n.perpendicular(e1); + e2.normalize(); + } + /* + * Angles of the middle and end points in the basis (e1, e2), in the [0 … 2π) range. If the middle point + * comes first when turning in the direction of positive angles, the arc turns in that direction. + * Otherwise it turns in the direction of negative angles. + */ + double am = angle(middle, centre, e1, e2); + double ae = (chord == 0) ? 2*Math.PI : angle(end, centre, e1, e2); + if (am > ae) { + am -= 2*Math.PI; + ae -= 2*Math.PI; + } + final DoubleFunction<Vector<?>> function = (a) -> centre.copy(DataType.DOUBLE) + .add(e1.copy(DataType.DOUBLE).scale(radius * Math.cos(a))) + .add(e2.copy(DataType.DOUBLE).scale(radius * Math.sin(a))); + return new Piece(function, + keep ? new double[] {0, am, ae} : new double[] {0, ae}, + keep ? new ReadOnly.Vector<?>[] {start, middle, end} : new ReadOnly.Vector<?>[] {start, end}); + } + + /** + * Returns the angle of the given position around the given centre, in the basis (e1, e2), in the [0 … 2π) range. + */ + private static double angle(final ReadOnly.Vector<?> p, final Vector<?> centre, final Vector<?> e1, final Vector<?> e2) { + final Vector<?> v = p.copy(DataType.DOUBLE).subtract(centre); + final double a = Math.atan2(v.dot(e2), v.dot(e1)); + return (a < 0) ? a + 2*Math.PI : a; + } + + /** + * Returns the distance between two positions. + */ + private static double distance(final ReadOnly.Vector<?> a, final ReadOnly.Vector<?> b) { + return a.copy(DataType.DOUBLE).subtract(b).length(); + } + + /** + * Returns the distance between a position and the segment from {@code a} to {@code b}. + */ + private static double distance(final ReadOnly.Vector<?> p, final ReadOnly.Vector<?> a, final ReadOnly.Vector<?> b) { + final Vector<?> ab = b.copy(DataType.DOUBLE).subtract(a); + final Vector<?> ap = p.copy(DataType.DOUBLE).subtract(a); + final double square = ab.lengthSquare(); + final double t = (square > 0) ? Math.max(0, Math.min(1, ap.dot(ab) / square)) : 0; + return ap.subtract(ab.scale(t)).length(); + } + + /** + * Returns the length of the given polyline. + */ + private static double length(final List<ReadOnly.Vector<?>> points) { + double length = 0; + for (int i = 1; i < points.size(); i++) { + length += distance(points.get(i - 1), points.get(i)); + } + return length; + } + + /** + * Returns the diagonal of the box containing the given pieces, estimated from a few positions of each interval. + */ + private static double extent(final List<Piece> pieces) { + GeometryEnvelope box = null; + for (final Piece piece : pieces) { + final double[] t = piece.parameters(); + for (int i = 0; i < t.length; i++) { + final int samples = (i == 0) ? 1 : 8; + for (int k = samples; --k >= 0;) { + final ReadOnly.Vector<?> p = (k == 0) ? piece.position(i) + : piece.function().apply(t[i] + (t[i-1] - t[i]) * k / samples); + if (box == null) { + box = new GeometryEnvelope(p, p); + } else { + box.add(p); + } + } + } + } + return diagonal(box); + } + + /** + * Returns whether the given distance is in percent of the curve length. + */ + private static boolean isPercent(final Quantity<?> q) { + return q != null && Curves.isPercent(q); + } + + /** + * Converts the given spacing or offset to the unit of the coordinates. + * + * @param q the spacing or the offset, or {@code null}. + * @param unit the unit of the coordinates. + * @param length the length of the curve, used if {@code q} is in percent. + * @return the distance in the unit of the coordinates, or 0 for no limit. + */ + private static double distance(final Quantity<?> q, final Unit<?> unit, final double length) { + if (q == null) { + return 0; + } + final double value = isPercent(q) ? q.getValue().doubleValue() / 100 * length : Curves.convert(q, unit); + if (!(value >= 0) || Double.isInfinite(value)) { + throw new IllegalArgumentException("The spacing and the offset of a linear approximation must be" + + " positive or zero, but got " + q + '.'); + } + return value; + } + + /** + * Returns a polyline approximating a curve made of the given parts, each part ending where the next one starts. + * The polylines of the parts are joined, without repeating the shared points. + * + * @param curve the curve made of the given parts. + * @param parts the parts of the curve, in order. + * @param spacing maximal distance along the curve between two points of the polyline, or zero or {@code null} + * for no limit. May be in percent of the length of the whole curve. + * @param offset maximal distance between the curve and the polyline, or zero or {@code null} for no limit. + * May be in percent of the length of the whole curve. + * @return the polyline approximating the curve. + */ + public static LineString concatenate(final Curve curve, final List<? extends Curve> parts, + Quantity<?> spacing, Quantity<?> offset) + { + if (isPercent(spacing) || isPercent(offset)) { + final Unit<?> unit = Geometries.getLinearUnit(curve); + double length = 0; + for (final Curve part : parts) { + final Quantity<?> precision = Quantities.create(LENGTH_OFFSET * diagonal(part.getEnvelope()), unit); + length += length(positions(part.asLine(null, precision))); + } + spacing = Quantities.create(distance(spacing, unit, length), unit); + offset = Quantities.create(distance(offset, unit, length), unit); + } + final var points = new ArrayList<ReadOnly.Vector<?>>(); + for (final Curve part : parts) { + final List<ReadOnly.Vector<?>> positions = positions(part.asLine(spacing, offset)); + points.addAll(points.isEmpty() ? positions : positions.subList(Math.min(1, positions.size()), positions.size())); + } + return lineString(curve.getCoordinateReferenceSystem(), points); + } + + /** + * Returns the polyline approximating the given curve, with its points in reverse order. + * + * @param curve the reversed curve. + * @param base the curve which is reversed. + * @param spacing maximal distance along the curve between two points of the polyline, or zero or {@code null}. + * @param offset maximal distance between the curve and the polyline, or zero or {@code null}. + * @return the polyline approximating the reversed curve. + */ + public static LineString reverse(final Curve curve, final Curve base, final Quantity<?> spacing, final Quantity<?> offset) { + final List<ReadOnly.Vector<?>> points = positions(base.asLine(spacing, offset)); + Collections.reverse(points); + return lineString(curve.getCoordinateReferenceSystem(), points); + } + + /** + * Returns the positions of the given line string. + */ + private static List<ReadOnly.Vector<?>> positions(final LineString line) { + final DataPoints points = line.getDataPoints(); + final var positions = new ArrayList<ReadOnly.Vector<?>>(points.size()); + for (int i = 0; i < points.size(); i++) { + positions.add(points.getPosition(i).copy(DataType.DOUBLE)); + } + return positions; + } + + /** + * Returns the length of the diagonal of the given envelope, or 0 if there is no envelope. + */ + private static double diagonal(final Envelope envelope) { + double sum = 0; + if (envelope != null) { + for (int i = 0; i < envelope.getDimension(); i++) { + final double span = envelope.getSpan(i); + if (Double.isFinite(span)) { + sum += span * span; + } + } + } + return Math.sqrt(sum); + } + + /** + * Returns a line string made of the given positions. + * + * @param crs reference system of the positions. + * @param points positions of the line string. + * @return the line string. + */ + public static LineString lineString(final CoordinateReferenceSystem crs, final List<? extends ReadOnly.Vector<?>> points) { + final GeometryFactory factory = GeometryFactory.DEFAULT; + return factory.createLineString(factory.createDataPoints(NDArrays.of(points, SampleSystem.of(crs), DataType.DOUBLE))); + } + + /** + * Returns the given polyline with points added on its segments, so that consecutive points are at most at + * the given spacing. The attributes of the added points are interpolated between the ends of their segment. + * + * @param curve the polyline, for its unit of coordinates. + * @param points the points of the polyline. + * @param spacing maximal distance between two consecutive points, or zero or {@code null} for no limit. + * May be in percent of the polyline length. + * @return the points of the densified polyline, or {@code points} if no point needs to be added. + */ + public static DataPoints densify(final Curve curve, final DataPoints points, final Quantity<?> spacing) { + final int n = points.size(); + final double[] segments = new double[Math.max(n - 1, 0)]; + double length = 0; + for (int i = 0; i < segments.length; i++) { + segments[i] = distance(points.getPosition(i), points.getPosition(i + 1)); + length += segments[i]; + } + final double s = distance(spacing, Geometries.getLinearUnit(curve), length); + if (s == 0) { + return points; + } + final var result = new ArrayList<Point>(); + for (int i = 0; i < segments.length; i++) { + final int parts = Math.max(1, (int) Math.ceil(segments[i] / s)); + for (int k = 0; k < parts; k++) { + result.add((k == 0) ? points.getPoint(i) : points.getInterpolatedPoint(i, i + 1, k / (double) parts)); + } + } + if (n != 0) { + result.add(points.getPoint(n - 1)); + } + if (result.size() == n) { + return points; + } + final DataPointsType type = points.getType(); + final DataPoints copy = GeometryFactory.DEFAULT.createDataPoints(type, result.size()); + for (int i = 0; i < result.size(); i++) { + for (final String name : type.getAttributeNames()) { + copy.setAttribute(i, name, result.get(i).getAttribute(name)); + } + } + return copy; + } +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArc.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArc.java index 59aa3b46fe..248465d222 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArc.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultArc.java @@ -16,10 +16,16 @@ */ package org.apache.sis.geometries.internal.shared; +import java.util.ArrayList; import java.util.List; +import javax.measure.Quantity; +import org.opengis.geometry.Envelope; import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.curve.Arc; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.maths.Array; +import org.apache.sis.maths.DataType; import org.apache.sis.maths.Vector; @@ -67,4 +73,87 @@ public non-sealed class DefaultArc extends DefaultConic implements Arc { return radius; } + /** + * Returns the smallest envelope containing this chain of arcs. That envelope contains the data points, + * plus each point where an arc reaches an extreme along an axis. + * + * @throws IllegalStateException if there is not one control point per arc, or if the arc is a half turn + * with no radius vector in more than two dimensions. + */ + @Override + public Envelope getEnvelope() { + if (points.isEmpty()) { + final var envelope = new GeometryEnvelope(getCoordinateReferenceSystem()); + envelope.setToNaN(); + return envelope; + } + final var envelope = new GeometryEnvelope(points.getEnvelope()); + for (final Vector<?>[] arc : circularArcs()) { + Curves.addArc(envelope, arc[0], arc[1], arc[2]); + } + return envelope; + } + + /** + * Returns a polyline whose points are on this chain of arcs, including all its data points. + * + * @throws IllegalStateException for the same reasons as {@link #getEnvelope()}. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + final var pieces = new ArrayList<CurveLinearizer.Piece>(); + for (final Vector<?>[] arc : circularArcs()) { + pieces.add(CurveLinearizer.circularArc(arc[0], arc[1], arc[2], false)); + } + if (pieces.isEmpty()) { + return GeometryFactory.DEFAULT.createLineString(points); + } + return CurveLinearizer.linearize(this, spacing, offset, pieces); + } + + /** + * Returns the start point, a middle point and the end point of each arc. The middle point is at the + * radius vector from the centre if there is one. Otherwise it is the middle of the shorter arc, + * since an arc is less than a half turn. + * + * @throws IllegalStateException if there is not one control point per arc, or if the arc is a half turn + * with no radius vector in more than two dimensions. + */ + private List<Vector<?>[]> circularArcs() { + final int count = getNumArc(); + if (count == 0) { + return List.of(); + } + if (controlPoints == null || controlPoints.getLength() < count) { + throw new IllegalStateException("A chain of " + count + " arcs needs one centre per arc as control points."); + } + final var arcs = new ArrayList<Vector<?>[]>(count); + for (int i = 0; i < count; i++) { + final Vector<?> start = points.getPosition(i).copy(DataType.DOUBLE); + final Vector<?> end = points.getPosition(i + 1).copy(DataType.DOUBLE); + final Vector<?> centre = controlPoints.get(i).copy(DataType.DOUBLE); + final Vector<?> middle; + if (i < radius.size()) { + middle = centre.copy(DataType.DOUBLE).add(radius.get(i)); + } else { + final double r = start.copy(DataType.DOUBLE).subtract(centre).length(); + final Vector<?> bisector = start.copy(DataType.DOUBLE).subtract(centre).add(end).subtract(centre); + final double length = bisector.length(); + if (length > r * 1E-12) { + middle = centre.copy(DataType.DOUBLE).add(bisector.scale(r / length)); + } else if (centre.getDimension() == 2) { + // Half turn: counterclockwise by default. + middle = centre.copy(DataType.DOUBLE); + middle.set(0, centre.get(0) - (start.get(1) - centre.get(1))); + middle.set(1, centre.get(1) + (start.get(0) - centre.get(0))); + } else { + throw new IllegalStateException("Arc " + i + " is a half turn, which plane is not determined" + + " without a radius vector in " + centre.getDimension() + " dimensions."); + } + } + arcs.add(new Vector<?>[] {start, middle, end}); + } + return arcs; + } + } 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 53eb4a2a11..cceec44f38 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 @@ -16,9 +16,12 @@ */ package org.apache.sis.geometries.internal.shared; +import java.util.List; +import javax.measure.Quantity; import org.opengis.geometry.Envelope; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.curve.ArcByBulge; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometry.GeneralEnvelope; import org.apache.sis.maths.DataType; import org.apache.sis.maths.Vector; @@ -86,11 +89,33 @@ public non-sealed class DefaultArcByBulge extends AbstractGeometry implements Ar @Override public Envelope getEnvelope() { final var envelope = new GeneralEnvelope(points.getEnvelope()); - if (bulge == 0) { - return envelope; // A straight segment. + if (bulge != 0) { + final Vector<?> start = points.getPosition(0).copy(DataType.DOUBLE); + final Vector<?> end = points.getPosition(1).copy(DataType.DOUBLE); + Curves.addArc(envelope, start, apex(start, end), end); } + return envelope; + } + + /** + * Returns a polyline whose points are on this arc, from its start point to its end point. + * + * @throws IllegalStateException if the two end points are identical, or if the normal is parallel to the chord. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { final Vector<?> start = points.getPosition(0).copy(DataType.DOUBLE); final Vector<?> end = points.getPosition(1).copy(DataType.DOUBLE); + final Vector<?> apex = (bulge != 0) ? apex(start, end) : start.copy(DataType.DOUBLE).lerp(end, 0.5); + return CurveLinearizer.linearize(this, spacing, offset, List.of(CurveLinearizer.circularArc(start, apex, end, false))); + } + + /** + * Returns the point of the arc at the bulge distance from the middle of the chord. + * + * @throws IllegalStateException if the two end points are identical, or if the normal is parallel to the chord. + */ + private Vector<?> apex(final Vector<?> start, final Vector<?> end) { 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(); @@ -106,10 +131,7 @@ public non-sealed class DefaultArcByBulge extends AbstractGeometry implements Ar 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; + return start.copy(DataType.DOUBLE).lerp(end, 0.5).add(n.scale(bulge / perpendicular)); } } 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 eb65f8e921..8f1c206938 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,13 +16,18 @@ */ package org.apache.sis.geometries.internal.shared; +import java.util.List; +import java.util.function.DoubleFunction; import javax.measure.IncommensurableException; +import javax.measure.Quantity; 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.geometries.curve.LineString; import org.apache.sis.maths.DataType; +import org.apache.sis.maths.ReadOnly; import org.apache.sis.maths.Vector; @@ -135,6 +140,29 @@ public non-sealed class DefaultArcByCenterPoint extends AbstractGeometry impleme return envelope; } + /** + * Returns a polyline whose points are on this arc, from its start angle to its end angle. + * + * @throws UnsupportedOperationException if the radius unit cannot be converted to the unit of an axis. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + final ReadOnly.Vector<?> c = center.getPosition(); + final double cx = c.get(0); + final double cy = c.get(1); + final double rx = radius(0); + final double ry = radius(1); + final DoubleFunction<Vector<?>> function = (angle) -> { + final double t = Math.toRadians(angle); + final Vector<?> position = c.copy(DataType.DOUBLE); + position.set(0, cx + rx * Math.cos(t)); + position.set(1, cy + ry * Math.sin(t)); + return position; + }; + return CurveLinearizer.linearize(this, spacing, offset, + List.of(new CurveLinearizer.Piece(function, new double[] {startAngle, endAngle}, null))); + } + /** * 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. 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 0d9eb3be22..099c50ba8d 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 @@ -23,10 +23,12 @@ import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.Point; import org.apache.sis.geometries.curve.BSplineCurve; import org.apache.sis.geometries.curve.FunctionArc; import org.apache.sis.geometries.curve.KnotType; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometries.curve.SplineCurveForm; import org.apache.sis.maths.Array; import org.apache.sis.measure.Quantities; @@ -222,4 +224,15 @@ public non-sealed class DefaultBSplineCurve extends AbstractGeometry implements final double[] kn = (knots.length == count + degree + 1) ? knots : BSplineEvaluator.clampedKnots(count, degree, knots); return new BSplineEvaluator(ctrl, null, kn, degree); } + + /** + * Returns a polyline whose points are on this curve, including the positions at the knots. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + if (isEmpty()) { + return GeometryFactory.DEFAULT.createLineString(points); + } + return evaluator().asLine(this, spacing, offset); + } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBezier.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBezier.java index af59ebfc97..62752c90cf 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBezier.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultBezier.java @@ -19,9 +19,11 @@ package org.apache.sis.geometries.internal.shared; import javax.measure.Quantity; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.Point; import org.apache.sis.geometries.curve.Bezier; import org.apache.sis.geometries.curve.KnotType; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometries.curve.SplineCurveForm; import org.apache.sis.maths.Array; import org.apache.sis.maths.Vector; @@ -91,4 +93,14 @@ public non-sealed class DefaultBezier extends DefaultPolynomialSpline implements return new BSplineEvaluator(ctrl, null, BSplineEvaluator.bezierKnots(ctrl.size(), degree), degree); } + /** + * Returns a polyline whose points are on this curve, including the positions at the knots. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + if (isEmpty()) { + return GeometryFactory.DEFAULT.createLineString(points); + } + return evaluator().asLine(this, spacing, offset); + } } 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 7c5fb9dbf1..f4dd587056 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 @@ -16,9 +16,13 @@ */ package org.apache.sis.geometries.internal.shared; +import java.util.ArrayList; +import javax.measure.Quantity; import org.opengis.geometry.Envelope; import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.curve.CircularString; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometry.GeneralEnvelope; import org.apache.sis.maths.DataType; import org.apache.sis.maths.Vector; @@ -72,4 +76,22 @@ public non-sealed class DefaultCircularString extends AbstractGeometry implement return envelope; } + /** + * Returns a polyline whose points are on this circular string, including all its data points. + * + * @throws IllegalStateException if an arc is a full circle in more than two dimensions. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + final var pieces = new ArrayList<CurveLinearizer.Piece>(); + for (int i = 2; i < points.size(); i += 2) { + pieces.add(CurveLinearizer.circularArc(points.getPosition(i - 2).copy(DataType.DOUBLE), + points.getPosition(i - 1).copy(DataType.DOUBLE), points.getPosition(i).copy(DataType.DOUBLE), true)); + } + if (pieces.isEmpty()) { + return GeometryFactory.DEFAULT.createLineString(points); + } + return CurveLinearizer.linearize(this, spacing, offset, pieces); + } + } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCompoundCurve.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCompoundCurve.java index facd052bdb..ee5a1c4ab2 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCompoundCurve.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultCompoundCurve.java @@ -17,11 +17,14 @@ package org.apache.sis.geometries.internal.shared; import java.util.ArrayList; +import java.util.List; +import javax.measure.Quantity; import org.opengis.geometry.Envelope; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.Curve; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.curve.CompoundCurve; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometries.Geometries; import org.apache.sis.maths.Array; import org.apache.sis.maths.DataType; @@ -133,4 +136,12 @@ public non-sealed class DefaultCompoundCurve extends AbstractGeometry implements return NDArrays.concatenate(arrays.toArray(Array[]::new)); } + /** + * Returns a polyline whose points are on this curve, made of the polylines of the component curves. + * A spacing or an offset in percent is relative to the length of this whole curve. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + return CurveLinearizer.concatenate(this, List.of(curves), spacing, offset); + } } 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 ae82eb2284..8b7b5953f2 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 @@ -17,16 +17,19 @@ package org.apache.sis.geometries.internal.shared; import java.util.ArrayList; +import java.util.List; import java.util.Objects; +import javax.measure.Quantity; 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.GeometryFactory; import org.apache.sis.geometries.curve.Conic; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.maths.Array; import org.apache.sis.maths.DataType; import org.apache.sis.maths.Vector; -import org.apache.sis.maths.Vectors; /** @@ -115,239 +118,61 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { */ @Override public Envelope getEnvelope() { - final var envelope = new GeometryEnvelope(points.getEnvelope()); if (points.isEmpty()) { + final var envelope = new GeometryEnvelope(getCoordinateReferenceSystem()); + envelope.setToNaN(); 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)); + final var envelope = new GeometryEnvelope(points.getEnvelope()); + for (final DefaultConicArc arc : arcs()) { + arc.addExtremes(envelope); } 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. + * Returns a polyline whose points are on this conic, including all its data points. + * Each arc is computed in the plane of its five data points, as in {@link #getEnvelope()}. + * + * @throws IllegalStateException for the same reasons as {@link #getEnvelope()}. */ - private static void addArc(final GeometryEnvelope 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."); - } - } - } + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + final var pieces = new ArrayList<CurveLinearizer.Piece>(); + for (final DefaultConicArc arc : arcs()) { + pieces.add(arc.piece()); } - /* - * 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(position); - } - } + if (pieces.isEmpty()) { + return GeometryFactory.DEFAULT.createLineString(points); } + return CurveLinearizer.linearize(this, spacing, offset, pieces); } /** - * Returns the real roots of q₂·s² + q₁·s + q₀ = 0. + * Returns the arcs of this conic, made of five 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. */ - 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]; + private List<DefaultConicArc> arcs() { + if (points.isEmpty()) { + return List.of(); } - final double discriminant = q1*q1 - 4*q2*q0; - if (discriminant < 0) { - return new double[0]; + final var positions = new ArrayList<Vector<?>>(); + for (int i = 0; i < points.size(); i++) { + positions.add(points.getPosition(i).copy(DataType.DOUBLE)); } - 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]; - } - } + if (cycle && positions.size() % 4 == 0) { + positions.add(positions.get(0)); // The first data point closes the curve. } - return det; + 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 + '.'); + } + 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))); + } + return arcs; } - } 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/DefaultConicArc.java similarity index 59% copy from incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConic.java copy to incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultConicArc.java index ae82eb2284..0bc3273ece 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/DefaultConicArc.java @@ -16,145 +16,74 @@ */ 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.maths.Array; +import org.apache.sis.geometries.curve.ConicArc; import org.apache.sis.maths.DataType; import org.apache.sis.maths.Vector; import org.apache.sis.maths.Vectors; /** - * A chain of conic section arcs, each of them determined by five data points and drawn in the - * tangent plane at one control point. + * A conic arc going through five 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. * * @author Johann Sorel (Geomatys) */ -public non-sealed class DefaultConic extends AbstractGeometry implements Conic { - +public final class DefaultConicArc implements ConicArc { /** - * Points lying on the conic, four of them being consumed by each arc plus the one shared with - * the previous arc. + * The five data points of the arc. */ - protected final DataPoints points; + private final Vector<?>[] p; /** - * Centers of the exponential maps in which the arcs are constructed, one per arc. - * May be {@code null} if this conic declares no control point. + * Whether the points are collinear or identical, in which case the arc is the polyline joining them + * and the other fields are not set. */ - protected final Array controlPoints; + private final boolean linear; /** - * Whether this conic closes on itself. + * Orthonormal basis (e1, e2) of the plane of the points, with the centroid as origin. The coordinates + * in that plane are divided by the scale, which is the distance of the farthest point from the origin. */ - protected final boolean cycle; + private Vector<?> origin, e1, e2; + private double scale; /** - * Creates a conic from the given data points and control points. - * - * @param points points lying on the conic. At least five of them, the first one also - * acting as the last one when {@code cycle} is {@code true}. - * @param controlPoints centers of the exponential maps of the arcs, or {@code null} if none. - * @param cycle whether this conic closes on itself. + * Coefficients of the conic A·a² + B·ab + C·b² + D·a + E·b + F = 0 in the plane. */ - public DefaultConic(final DataPoints points, final Array controlPoints, final boolean cycle) { - this.points = Objects.requireNonNull(points); - this.controlPoints = controlPoints; - this.cycle = cycle; - } - - @Override - public DataPoints getDataPoints() { - return points; - } - - @Override - public Array getControlPoints() { - return controlPoints; - } - - @Override - public boolean isCycle() { - return cycle; - } - - @Override - public CoordinateReferenceSystem getCoordinateReferenceSystem() { - return points.getCoordinateReferenceSystem(); - } - - @Override - public void setCoordinateReferenceSystem(CoordinateReferenceSystem crs) throws IllegalArgumentException { - points.setCoordinateReferenceSystem(crs); - } + private double A, B, C, D, E, F; - @Override - public DataPointsType getDataPointsType() { - return points.getType(); - } + /** + * Coordinates in the plane of the start point p₀. + */ + private double a0, b0; - @Override - public boolean isEmpty() { - return points.isEmpty(); - } + /** + * 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]; /** - * 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. + * Whether the directions from p₀ to the next data points increase. */ - @Override - public Envelope getEnvelope() { - final var envelope = new GeometryEnvelope(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; - } + private boolean increasing; /** - * 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. + * Computes the conic through the given points. + * + * @param p the five data points of the arc, in their order along the arc. + * @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. */ - private static void addArc(final GeometryEnvelope envelope, final Vector<?>[] p) { + DefaultConicArc(final Vector<?>[] p) { + this.p = 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); + 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; + double width = 0; for (final Vector<?> q : p) { final Vector<?> v = q.copy(DataType.DOUBLE).subtract(origin); final double length = v.length(); @@ -164,7 +93,8 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { } } if (scale == 0) { - return; // All points are identical. + linear = true; // All points are identical. + return; } e1.scale(1 / scale); for (final Vector<?> q : p) { @@ -176,8 +106,10 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { } } if (!(width > scale * 1E-10)) { - return; // Collinear points: a straight segment. + linear = true; // Collinear points: a straight segment. + return; } + linear = false; 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]; @@ -190,7 +122,7 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { } } /* - * Coefficients of the conic A·a² + B·ab + C·b² + D·a + E·b + F = 0 through the five points, + * 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. */ final Vector<?> coef = Vectors.createDouble(6); @@ -210,7 +142,7 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { 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); + 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."); @@ -220,19 +152,18 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { * 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]; + 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++) { r[k] = relative(direction(a[k] - a0, b[k] - b0), tangent); } - final boolean increasing = r[2] < r[4]; + 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. @@ -244,21 +175,46 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { 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) { + if (onArc(relative(direction(Math.cos(twice / 2), Math.sin(twice / 2)), tangent))) { throw new IllegalStateException("A conic arc goes through a point at infinity," + " from one branch of a hyperbola to the other."); } } } } + } + + /** + * 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]; + } + + /** + * Returns the position of the given coordinates in the plane. + */ + private Vector<?> position(final double xa, final double xb) { + return origin.copy(DataType.DOUBLE) + .add(e1.copy(DataType.DOUBLE).scale(scale * xa)) + .add(e2.copy(DataType.DOUBLE).scale(scale * xb)); + } + + /** + * Expands the given envelope for containing the extremes of this arc along each axis. + * The envelope shall already contain the data points. + */ + void addExtremes(final GeometryEnvelope envelope) { + if (linear) { + return; + } /* * 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++) { + for (int i = 0; i < origin.getDimension(); 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; @@ -277,17 +233,40 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { 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(position); + if (onArc(relative(direction(xa - a0, xb - b0), tangent))) { + envelope.add(position(xa, xb)); } } } } + /** + * Returns this arc as a piece of curve. The parameter is the angle between the tangent at p₀ and the + * direction from p₀ to the position, which is the second intersection of the conic with the line through p₀ + * in that direction. If the points are collinear, the parameter is the index of the data points. + */ + CurveLinearizer.Piece piece() { + 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); + } + final double[] t = new double[5]; + for (int k = 1; k < 5; k++) { + t[k] = increasing ? r[k] : Math.PI - r[k]; + } + return new CurveLinearizer.Piece((angle) -> { + // Q(p₀ + s·d) = s·(∇Q(p₀)·d) + s²·q(d) = 0 since p₀ is on the conic. + final double phi = tangent + (increasing ? angle : -angle); + final double da = Math.cos(phi), db = Math.sin(phi); + final double q = A*da*da + B*da*db + C*db*db; + final double g = (2*A*a0 + B*b0 + D) * da + (B*a0 + 2*C*b0 + E) * db; + final double s = (q != 0) ? -g / q : 0; + return position(a0 + s * da, b0 + s * db); + }, t, p); + } + /** * Returns the real roots of q₂·s² + q₁·s + q₀ = 0. */ @@ -349,5 +328,4 @@ public non-sealed class DefaultConic extends AbstractGeometry implements Conic { } return det; } - } 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 1082193172..6b216cacf7 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,8 +16,11 @@ */ 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; @@ -40,4 +43,21 @@ public non-sealed class DefaultEllipticArc extends DefaultConic implements Ellip super(points, controlPoints, cycle); } + /** + * 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. + */ + @Override + public Envelope getEnvelope() { + throw new UnsupportedOperationException("Not supported yet."); + } + + /** + * Not yet supported, for the same reason as {@link #getEnvelope()}. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + throw new UnsupportedOperationException("Not supported yet."); + } + } 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 90de4c52ed..b74cb7ab43 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 @@ -23,7 +23,9 @@ import org.apache.sis.geometries.BBox; import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.Geometries; +import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.Point; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometries.curve.NurbCurve; import org.apache.sis.geometries.operation.DeBoorAlgorithm; import org.apache.sis.maths.Array; @@ -335,4 +337,15 @@ public non-sealed class DefaultNurbCurve extends AbstractGeometry implements Nur static Vector<?> fromHomogeneous(final Vector<?> h, final int spatialDim) { return h.shrink(spatialDim).scale(1.0 / h.get(spatialDim)); } + + /** + * Returns a polyline whose points are on this curve, including the positions at the knots. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + if (isEmpty()) { + return GeometryFactory.DEFAULT.createLineString(points); + } + return new BSplineEvaluator(points, weights, knots, degree).asLine(this, spacing, offset); + } } 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 1dd945c9ce..bc1d1388c2 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 @@ -28,6 +28,7 @@ import org.apache.sis.geometries.Orientable; import org.apache.sis.geometries.Point; import org.apache.sis.geometries.Primitive; import org.apache.sis.geometries.DataPointsType; +import org.apache.sis.geometries.curve.LineString; import org.apache.sis.maths.Array; import org.apache.sis.measure.Quantities; import org.apache.sis.measure.Units; @@ -154,4 +155,12 @@ public non-sealed class DefaultReversedCurve extends AbstractGeometry implements public DataPoints getDataPoints() { return base.getDataPoints(); } + + /** + * Returns the polyline of the base curve with its points in reverse order. + */ + @Override + public LineString asLine(final Quantity<?> spacing, final Quantity<?> offset) { + return CurveLinearizer.reverse(this, base, spacing, offset); + } } diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/internal/shared/CurveLinearizerTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/internal/shared/CurveLinearizerTest.java new file mode 100644 index 0000000000..42eab9b13d --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/internal/shared/CurveLinearizerTest.java @@ -0,0 +1,266 @@ +/* + * 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 javax.measure.Quantity; +import org.apache.sis.geometries.Curve; +import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.geometries.curve.Arc; +import org.apache.sis.geometries.curve.LineString; +import org.apache.sis.maths.NDArrays; +import org.apache.sis.maths.Vector; +import org.apache.sis.maths.Vectors; +import org.apache.sis.measure.Quantities; +import org.apache.sis.measure.Units; + +// Test dependencies +import static org.junit.jupiter.api.Assertions.*; +import org.junit.jupiter.api.Test; +import org.apache.sis.geometries.GeometryTest; + + +/** + * Tests {@link Curve#asLine(Quantity, Quantity)} on the curves supported by {@link CurveLinearizer}. + * + * @author Johann Sorel (Geomatys) + */ +public class CurveLinearizerTest extends GeometryTest { + + private static final GeometryFactory F = GeometryFactory.DEFAULT; + + /** + * Returns a distance in the unit of the axes of {@link #CRS_2D}. + */ + private static Quantity<?> degrees(final double value) { + return Quantities.create(value, Units.DEGREE); + } + + /** + * Returns the positions of the given line string. + */ + private static double[][] positions(final LineString line) { + final DataPoints points = line.getDataPoints(); + final double[][] positions = new double[points.size()][]; + for (int i = 0; i < positions.length; i++) { + positions[i] = points.getPosition(i).toArrayDouble(); + } + return positions; + } + + private static void assertPosition(final double x, final double y, final double[] position) { + assertEquals(x, position[0], TOLERANCE, "x"); + assertEquals(y, position[1], TOLERANCE, "y"); + } + + /** + * Verifies that the points are on the circle of the given centre and radius, and that the polyline meets + * the given criteria. The arc length and the distance to the arc of each chord are computed exactly. + * + * @return the number of points. + */ + private static int assertOnCircle(final double[][] points, final double cx, final double cy, final double r, + final double spacing, final double offset) + { + for (final double[] p : points) { + assertEquals(r, Math.hypot(p[0] - cx, p[1] - cy), TOLERANCE, "Point not on the circle."); + } + for (int i = 1; i < points.length; i++) { + final double chord = Math.hypot(points[i][0] - points[i-1][0], points[i][1] - points[i-1][1]); + if (spacing > 0) { + assertTrue(2 * r * Math.asin(chord / (2*r)) <= spacing * (1 + 1E-6), "Arc longer than the spacing."); + } + if (offset > 0) { + assertTrue(r - Math.sqrt(r*r - chord*chord/4) <= offset * (1 + 1E-9), "Chord farther than the offset."); + } + } + return points.length; + } + + @Test + public void testArcByCenterPoint() { + final var arc = F.createArcByCenterPoint(F.createPoint(CRS_2D, 10, 20), 5, null, 0, 90); + double[][] points = positions(arc.asLine(null, null)); + assertPosition(15, 20, points[0]); + assertPosition(10, 25, points[points.length - 1]); + assertTrue(assertOnCircle(points, 10, 20, 5, 0, 0) > 5, "Default precision too coarse."); + for (int i = 1; i < points.length; i++) { + assertTrue(points[i][0] < points[i-1][0], "Points not in their order along the arc."); + } + // Offset: the chords are at most at 0.01 from the arc, which needs at least 8 chords for a quarter of circle. + points = positions(arc.asLine(null, degrees(0.01))); + final int count = assertOnCircle(points, 10, 20, 5, 0, 0.01); + assertTrue(count >= 9 && count <= 17, "Unexpected number of points: " + count); + // Spacing: the arcs between points are at most 0.5 long. + points = positions(arc.asLine(degrees(0.5), null)); + assertOnCircle(points, 10, 20, 5, 0.5, 0); + /* + * Spacing slightly less than one eighth of the arc. The sum of the chords underestimates the length of + * each eighth by about 10⁻⁴ of it, so the arc length shall be extrapolated for detecting the excess. + */ + final double eighth = Math.PI * 5 / 2 / 8 * (1 - 5E-5); + assertTrue(assertOnCircle(positions(arc.asLine(degrees(eighth), null)), 10, 20, 5, eighth, 0) > 9); + // Spacing in percent of the arc length. + points = positions(arc.asLine(Quantities.create(10, Units.PERCENT), null)); + assertOnCircle(points, 10, 20, 5, 0.1 * Math.PI * 5 / 2, 0); + // Clockwise arc. + points = positions(F.createArcByCenterPoint(F.createPoint(CRS_2D, 10, 20), 5, null, 90, 0).asLine(null, null)); + assertPosition(10, 25, points[0]); + assertPosition(15, 20, points[points.length - 1]); + } + + @Test + public void testCircularString() { + final var curve = F.createCircularString(F.createDataPoints(CRS_2D, 0,0, 5,5, 10,0, 5,-5, 0,0)); + final double[][] points = positions(curve.asLine(null, degrees(0.05))); + assertOnCircle(points, 5, 0, 5, 0, 0.05); + assertPosition(0, 0, points[0]); + assertPosition(0, 0, points[points.length - 1]); + // The data points are in the polyline, in their order. + int next = 1; + final double[][] data = {{5,5}, {10,0}, {5,-5}}; + for (final double[] p : points) { + if (next <= data.length && Math.abs(p[0] - data[next-1][0]) < TOLERANCE && Math.abs(p[1] - data[next-1][1]) < TOLERANCE) { + next++; + } + } + assertEquals(data.length + 1, next, "Data point missing in the polyline."); + // Full circle given by identical start and end points. + assertOnCircle(positions(F.createCircularString(F.createDataPoints(CRS_2D, 0,0, 10,0, 0,0)).asLine(null, null)), 5, 0, 5, 0, 0); + } + + @Test + public void testArcByBulge() { + final Vector<?> normal = Vectors.createDouble(2).set(new double[] {0, 1}); + final var arc = F.createArcByBulge(F.createDataPoints(CRS_2D, 0,0, 10,0), 2, normal); + final double[][] points = positions(arc.asLine(degrees(1), null)); + assertOnCircle(points, 5, -5.25, 7.25, 1, 0); + assertPosition(0, 0, points[0]); + assertPosition(10, 0, points[points.length - 1]); + for (final double[] p : points) { + assertTrue(p[1] >= -TOLERANCE && p[1] <= 2 + TOLERANCE, "Point outside the arc."); + } + // A null bulge is a straight segment. + assertEquals(2, positions(F.createArcByBulge(F.createDataPoints(CRS_2D, 0,0, 10,0), 0, normal).asLine(null, null)).length); + } + + @Test + public void testArc() { + final Arc arc = F.createArc(F.createDataPoints(CRS_2D, 5,0, 0,5, -5,0), + NDArrays.of(CRS_2D, new double[] {0,0, 0,0}), null, false); + final double[][] points = positions(arc.asLine(null, degrees(0.01))); + assertOnCircle(points, 0, 0, 5, 0, 0.01); + for (final double[] p : points) { + assertTrue(p[1] >= -TOLERANCE, "Point on the wrong side of the circle."); + } + // The envelope of an arc uses its centres, not the five points of a general conic. + final var envelope = arc.getEnvelope(); + assertEquals(-5, envelope.getMinimum(0), TOLERANCE); + assertEquals( 0, envelope.getMinimum(1), TOLERANCE); + assertEquals( 5, envelope.getMaximum(0), TOLERANCE); + assertEquals( 5, envelope.getMaximum(1), TOLERANCE); + } + + @Test + public void testConic() { + final double[] t = {-0.3, 0.4, 1.2, 2.0, 2.6}; + final double[] coordinates = new double[10]; + for (int i = 0; i < t.length; i++) { + coordinates[2*i] = 2 * Math.cos(t[i]); + coordinates[2*i + 1] = Math.sin(t[i]); + } + final var conic = F.createConic(F.createDataPoints(CRS_2D, coordinates), null, false); + final double[][] points = positions(conic.asLine(null, degrees(0.001))); + double previous = Double.NEGATIVE_INFINITY; + for (final double[] p : points) { + 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; + } + assertPosition(coordinates[0], coordinates[1], points[0]); + assertPosition(coordinates[8], coordinates[9], points[points.length - 1]); + assertTrue(points.length > 20, "Offset not met."); + } + + @Test + public void testNurbCurve() { + // A quarter of the unit circle, as a rational quadratic curve. + final var curve = F.createNurbCurve(F.createDataPoints(CRS_2D, 1,0, 1,1, 0,1), + new double[] {1, Math.sqrt(0.5), 1}, new double[] {0, 0, 0, 1, 1, 1}, 2); + final double[][] points = positions(curve.asLine(degrees(0.1), degrees(0.001))); + assertOnCircle(points, 0, 0, 1, 0.1, 0.001); + assertPosition(1, 0, points[0]); + assertPosition(0, 1, points[points.length - 1]); + } + + @Test + public void testBezier() { + // Quadratic Bézier curve on the parabola y = x². + final var curve = F.createBezier(F.createDataPoints(CRS_2D, -1,1, 1,1), + NDArrays.of(CRS_2D, new double[] {-1,1, 0,-1, 1,1}), new double[] {0, 1}, 2, null, null, null, null); + final double[][] points = positions(curve.asLine(null, null)); + assertTrue(points.length > 5); + for (final double[] p : points) { + assertEquals(p[0] * p[0], p[1], TOLERANCE, "Point not on the parabola."); + } + assertPosition(-1, 1, points[0]); + assertPosition( 1, 1, points[points.length - 1]); + } + + @Test + public void testCompoundAndReversed() { + final var line = F.createLineString(F.createDataPoints(CRS_2D, 0,0, 5,0)); + final var arc = F.createCircularString(F.createDataPoints(CRS_2D, 5,0, 10,5, 15,0)); + final var compound = F.createCompoundCurve(line, arc); + double[][] points = positions(compound.asLine(degrees(1), null)); + assertPosition(0, 0, points[0]); + assertPosition(15, 0, points[points.length - 1]); + for (int i = 1; i < points.length; i++) { + final double step = Math.hypot(points[i][0] - points[i-1][0], points[i][1] - points[i-1][1]); + assertTrue(step > TOLERANCE, "Repeated point."); + assertTrue(step <= 1 + TOLERANCE, "Points farther than the spacing."); + } + // The reversed curve has the same points in reverse order. + final double[][] reversed = positions(F.createReversed(compound).asLine(degrees(1), null)); + assertEquals(points.length, reversed.length); + for (int i = 0; i < points.length; i++) { + assertArrayEquals(points[i], reversed[reversed.length - 1 - i], TOLERANCE); + } + } + + @Test + public void testLineString() { + final var line = F.createLineString(F.createDataPoints(CRS_2D, 0,0, 10,0)); + assertSame(line, line.asLine(null, null)); + assertSame(line, line.asLine(null, degrees(1))); + final double[][] points = positions(line.asLine(degrees(3), null)); + assertEquals(5, points.length); + for (int i = 0; i < points.length; i++) { + assertPosition(2.5 * i, 0, points[i]); + } + } + + @Test + public void testInvalid() { + final var arc = F.createArcByCenterPoint(F.createPoint(CRS_2D, 10, 20), 5, null, 0, 90); + assertThrows(IllegalArgumentException.class, () -> arc.asLine(degrees(-1), null)); + final Vector<?> origin = Vectors.createDouble(CRS_2D).set(new double[] {0, 0}); + final Vector<?> direction = Vectors.createDouble(CRS_2D).set(new double[] {1, 0}); + assertThrows(UnsupportedOperationException.class, () -> F.createRay(origin, direction).asLine(null, null)); + } +}
