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 adcc259f4d37d4b7b8c9070dff230cb5d218c77d Author: jsorel <[email protected]> AuthorDate: Tue Oct 13 15:20:06 2026 +0200 feat(Geometry): implement BSpline.subCurve --- .../internal/shared/BSplineEvaluator.java | 379 +++++++++++++++++++++ .../sis/geometries/internal/shared/Curves.java | 43 +-- .../internal/shared/DefaultBSplineCurve.java | 44 +++ .../geometries/internal/shared/DefaultBezier.java | 36 ++ .../internal/shared/DefaultNurbCurve.java | 28 +- .../sis/geometries/curve/BSplineCurveTest.java | 54 ++- .../apache/sis/geometries/curve/BezierTest.java | 37 ++ .../sis/geometries/curve/LineStringTest.java | 33 -- .../apache/sis/geometries/curve/NurbCurveTest.java | 61 ++++ 9 files changed, 638 insertions(+), 77 deletions(-) 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 new file mode 100644 index 0000000000..5a2b2238e2 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/BSplineEvaluator.java @@ -0,0 +1,379 @@ +/* + * 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.Arrays; +import java.util.List; +import javax.measure.Quantity; +import javax.measure.Unit; +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.Point; +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; + + +/** + * Arc length parameterization of a B-spline, optionally rational, in the coordinate space of its + * control points. A Bézier curve is the B-spline of its control points with clamped knots. + * + * <p>The arc length is integrated over each knot span by adaptive Gauss-Legendre quadrature of the + * speed, computed from the analytic derivatives of the basis functions. The parameter at a given + * arc length is then found by Newton iterations safeguarded by bisection.</p> + * + * @author Johann Sorel (Geomatys) + */ +final class BSplineEvaluator { + /** + * The quadrature used for integrating the speed into an arc length. + */ + private static final GaussLegendreQuadrature QUADRATURE = new GaussLegendreQuadrature(10); + + /** + * Relative accuracy of the integrated arc length. + */ + private static final double TOLERANCE = 1E-13; + + /** + * Maximal number of interval halvings when integrating the arc length. + */ + private static final int MAX_DEPTH = 20; + + private final DataPoints controls; + private final double[][] positions; + private final double[] weights; + private final double[] knots; + private final int degree; + + /** + * Index of the last control point. + */ + private final int last; + + /** + * Intervals of the parameter over which the arc length has been integrated, in increasing order. + * Each entry is {start, end, arc length at start, arc length of the interval, knot span index}. + */ + private final List<double[]> intervals = new ArrayList<>(); + + /** + * Total arc length. + */ + private double length; + + /** + * Creates an evaluator for the given B-spline. + * + * @param controls control points, with their attributes. + * @param weights weight of each control point, or {@code null} for a non-rational curve. + * @param knots knot vector, of length {@code controls.size() + degree + 1}. + * @param degree degree of the basis functions. + * @throws IllegalArgumentException if the arguments do not define a B-spline. + */ + BSplineEvaluator(final DataPoints controls, final double[] weights, final double[] knots, final int degree) { + final int count = controls.size(); + if (degree < 0 || count < degree + 1) { + throw new IllegalArgumentException("A curve of degree " + degree + " needs at least " + (degree + 1) + + " control points, but got " + count + '.'); + } + if (knots.length != count + degree + 1) { + throw new IllegalArgumentException("A curve of degree " + degree + " with " + count + " control points needs " + + (count + degree + 1) + " knots, but got " + knots.length + '.'); + } + if (weights != null && weights.length != count) { + throw new IllegalArgumentException("Expected " + count + " weights, but got " + weights.length + '.'); + } + this.controls = controls; + this.knots = knots; + this.degree = degree; + this.last = count - 1; + this.weights = (weights != null) ? weights : filled(count, 1); + positions = new double[count][]; + for (int i = 0; i < count; i++) { + positions[i] = controls.getPosition(i).toArrayDouble(); + } + for (int k = degree; k <= last; k++) { + final double a = knots[k], b = knots[k+1]; + if (b > a) { + final int span = k; + QUADRATURE.integrate((u) -> speed(u, span), a, b, TOLERANCE, MAX_DEPTH, (start, end, integral) -> { + intervals.add(new double[] {start, end, length, integral, span}); + length += integral; + }); + } + } + } + + private static double[] filled(final int count, final double value) { + final double[] array = new double[count]; + Arrays.fill(array, value); + return array; + } + + /** + * Returns a knot vector clamped at both ends, with the given knots in between. + * + * @param count number of control points. + * @param degree degree of the basis functions. + * @param breakpoints the distinct knot values, from the first to the last one, or an empty array + * for uniformly distributed knots between 0 and 1. + * @return the knot vector, of length {@code count + degree + 1}. + * @throws IllegalArgumentException if the number of breakpoints does not match the number of control points. + */ + static double[] clampedKnots(final int count, final int degree, double[] breakpoints) { + final int spans = count - degree; + if (spans < 1) { + throw new IllegalArgumentException("A curve of degree " + degree + " needs at least " + (degree + 1) + + " control points, but got " + count + '.'); + } + if (breakpoints.length == 0) { + breakpoints = new double[spans + 1]; + for (int i = 0; i <= spans; i++) { + breakpoints[i] = i / (double) spans; + } + } else if (breakpoints.length != spans + 1) { + throw new IllegalArgumentException("A clamped curve of degree " + degree + " with " + count + + " control points needs " + (spans + 1) + " distinct knots, but got " + breakpoints.length + '.'); + } + final double[] knots = new double[count + degree + 1]; + Arrays.fill(knots, 0, degree + 1, breakpoints[0]); + System.arraycopy(breakpoints, 1, knots, degree + 1, spans - 1); + Arrays.fill(knots, count, knots.length, breakpoints[spans]); + return knots; + } + + /** + * Returns the knot vector of a piecewise Bézier curve: each piece is a span of the parameter, + * delimited by knots of multiplicity {@code degree}. + * + * @param count number of control points, {@code degree} × pieces + 1. + * @param degree degree of the Bézier pieces. + * @return the knot vector, of length {@code count + degree + 1}. + * @throws IllegalArgumentException if the control points do not make a whole number of pieces. + */ + static double[] bezierKnots(final int count, final int degree) { + if (degree < 1 || count < degree + 1 || (count - 1) % degree != 0) { + throw new IllegalArgumentException("A Bézier curve of degree " + degree + " needs a multiple of " + + degree + " plus one control points, but got " + count + '.'); + } + final int pieces = (count - 1) / degree; + final double[] knots = new double[count + degree + 1]; + for (int i = 0; i < knots.length; i++) { + knots[i] = Math.min(pieces, Math.max(0, (i - 1) / degree)); + } + return knots; + } + + /** + * Returns the given control points as a sequence in the given reference system. + * The control points of some curves are only positions, sometimes without reference system. + * + * @param points the control points. + * @param crs the reference system of the curve. + * @return the control points as a sequence of positions. + */ + static DataPoints asDataPoints(final Array points, final CoordinateReferenceSystem crs) { + if (points.getCoordinateReferenceSystem() != null) { + return GeometryFactory.DEFAULT.createDataPoints(points); + } + final Array copy = NDArrays.of(SampleSystem.of(crs), DataType.DOUBLE, points.getLength()); + for (long i = 0, n = points.getLength(); i < n; i++) { + copy.set(i, points.get(i)); + } + return GeometryFactory.DEFAULT.createDataPoints(copy); + } + + /** + * Returns the total arc length, in the unit of the coordinates. + */ + double length() { + return length; + } + + /** + * Returns the point at the given distance along the curve. + * + * @param dist the distance, in {@code unit} or in percent of the length. + * @param unit unit of the coordinates. + * @return the point at the given distance. + * @throws IllegalArgumentException if the distance is outside the curve. + */ + Point subCurve(final Quantity<?> dist, final Unit<?> unit) { + final double s = Curves.toArcLength(dist, length, unit); + if (intervals.isEmpty()) { + return pointAt(knots[degree], degree); // Degenerate curve with no span. + } + double[] interval = intervals.get(intervals.size() - 1); + for (final double[] candidate : intervals) { + if (candidate[2] + candidate[3] >= s) { + interval = candidate; + break; + } + } + final double a = interval[0], c0 = interval[2], len = interval[3]; + final int k = (int) interval[4]; + double u = a; + if (len > 0) { + double lo = a, hi = interval[1]; + u = a + (hi - a) * Math.min(1, (s - c0) / len); + for (int i = 0; i < 60; i++) { + final double f = c0 + quadrature(k, a, u) - s; + if (Math.abs(f) <= TOLERANCE * length) break; + if (f > 0) hi = u; else lo = u; + final double speed = speed(u, k); + double next = (speed > 0) ? u - f / speed : Double.NaN; + if (!(next > lo && next < hi)) { + next = (lo + hi) / 2; // Bisection when Newton leaves the bracket. + } + if (next == u) break; + u = next; + } + } + return pointAt(u, k); + } + + /** + * Integrates the speed over [a … b] in span {@code k}, which gives the arc length. + */ + private double quadrature(final int k, final double a, final double b) { + return QUADRATURE.integrate((u) -> speed(u, k), a, b); + } + + /** + * Returns the norm of the derivative of the curve at parameter {@code u} in span {@code k}. + */ + private double speed(final double u, final int k) { + final double[] N = new double[degree + 1], dN = new double[degree + 1]; + basis(u, k, N, dN); + final int dim = positions[0].length; + final double[] A = new double[dim], dA = new double[dim]; + double W = 0, dW = 0; + for (int r = 0; r <= degree; r++) { + final int i = k - degree + r; + final double w = weights[i]; + W += N[r] * w; + dW += dN[r] * w; + for (int c = 0; c < dim; c++) { + A[c] += N[r] * w * positions[i][c]; + dA[c] += dN[r] * w * positions[i][c]; + } + } + double sum = 0; + for (int c = 0; c < dim; c++) { + final double d = (dA[c] - dW * A[c] / W) / W; // Quotient rule on A/W. + sum += d * d; + } + return Math.sqrt(sum); + } + + /** + * Computes the {@code degree + 1} basis functions which are non-zero in span {@code k}, + * and optionally their derivatives. + * + * @param u the parameter, in span {@code k}. + * @param k index of the knot span. + * @param N where to store the values of the basis functions {@code k - degree} to {@code k}. + * @param dN where to store their derivatives, or {@code null} if not needed. + */ + private void basis(final double u, final int k, final double[] N, final double[] dN) { + final double[] left = new double[degree + 1], right = new double[degree + 1]; + double[] lower = null; // Basis functions of degree − 1. + N[0] = 1; + for (int j = 1; j <= degree; j++) { + if (j == degree) { + lower = Arrays.copyOf(N, degree); + } + left[j] = u - knots[k + 1 - j]; + right[j] = knots[k + j] - u; + double saved = 0; + for (int r = 0; r < j; r++) { + final double denom = right[r + 1] + left[j - r]; + final double temp = (denom != 0) ? N[r] / denom : 0; + N[r] = saved + right[r + 1] * temp; + saved = left[j - r] * temp; + } + N[j] = saved; + } + if (dN != null) { + for (int r = 0; r <= degree; r++) { + double d = 0; + if (degree > 0) { + final int i = k - degree + r; + if (r >= 1) { + final double denom = knots[i + degree] - knots[i]; + if (denom != 0) d += lower[r - 1] / denom; + } + if (r < degree) { + final double denom = knots[i + degree + 1] - knots[i + 1]; + if (denom != 0) d -= lower[r] / denom; + } + } + dN[r] = degree * d; + } + } + } + + /** + * Returns the point of the curve at parameter {@code u} in span {@code k}, with all the attributes + * of the control points. Attributes interpolated linearly are blended with the same rational basis + * functions as the position. Otherwise the value is the one of the control point having the largest + * weight, or of the first control point of the span for leading values. + */ + private Point pointAt(final double u, final int k) { + final double[] R = new double[degree + 1]; + basis(u, k, R, null); + double W = 0; + int nearest = 0; + for (int r = 0; r <= degree; r++) { + R[r] *= weights[k - degree + r]; + W += R[r]; + } + for (int r = 0; r <= degree; r++) { + R[r] /= W; + if (R[r] > R[nearest]) nearest = r; + } + final int first = k - degree; + final int leading = (u >= knots[last + 1]) ? last : first; + final int dominant = first + nearest; + final DataPointsType type = controls.getType(); + final Point point = GeometryFactory.DEFAULT.createPoint(type); + for (final String name : type.getAttributeNames()) { + final Vector<?> value; + switch (type.getAttributeInterpolation(name)) { + case LINEAR: + case TRANSLATE: { + value = controls.getAttribute(first, name).scale(R[0]); + for (int r = 1; r <= degree; r++) { + value.add(controls.getAttribute(first + r, name).scale(R[r])); + } + break; + } + case LEADING: value = controls.getAttribute(leading, name); break; + default: value = controls.getAttribute(dominant, name); break; + } + point.setAttribute(name, value); + } + return point; + } +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java index 45462b3f3e..aad6862ed1 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/Curves.java @@ -16,19 +16,11 @@ */ package org.apache.sis.geometries.internal.shared; -import java.util.LinkedHashMap; -import java.util.Map; import javax.measure.IncommensurableException; import javax.measure.Quantity; import javax.measure.Unit; -import org.apache.sis.geometries.AttributeInterpolation; 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.maths.Array; -import org.apache.sis.maths.NDArrays; -import org.apache.sis.maths.Vector; import org.apache.sis.measure.Units; @@ -120,41 +112,10 @@ public final class Curves { final double length = segments[i]; if (s <= length || i == segments.length - 1) { final double ratio = (length > 0) ? Math.min(s / length, 1) : 0; - return interpolate(points, i, i+1, ratio); + return points.interpolatePoint(i, i+1, ratio); } s -= length; } - return interpolate(points, 0, 0, 0); // A single position, at distance zero. - } - - /** - * Returns a new point on the segment joining two points of the given sequence. - * The position is linearly interpolated, and each other attribute is computed as declared by - * its {@linkplain DataPointsType#getAttributeInterpolation(String) attribute interpolation}. - * - * @param points the sequence holding the two points. - * @param i0 index of the point at ratio 0, the start of the segment. - * @param i1 index of the point at ratio 1, the end of the segment. - * @param ratio position between the two points, from 0 to 1. - * @return a new point holding the computed attributes, independent of the given sequence. - */ - public static Point interpolate(final DataPoints points, final int i0, final int i1, final double ratio) { - final DataPointsType type = points.getType(); - final Map<String,Array> attributes = new LinkedHashMap<>(); - for (final String name : type.getAttributeNames()) { - final AttributeInterpolation method = DataPointsType.ATT_POSITION.equals(name) - ? AttributeInterpolation.LINEAR : type.getAttributeInterpolation(name); - final Vector<?> value; - switch ((method != null) ? method : AttributeInterpolation.NEAREST) { - case LINEAR: - case TRANSLATE: value = points.getAttribute(i0, name).lerp(points.getAttribute(i1, name), ratio); break; - case LEADING: value = points.getAttribute(ratio < 1 ? i0 : i1, name); break; - default: value = points.getAttribute(ratio <= 0.5 ? i0 : i1, name); break; - } - final Array array = NDArrays.of(type.getAttributeSystem(name), type.getAttributeType(name), 1); - array.set(0, value); - attributes.put(name, array); - } - return GeometryFactory.DEFAULT.createPoint(GeometryFactory.DEFAULT.createDataPoints(attributes)); + return points.interpolatePoint(0, 0, 0); // A single position, at distance zero. } } 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 4c0817756c..661ab93fe1 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 @@ -17,15 +17,19 @@ package org.apache.sis.geometries.internal.shared; 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.Geometries; +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.SplineCurveForm; import org.apache.sis.maths.Array; +import org.apache.sis.measure.Quantities; /** @@ -174,4 +178,44 @@ public non-sealed class DefaultBSplineCurve extends AbstractGeometry implements throw new UnsupportedOperationException("Not supported yet."); } + /** + * Returns the arc length of this curve, integrated in the coordinate space of its control points. + * + * @throws UnsupportedOperationException if this curve is rational. + */ + @Override + public Quantity<?> getLength() { + final double length = isEmpty() ? 0 : evaluator().length(); + return Quantities.create(length, Geometries.getLinearUnit(this)); + } + + /** + * Returns the point at the given distance along this curve, in the unit of the coordinates or in + * percent of the length. + * + * @throws UnsupportedOperationException if this curve is rational. + */ + @Override + public Point subCurve(final Quantity<?> dist) { + if (isEmpty()) { + return null; + } + return evaluator().subCurve(dist, Geometries.getLinearUnit(this)); + } + + /** + * Returns the arc length parameterization of this curve. The control points are the data points + * if no control point is given. The knots are either the complete knot vector, or the distinct + * knot values of a curve clamped at both ends. + */ + private BSplineEvaluator evaluator() { + if (rational) { + throw new UnsupportedOperationException("Rational b-spline curves are not supported, use a NURBS curve."); + } + final DataPoints ctrl = (controlPoints != null) + ? BSplineEvaluator.asDataPoints(controlPoints, getCoordinateReferenceSystem()) : points; + final int count = ctrl.size(); + final double[] kn = (knots.length == count + degree + 1) ? knots : BSplineEvaluator.clampedKnots(count, degree, knots); + return new BSplineEvaluator(ctrl, null, kn, degree); + } } 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 9e88f0d055..af59ebfc97 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 @@ -16,12 +16,16 @@ */ 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.Point; import org.apache.sis.geometries.curve.Bezier; import org.apache.sis.geometries.curve.KnotType; import org.apache.sis.geometries.curve.SplineCurveForm; import org.apache.sis.maths.Array; import org.apache.sis.maths.Vector; +import org.apache.sis.measure.Quantities; /** @@ -55,4 +59,36 @@ public non-sealed class DefaultBezier extends DefaultPolynomialSpline implements derivativeAtStart, derivativeAtEnd, Math.max(0, degree - 1)); } + /** + * Returns the arc length of this curve, integrated in the coordinate space of its control points. + */ + @Override + public Quantity<?> getLength() { + final double length = isEmpty() ? 0 : evaluator().length(); + return Quantities.create(length, Geometries.getLinearUnit(this)); + } + + /** + * Returns the point at the given distance along this curve, in the unit of the coordinates or in + * percent of the length. + */ + @Override + public Point subCurve(final Quantity<?> dist) { + if (isEmpty()) { + return null; + } + return evaluator().subCurve(dist, Geometries.getLinearUnit(this)); + } + + /** + * Returns the arc length parameterization of this curve, as the B-spline whose knots of multiplicity + * {@code degree} delimit the Bézier pieces. The control points are the data points if no control + * point is given. + */ + private BSplineEvaluator evaluator() { + final DataPoints ctrl = (controlPoints != null) + ? BSplineEvaluator.asDataPoints(controlPoints, getCoordinateReferenceSystem()) : points; + return new BSplineEvaluator(ctrl, null, BSplineEvaluator.bezierKnots(ctrl.size(), degree), degree); + } + } 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 4f4b19e07e..c235e5f91a 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 @@ -18,15 +18,19 @@ package org.apache.sis.geometries.internal.shared; import java.util.Arrays; import java.util.List; +import javax.measure.Quantity; 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.Point; import org.apache.sis.geometries.curve.NurbCurve; import org.apache.sis.geometries.operation.DeBoorAlgorithm; import org.apache.sis.maths.Array; import org.apache.sis.maths.NDArrays; import org.apache.sis.maths.Vector; import org.apache.sis.maths.Vectors; +import org.apache.sis.measure.Quantities; /** @@ -85,7 +89,7 @@ public non-sealed class DefaultNurbCurve extends AbstractGeometry implements Nur * Point C(u) of the curve at the given parameter. */ public Vector<?> evaluate(final double u) { - final int spatialDim = controlPoints.length; + final int spatialDim = controlPoints[0].getDimension(); final Vector<?>[] hom = new Vector<?>[nbCtrlPts]; for (int i = 0; i < nbCtrlPts; i++) { hom[i] = toHomogeneous(getControlPoint(i), weights[i]); @@ -137,6 +141,28 @@ public non-sealed class DefaultNurbCurve extends AbstractGeometry implements Nur return Vectors.createDouble(result.length).set(result); } + /** + * Returns the arc length of this curve, integrated in the coordinate space of its control points. + */ + @Override + public Quantity<?> getLength() { + final double length = isEmpty() ? 0 : new BSplineEvaluator(points, weights, knots, degree).length(); + return Quantities.create(length, Geometries.getLinearUnit(this)); + } + + /** + * Returns the point at the given distance along this curve, in the unit of the coordinates or in + * percent of the length. Attributes are blended with the basis functions when interpolated linearly, + * or taken from the control point of largest weight. + */ + @Override + public Point subCurve(final Quantity<?> dist) { + if (isEmpty()) { + return null; + } + return new BSplineEvaluator(points, weights, knots, degree).subCurve(dist, Geometries.getLinearUnit(this)); + } + public double domainStart() { return knots[degree]; } 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 2c1c1a790d..b950a83da9 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 @@ -17,9 +17,12 @@ package org.apache.sis.geometries.curve; import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.measure.Quantities; +import org.apache.sis.measure.Units; // Test dependencies import static org.junit.jupiter.api.Assertions.assertEquals; +import static org.junit.jupiter.api.Assertions.assertThrows; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -84,9 +87,12 @@ public class BSplineCurveTest extends GeometryTest { public void testGetDataPoints() { } + /** + * A b-spline of degree 1 is the polyline of its control points. + */ @Test - @Disabled("Not implemented yet.") public void testGetLength() { + assertEquals(20, createPolyline().getLength().getValue().doubleValue(), 1E-9); } @Test @@ -214,6 +220,50 @@ public class BSplineCurveTest extends GeometryTest { public void testSubCurve_Length_Length() { } + /** + * Creates a cubic b-spline clamped at both ends without interior knot, which is the Bézier curve + * of its 4 control points. The control polygon is symmetric around x = 5, so is the curve. + */ + private static BSplineCurve createArch() { + return GeometryFactory.DEFAULT.createBSplineCurve( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 0,10, 10,10, 10,0), null, null, 3, + null, null, false); + } + + /** + * Creates a b-spline of degree 1 given by its distinct knot values, + * which is the polyline (0 0) (10 0) (10 10) of its control points. + */ + private static BSplineCurve createPolyline() { + return GeometryFactory.DEFAULT.createBSplineCurve( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 10,0, 10,10), null, new double[] {0, 1, 2}, 1, + null, null, false); + } + + @Test + public void testSubCurve_Length() { + final BSplineCurve arch = createArch(); + assertPositionEquals( 0, 0, arch.subCurve(Quantities.create( 0, Units.PERCENT))); + assertPositionEquals( 5, 7.5, arch.subCurve(Quantities.create( 50, Units.PERCENT))); + assertPositionEquals(10, 0, arch.subCurve(Quantities.create(100, Units.PERCENT))); + + final BSplineCurve polyline = createPolyline(); + assertPositionEquals( 5, 0, polyline.subCurve(Quantities.create(25, Units.PERCENT))); + assertPositionEquals(10, 5, polyline.subCurve(Quantities.create(75, Units.PERCENT))); + assertPositionEquals(10, 8, polyline.subCurve(Quantities.create(18, Units.DEGREE))); + } + + /** + * Rational b-splines are not supported: NURBS curves shall be used instead. + */ + @Test + public void testSubCurve_LengthRational() { + final BSplineCurve curve = GeometryFactory.DEFAULT.createBSplineCurve( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 0,10, 10,10, 10,0), null, null, 3, + null, null, true); + assertThrows(UnsupportedOperationException.class, () -> curve.subCurve(Quantities.create(50, Units.PERCENT))); + } + @Test @Disabled("Not implemented yet.") public void testGetOrientationSign() { @@ -300,8 +350,8 @@ public class BSplineCurveTest extends GeometryTest { } @Test - @Disabled("Not implemented yet.") public void testGetRepresentativePoint() { + assertPositionEquals(5, 7.5, createArch().getRepresentativePoint()); } @Test diff --git a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BezierTest.java b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BezierTest.java index df9521bea2..3ae6355e78 100644 --- a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BezierTest.java +++ b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/curve/BezierTest.java @@ -16,7 +16,13 @@ */ package org.apache.sis.geometries.curve; +import org.apache.sis.geometries.GeometryFactory; +import org.apache.sis.maths.NDArrays; +import org.apache.sis.measure.Quantities; +import org.apache.sis.measure.Units; + // Test dependencies +import static org.junit.jupiter.api.Assertions.assertThrows; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -32,4 +38,35 @@ public class BezierTest extends GeometryTest { @Disabled("Not implemented yet.") public void testGetInterpolation() { } + + /** + * Two cubic pieces joined at (10 0): an arch above the x axis from (0 0), then its point + * reflection through (10 0) below the axis up to (20 0). Each piece is symmetric, so its + * middle is at parameter ½, and both have the same length, so the curve middle is the junction. + */ + @Test + public void testSubCurve_Length() { + final Bezier curve = GeometryFactory.DEFAULT.createBezier( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 10,0, 20,0), + NDArrays.of(CRS_2D, new double[] {0,0, 0,10, 10,10, 10,0, 10,-10, 20,-10, 20,0}), + new double[] {0, 1}, 3, null, null, null, null); + assertPositionEquals( 0, 0, curve.subCurve(Quantities.create( 0, Units.PERCENT))); + assertPositionEquals( 5, 7.5, curve.subCurve(Quantities.create( 25, Units.PERCENT))); + assertPositionEquals(10, 0, curve.subCurve(Quantities.create( 50, Units.PERCENT))); + assertPositionEquals(15, -7.5, curve.subCurve(Quantities.create( 75, Units.PERCENT))); + assertPositionEquals(20, 0, curve.subCurve(Quantities.create(100, Units.PERCENT))); + assertPositionEquals(10, 0, curve.getRepresentativePoint()); + } + + /** + * The control points shall make a whole number of pieces of the curve degree. + */ + @Test + public void testSubCurve_LengthInvalid() { + final Bezier curve = GeometryFactory.DEFAULT.createBezier( + GeometryFactory.DEFAULT.createDataPoints(CRS_2D, 0,0, 10,0), + NDArrays.of(CRS_2D, new double[] {0,0, 0,10, 10,10, 10,5, 10,0}), + new double[] {0, 1}, 3, null, null, null, null); + assertThrows(IllegalArgumentException.class, () -> curve.subCurve(Quantities.create(50, Units.PERCENT))); + } } 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 7146cab4d7..1e8cd52b46 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 @@ -17,24 +17,18 @@ package org.apache.sis.geometries.curve; import java.util.LinkedHashMap; -import java.util.List; import java.util.Map; import java.util.Set; import javax.measure.Quantity; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.geometries.AttributeInterpolation; -import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.DataPointsType; import org.apache.sis.geometries.Geometry; import org.apache.sis.geometries.GeometryCollection; import org.apache.sis.geometries.GeometryFactory; import org.apache.sis.geometries.Point; 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; import org.apache.sis.measure.Quantities; import org.apache.sis.measure.Units; @@ -326,33 +320,6 @@ public abstract class LineStringTest extends GeometryTest { assertEquals(30, point.getAttribute("leading").get(0), TOLERANCE, "The end point has its own value."); } - /** - * Returns a view of the given sequence declaring the given attribute interpolations. - * Attributes absent from the map keep the interpolation declared by the sequence. - */ - private static DataPoints withInterpolations(final DataPoints points, final Map<String,AttributeInterpolation> methods) { - final DataPointsType base = points.getType(); - final DataPointsType type = new DataPointsType() { - @Override public List<String> getAttributeNames() {return base.getAttributeNames();} - @Override public SampleSystem getAttributeSystem(String name) {return base.getAttributeSystem(name);} - @Override public DataType getAttributeType(String name) {return base.getAttributeType(name);} - @Override public AttributeInterpolation getAttributeInterpolation(String name) { - return methods.getOrDefault(name, base.getAttributeInterpolation(name)); - } - }; - return new DataPoints() { - @Override public CoordinateReferenceSystem getCoordinateReferenceSystem() {return points.getCoordinateReferenceSystem();} - @Override public void setCoordinateReferenceSystem(CoordinateReferenceSystem crs) {points.setCoordinateReferenceSystem(crs);} - @Override public DataPointsType getType() {return type;} - @Override public int size() {return points.size();} - @Override public Point getPoint(int index) {return points.getPoint(index);} - @Override public Vector<?> getPosition(int index) {return points.getPosition(index);} - @Override public void setPosition(int index, ReadOnly.Vector<?> value) {points.setPosition(index, value);} - @Override public Vector<?> getAttribute(int index, String name) {return points.getAttribute(index, name);} - @Override public void setAttribute(int index, String name, ReadOnly.Vector<?> value) {points.setAttribute(index, name, value);} - }; - } - /** * A distance outside the curve, or which cannot be expressed in the unit of the axes, is rejected. */ 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 ae52a6b4d5..1ecd09cc7b 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 @@ -16,7 +16,20 @@ */ package org.apache.sis.geometries.curve; +import java.util.LinkedHashMap; +import java.util.Map; +import org.apache.sis.geometries.AttributeInterpolation; +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.maths.Array; +import org.apache.sis.maths.NDArrays; +import org.apache.sis.measure.Quantities; +import org.apache.sis.measure.Units; + // Test dependencies +import static org.junit.jupiter.api.Assertions.assertEquals; import org.junit.jupiter.api.Disabled; import org.junit.jupiter.api.Test; import org.apache.sis.geometries.GeometryTest; @@ -82,4 +95,52 @@ public class NurbCurveTest extends GeometryTest { @Disabled("Not implemented yet.") public void testGetSegment() { } + + /** + * Creates the quarter of the circle of radius 10 centered on (0 0), from (10 0) to (0 10), + * as a rational quadratic curve. Its length is 5π. + */ + private static NurbCurve createQuarterCircle(final DataPoints points) { + return GeometryFactory.DEFAULT.createNurbCurve(points, + new double[] {1, Math.sqrt(0.5), 1}, new double[] {0,0,0, 1,1,1}, 2); + } + + @Test + public void testGetLength() { + final NurbCurve arc = createQuarterCircle(GeometryFactory.DEFAULT.createDataPoints(CRS_2D, new double[] {10,0, 10,10, 0,10})); + assertEquals(5 * Math.PI, arc.getLength().getValue().doubleValue(), 1E-9); + } + + /** + * The point at a distance along the quarter circle is at the angle of that distance divided by the radius. + * The rational parameter is not proportional to the angle, so this verifies the arc length parameterization. + */ + @Test + public void testSubCurve_Length() { + final NurbCurve arc = createQuarterCircle(GeometryFactory.DEFAULT.createDataPoints(CRS_2D, new double[] {10,0, 10,10, 0,10})); + final double h = 10 * Math.sqrt(0.5); + assertPositionEquals(10, 0, arc.subCurve(Quantities.create(0, Units.PERCENT))); + assertPositionEquals( h, h, arc.subCurve(Quantities.create(50, Units.PERCENT))); + assertPositionEquals( 0, 10, arc.subCurve(Quantities.create(100, Units.PERCENT))); + assertPositionEquals(10 * Math.cos(Math.PI / 6), 5, arc.subCurve(Quantities.create(5 * Math.PI / 3, Units.DEGREE))); + assertPositionEquals(5, 10 * Math.sin(Math.PI / 3), arc.subCurve(Quantities.create(10 * Math.PI / 3, Units.DEGREE))); + assertPositionEquals( h, h, arc.getRepresentativePoint()); + } + + /** + * Attributes interpolated linearly are blended with the rational basis functions, + * while the others take the value of the control point of largest weight. + */ + @Test + public void testSubCurve_LengthWithAttributes() { + final Map<String,Array> attributes = new LinkedHashMap<>(); + attributes.put(DataPointsType.ATT_POSITION, NDArrays.of(CRS_2D, new double[] {10,0, 10,10, 0,10})); + attributes.put("nearest", NDArrays.of(1, new double[] {0, 0, 30})); + attributes.put("linear", NDArrays.of(1, new double[] {0, 0, 30})); + final NurbCurve arc = createQuarterCircle(withInterpolations(GeometryFactory.DEFAULT.createDataPoints(attributes), + Map.of("linear", AttributeInterpolation.LINEAR))); + final Point point = arc.subCurve(Quantities.create(50, Units.PERCENT)); + assertEquals(0, point.getAttribute("nearest").get(0), TOLERANCE, "Value of the control point of largest weight."); + assertEquals(30 * (1 - Math.sqrt(0.5)), point.getAttribute("linear").get(0), TOLERANCE, "Blended by the basis functions."); + } }
