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 85a2bc74776e67066a5f7f7728dd3a03c0f955dd Author: jsorel <[email protected]> AuthorDate: Tue Sep 8 10:26:07 2026 +0200 feat(Geometry): add an incomplete implementation of nurbs curve and surface --- .../org/apache/sis/geometries/GeometryType.java | 17 +- .../main/org/apache/sis/geometries/Surface.java | 2 +- .../org/apache/sis/geometries/curve/NurbCurve.java | 101 ++++++ .../internal/shared/DefaultNurbCurve.java | 311 +++++++++++++++++ .../internal/shared/DefaultNurbSurface.java | 383 +++++++++++++++++++++ .../sis/geometries/operation/Intersection.java | 362 ++++++++++++++++++- .../geometries/operation/SutherlandHodgman.java | 12 +- .../org/apache/sis/geometries/solid/Sphere.java | 3 +- .../sis/geometries/surface/BSplineSurface.java | 29 +- .../{curve/NURB.java => surface/NurbSurface.java} | 18 +- .../geometries/surface/ParametricCurveSurface.java | 3 +- .../main/org/apache/sis/maths/Array.java | 17 + .../main/org/apache/sis/maths/NDArrays.java | 4 + 13 files changed, 1232 insertions(+), 30 deletions(-) diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/GeometryType.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/GeometryType.java index e9b4d1e07a..8962c32e57 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/GeometryType.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/GeometryType.java @@ -21,6 +21,7 @@ import org.opengis.annotation.UML; /** + * NOTE : ISO:19107 define a short list, but in different part of the UML it refers to more accurate types. * * @author Johann Sorel (Geomatys) */ @@ -28,15 +29,23 @@ import org.opengis.annotation.UML; public enum GeometryType { EMPTY, GEOMETRY, + COLLECTION, + + //point types POINT, + + //curve types CURVE, - SURFACE, LINE, GEODESIC, - POLYGON, - COLLECTION, + RHUMB, SPLINECURVE, + + //surface types + SURFACE, + POLYGON, SPLINESURFACE, + + //solid types SPLINESOLID, - RHUMB } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Surface.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Surface.java index eb50d39ed9..3a4f306a2b 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Surface.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/Surface.java @@ -116,7 +116,7 @@ public interface Surface extends Orientable { } @UML(identifier="dataPoint", specification=ISO_19107) // section 6.4.25.8 - default List<DirectPosition> getDataPoints() { + default DataPoints getDataPoints() { //TODO throw new UnsupportedOperationException(); } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/NurbCurve.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/NurbCurve.java new file mode 100644 index 0000000000..33004cd862 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/NurbCurve.java @@ -0,0 +1,101 @@ +/* + * 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.curve; + +import org.apache.sis.geometries.AttributesType; +import org.apache.sis.geometries.CurveInterpolation; +import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.internal.shared.AbstractGeometry; +import static org.opengis.annotation.Specification.ISO_19107; +import org.opengis.annotation.UML; +import org.opengis.referencing.crs.CoordinateReferenceSystem; + + +/** + * + * @author Johann Sorel (Geomatys) + */ +@UML(identifier="NURB", specification=ISO_19107) // section 7.13.8 +public interface NurbCurve extends BSplineCurve { + + public static final String TYPE = "NURBS"; + + @Override + default String getGeometryType() { + return TYPE; + } + + @UML(identifier="interpolation", specification=ISO_19107) // section 7.1.2.2 + @Override + default CurveInterpolation getInterpolation() { + return CurveInterpolation.NURBS; + } + + @Override + default CoordinateReferenceSystem getCoordinateReferenceSystem() { + return getDataPoints().getCoordinateReferenceSystem(); + } + + @Override + default void setCoordinateReferenceSystem(CoordinateReferenceSystem crs) throws IllegalArgumentException { + getDataPoints().setCoordinateReferenceSystem(crs); + } + + @Override + default AttributesType getAttributesType() { + return getDataPoints().getAttributesType(); + } + + @Override + default boolean isEmpty() { + return getDataPoints().isEmpty(); + } + + @Override + default SplineCurveForm getCurveForm() { + return null; + } + + @Override + default KnotType getKnotSpec() { + return KnotType.NON_UNIFORM; + } + + @Override + default boolean isRational() { + return true; + } + + @Override + default Integer getNumArc() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + default FunctionArc getSegment(int idx) { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + default String asText() { + final StringBuilder sb = new StringBuilder("NURBS ("); + final DataPoints points = getDataPoints(); + AbstractGeometry.toText(sb, points); + sb.append(')'); + return sb.toString(); + } +} 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 new file mode 100644 index 0000000000..8267e1edb2 --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbCurve.java @@ -0,0 +1,311 @@ +/* + * 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.Arrays; +import java.util.List; +import org.apache.sis.geometries.AttributesType; +import org.apache.sis.geometries.BBox; +import org.apache.sis.geometries.DataPoints; +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; + +/** + * A curve defined by control points, weights, a knot vector and a degree. + * + * Depending on the given values it is a Bézier curve, a B-spline or a NURBS: + * the weights make it rational, and the knot vector multiplicities make it clamped or periodic. + * + * @author Johann Sorel (Geomatys) + */ +public final class DefaultNurbCurve extends AbstractGeometry implements NurbCurve { + + final DataPoints points; + final Array controlPointsArray; + final Vector<?>[] controlPoints; + final int nbCtrlPts; + final double[] weights; + final double[] knots; + final int degree; + + public DefaultNurbCurve(final DataPoints points, final double[] weights, final double[] knots, final int degree) { + this.points = points; + this.controlPointsArray = points.getAttributeArray(AttributesType.ATT_POSITION); + this.controlPoints = controlPointsArray.toArray(); + this.weights = weights; + this.knots = knots; + this.degree = degree; + this.nbCtrlPts = (int) this.controlPointsArray.getLength(); + } + + @Override + public int getDegree() { + return degree; + } + + @Override + public double[] getKnots() { + return knots.clone(); + } + + @Override + public DataPoints getDataPoints() { + return points; + } + + @Override + public Array getControlPoints() { + return controlPointsArray; + } + + Vector getControlPoint(int index) { + return Vectors.castOrWrap(controlPoints[index]); + } + + /** + * Point C(u) of the curve at the given parameter. + */ + public Vector<?> evaluate(final double u) { + final int spatialDim = controlPoints.length; + final Vector<?>[] hom = new Vector<?>[nbCtrlPts]; + for (int i = 0; i < nbCtrlPts; i++) { + hom[i] = toHomogeneous(getControlPoint(i), weights[i]); + } + final Vector<?> result = DeBoorAlgorithm.evaluate(u, hom, knots, degree); + return fromHomogeneous(result, spatialDim); + } + + /** + * First derivative of a curve (Bézier, B-spline or NURBS) at parameter u. + * + * Principle: the homogeneous curve Cw(u) = Σ N_i,p(u) Pw_i (weighted control points, dimension N+1) is an ordinary + * NON rational B-spline. Its derivative Cw'(u) is itself a B-spline of degree p-1 whose control points are computed + * directly (the classical knot-weighted finite difference formula, without resorting to numerical differentiation). + * + * The quotient rule is then applied in order to "divide" by the weight: + * + * C(u) = A(u) / w(u) C'(u) = (A'(u) - w'(u) C(u)) / w(u) + * + * where A(u) is the spatial part of Cw(u) and w(u) its last component (the homogeneous weight). + * If all weights are 1 (Bézier / non rational B-spline), then w'(u)=0 and we simply fall back on C'(u) = A'(u). + */ + public Vector<?> derivative(final double u) { + 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]); + } + + final Vector<?> Cw = DeBoorAlgorithm.evaluate(u, hom, knots, degree); + final double w = Cw.get(spatialDim); + final double[] Cu = new double[spatialDim]; + for (int c = 0; c < spatialDim; c++) { + Cu[c] = Cw.get(c) / w; + } + + if (degree == 0) { + return Vectors.createDouble(spatialDim); // constant curve -> null derivative + } + final Vector<?>[] derivCtrl = DeBoorAlgorithm.derivativeControlPoints(hom, knots, degree); + final double[] derivKnots = Arrays.copyOfRange(knots, 1, knots.length - 1); + final Vector<?> Cwp = DeBoorAlgorithm.evaluate(u, derivCtrl, derivKnots, degree - 1); + + final double wp = Cwp.get(spatialDim); + final double[] result = new double[spatialDim]; + for (int c = 0; c < spatialDim; c++) { + result[c] = (Cwp.get(c) - wp * Cu[c]) / w; + } + return Vectors.createDouble(result.length).set(result); + } + + public double domainStart() { + return knots[degree]; + } + + public double domainEnd() { + return knots[nbCtrlPts]; + } + + // ------------------------------------------------------------------ + // Subdivision (Boehm's algorithm: knot insertion) + // ------------------------------------------------------------------ + /** + * Returns an equivalent curve, of the same degree and describing exactly the same shape, in which the given knot + * has been inserted the given number of times. + * + * @param u knot value to insert + * @param times how many times to insert it + * @return the refined curve + */ + public DefaultNurbCurve insertKnot(final double u, final int times) { + Vector<?>[] ctrl = controlPoints; + double[] w8 = weights; + double[] kn = knots; + final int p = degree; + final int spatialDim = controlPoints[0].getDimension(); + + for (int t = 0; t < times; t++) { + final int n = ctrl.length - 1; + final int k = DeBoorAlgorithm.findKnotSpan(u, kn, p, n); + final int s = knotMultiplicity(kn, u); + + final Vector<?>[] Pw = new Vector<?>[ctrl.length]; + for (int i = 0; i < ctrl.length; i++) { + Pw[i] = toHomogeneous(ctrl[i], w8[i]); + } + + final Vector<?>[] Qw = new Vector<?>[Pw.length + 1]; + for (int i = 0; i <= k - p; i++) { + Qw[i] = Pw[i]; + } + for (int i = k - p + 1; i <= k - s; i++) { + final double alpha = (u - kn[i]) / (kn[i + p] - kn[i]); + Qw[i] = Pw[i - 1].copy().lerp(Pw[i], alpha); + } + for (int i = k - s; i <= n; i++) { + Qw[i + 1] = Pw[i]; + } + + final double[] newKnots = new double[kn.length + 1]; + System.arraycopy(kn, 0, newKnots, 0, k + 1); + newKnots[k + 1] = u; + System.arraycopy(kn, k + 1, newKnots, k + 2, kn.length - (k + 1)); + + final Vector<?>[] newCtrl = new Vector<?>[Qw.length]; + final double[] newWeights = new double[Qw.length]; + for (int i = 0; i < Qw.length; i++) { + final double w = Qw[i].get(spatialDim); + newWeights[i] = w; + final double[] arr = new double[spatialDim]; + for (int c = 0; c < spatialDim; c++) { + arr[c] = Qw[i].get(c) / w; + } + newCtrl[i] = Vectors.createDouble(arr.length).set(arr); + } + + ctrl = newCtrl; + w8 = newWeights; + kn = newKnots; + } + + final DataPoints ctrlDp = new ArrayDataPoints( + NDArrays.of(List.of(ctrl), + controlPointsArray.getSampleSystem(), + controlPointsArray.getDataType())); + return new DefaultNurbCurve(ctrlDp, w8, kn, p); + } + + /** + * Splits this curve in two halves at the given parameter. + * The knot is first inserted enough times for its multiplicity to reach the degree, + * after which the control points can simply be shared between the two halves. + * + * @param u parameter at which to split the curve + * @return the two halves, in parameter order + */ + public DefaultNurbCurve[] subdivide(final double u) { + final int p = degree; + final int s = knotMultiplicity(knots, u); + + final DefaultNurbCurve refined = (s < p) ? insertKnot(u, p - s) : this; + + final double[] kn = refined.knots; + final int lastIdx = lastIndexOf(kn, u); + final int a = lastIdx - p; + + final Vector<?>[] leftCtrl = Arrays.copyOfRange(refined.controlPoints, 0, a + 1); + final double[] leftWeights = Arrays.copyOfRange(refined.weights, 0, a + 1); + final double[] leftKnots = new double[lastIdx + 2]; + System.arraycopy(kn, 0, leftKnots, 0, lastIdx + 1); + leftKnots[lastIdx + 1] = u; + + final int n = refined.controlPoints.length - 1; + final Vector<?>[] rightCtrl = Arrays.copyOfRange(refined.controlPoints, a, n + 1); + final double[] rightWeights = Arrays.copyOfRange(refined.weights, a, n + 1); + final int rightStart = a + 1; + final double[] rightKnots = new double[kn.length - rightStart + 1]; + rightKnots[0] = u; + System.arraycopy(kn, rightStart, rightKnots, 1, kn.length - rightStart); + + final DataPoints leftDp = new ArrayDataPoints( + NDArrays.of(List.of(leftCtrl), + controlPointsArray.getSampleSystem(), + controlPointsArray.getDataType())); + final DataPoints rightDp = new ArrayDataPoints( + NDArrays.of(List.of(rightCtrl), + controlPointsArray.getSampleSystem(), + controlPointsArray.getDataType())); + + return new DefaultNurbCurve[]{ + new DefaultNurbCurve(leftDp, leftWeights, leftKnots, p), + new DefaultNurbCurve(rightDp, rightWeights, rightKnots, p) + }; + } + + /** + * Bounding box of the control points, which contains the curve (convex hull property). + * + * @return the bounding box of this curve + */ + @Override + public BBox getEnvelope() { + final BBox bbox = new BBox(controlPoints[0], controlPoints[0]); + for (final Vector<?> p : controlPoints) { + bbox.add(p); + } + return bbox; + } + + private static int knotMultiplicity(final double[] knots, final double u) { + int count = 0; + for (final double k : knots) { + if (k == u) { + count++; + } + } + return count; + } + + private static int lastIndexOf(final double[] knots, final double u) { + for (int i = knots.length - 1; i >= 0; i--) { + if (knots[i] == u) { + return i; + } + } + throw new IllegalArgumentException("u=" + u + " is not a knot of the vector"); + } + + // ------------------------------------------------------------------ + // Homogeneous coordinates + // + // A NURBS is evaluated as a non rational B-spline on the weighted points + // (x*w, …, w), which is what makes DeBoorAlgorithm usable as-is; the + // spatial point is recovered by dividing by that last component. + // ------------------------------------------------------------------ + static Vector<?> toHomogeneous(final Vector<?> p, final double w) { + return p.extend(1).scale(w); + } + + static Vector<?> fromHomogeneous(final Vector<?> h, final int spatialDim) { + return h.shrink(spatialDim).scale(1.0 / h.get(spatialDim)); + } +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbSurface.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbSurface.java new file mode 100644 index 0000000000..253d1498fa --- /dev/null +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/internal/shared/DefaultNurbSurface.java @@ -0,0 +1,383 @@ +/* + * 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.Arrays; +import java.util.List; +import org.apache.sis.geometries.BBox; +import org.apache.sis.geometries.Curve; +import org.apache.sis.geometries.DataPoints; +import org.apache.sis.geometries.curve.KnotType; +import org.apache.sis.geometries.operation.DeBoorAlgorithm; +import org.apache.sis.geometries.surface.BSplineSurfaceForm; +import org.apache.sis.geometries.surface.NurbSurface; +import org.apache.sis.maths.NDArrays; +import org.apache.sis.maths.Vector; +import org.apache.sis.maths.Vectors; +import org.opengis.geometry.DirectPosition; +import org.opengis.metadata.Identifier; + +/** + * NURBS surface: tensor product of two directions (u, v), each one with its own degree and its own knot vector. + * The control points form a grid controlPoints[i][j] (i = u direction, j = v direction). + * + * @author Johann Sorel (Geomatys) + */ +public final class DefaultNurbSurface extends AbstractGeometry implements NurbSurface { + + public final Vector<?>[][] controlPoints; + public final double[][] weights; + public final double[] knotsU; + public final double[] knotsV; + public final int degree; + + public DefaultNurbSurface(final Vector<?>[][] controlPoints, final double[][] weights, + final double[] knotsU, final double[] knotsV, final int degree) { + this.controlPoints = controlPoints; + this.weights = weights; + this.knotsU = knotsU; + this.knotsV = knotsV; + this.degree = degree; + } + + /** + * Evaluates the point S(u,v) of the surface at the given parameters. + * + * Principle: the formula S(u,v) = ΣΣ N_i,p(u) N_j,q(v) Pw_ij is computed as two passes of the "classical" De Boor + * algorithm in homogeneous coordinates: 1) for each row i, evaluate a curve along v -> one homogeneous point Q_i + * per row (this "flattens" the v direction) 2) then evaluate a curve along u using the Q_i as homogeneous control + * points and only dehomogenize at the very end — exactly the same principle as for curves, just applied twice. + * + * N_i,p(u) and N_j,q(v) vanish outside of their support, so only degreeU+1 rows and, within each of them, degreeV+1 + * columns actually contribute. They are selected first through {@link DeBoorAlgorithm#findKnotSpan}, then De Boor's + * recursion is run on that small window only, through {@link DeBoorAlgorithm#evaluateWindow} — instead of walking + * the whole grid. The cost is therefore O(degreeU × degreeV) instead of O(rows × cols), independently of the grid + * size. + */ + public Vector<?> evaluate(final double u, final double v) { + final int rows = controlPoints.length; + final int spatialDim = controlPoints[0][0].getDimension(); + + final int ku = DeBoorAlgorithm.findKnotSpan(u, knotsU, degree, rows - 1); + + final Vector<?>[] windowU = new Vector<?>[degree + 1]; + for (int t = 0; t <= degree; t++) { + final int i = ku - degree + t; + windowU[t] = evaluateRowAlongV(controlPoints[i], weights[i], v, knotsV, degree); + } + + final Vector<?> result = DeBoorAlgorithm.evaluateWindow(u, windowU, knotsU, degree, ku); + return DefaultNurbCurve.fromHomogeneous(result, spatialDim); + } + + /** + * Evaluates, for a single row of the grid, the homogeneous point along v (degreeV+1 window only). + */ + private static Vector<?> evaluateRowAlongV(final Vector<?>[] rowControlPoints, final double[] rowWeights, + final double v, final double[] knotsV, final int degreeV) { + final int cols = rowControlPoints.length; + final int kv = DeBoorAlgorithm.findKnotSpan(v, knotsV, degreeV, cols - 1); + + final Vector<?>[] windowV = new Vector<?>[degreeV + 1]; + for (int t = 0; t <= degreeV; t++) { + final int j = kv - degreeV + t; + windowV[t] = DefaultNurbCurve.toHomogeneous(rowControlPoints[j], rowWeights[j]); + } + return DeBoorAlgorithm.evaluateWindow(v, windowV, knotsV, degreeV, kv); + } + + /** + * Partial derivatives {Su, Sv} at point (u,v) — the two vectors tangent to the surface. + * + * Same principle as for curves (differentiation of the homogeneous surface then the quotient rule), applied in each + * direction: + * + * - For Su: the v direction is flattened first (as in {@link #evaluate}), which yields a homogeneous curve in u; it + * is then differentiated with the same formula as for curves. - For Sv: symmetrically, u is flattened first, then + * the differentiation happens in v. + * + * Unlike {@link #evaluate} (optimized to a local window), this method rebuilds the complete flattened curve in each + * direction, because the differentiation formula needs the whole knot vector in order to stay correct at the window + * boundaries. + */ + public Vector<?>[] derivative(final double u, final double v) { + final int rows = controlPoints.length; + final int cols = controlPoints[0].length; + final int spatialDim = controlPoints[0][0].getDimension(); + + // Flatten v -> a homogeneous curve along u (one value per row) + final Vector<?>[] alongU = new Vector<?>[rows]; + for (int i = 0; i < rows; i++) { + final Vector<?>[] rowHom = new Vector<?>[cols]; + for (int j = 0; j < cols; j++) { + rowHom[j] = DefaultNurbCurve.toHomogeneous(controlPoints[i][j], weights[i][j]); + } + alongU[i] = DeBoorAlgorithm.evaluate(v, rowHom, knotsV, degree); + } + final Vector<?> Sw = DeBoorAlgorithm.evaluate(u, alongU, knotsU, degree); + final double w = Sw.get(spatialDim); + final double[] Suv = new double[spatialDim]; + for (int c = 0; c < spatialDim; c++) { + Suv[c] = Sw.get(c) / w; + } + + final Vector<?> su = partialDerivative(u, alongU, knotsU, degree, w, Suv, spatialDim); + + // Flatten u -> a homogeneous curve along v (one value per column) + final Vector<?>[] alongV = new Vector<?>[cols]; + for (int j = 0; j < cols; j++) { + final Vector<?>[] colHom = new Vector<?>[rows]; + for (int i = 0; i < rows; i++) { + colHom[i] = DefaultNurbCurve.toHomogeneous(controlPoints[i][j], weights[i][j]); + } + alongV[j] = DeBoorAlgorithm.evaluate(u, colHom, knotsU, degree); + } + final Vector<?> sv = partialDerivative(v, alongV, knotsV, degree, w, Suv, spatialDim); + + return new Vector<?>[]{su, sv}; + } + + /** + * Differentiates the "flattened" homogeneous curve (already in homogeneous coordinates) and applies the quotient + * rule. + */ + private static Vector<?> partialDerivative(final double param, final Vector<?>[] flattenedHomogeneous, final double[] knots, final int degree, + final double w, final double[] surfacePoint, final int spatialDim) { + if (degree == 0) { + return Vectors.createDouble(spatialDim); + } + + final Vector<?>[] derivCtrl = DeBoorAlgorithm.derivativeControlPoints(flattenedHomogeneous, knots, degree); + final double[] derivKnots = Arrays.copyOfRange(knots, 1, knots.length - 1); + final Vector<?> deriv = DeBoorAlgorithm.evaluate(param, derivCtrl, derivKnots, degree - 1); + + final double wDeriv = deriv.get(spatialDim); + final double[] result = new double[spatialDim]; + for (int c = 0; c < spatialDim; c++) { + result[c] = (deriv.get(c) - wDeriv * surfacePoint[c]) / w; + } + return Vectors.createDouble(result.length).set(result); + } + + /** + * Unit normal at point (u,v) — defined for a 3D surface only. + */ + public Vector<?> normal(final double u, final double v) { + final Vector<?>[] d = derivative(u, v); + final Vector<?> n = d[0].cross(d[1]); + n.normalize(); + return n; + } + + public double domainStartU() { + return knotsU[degree]; + } + + public double domainEndU() { + return knotsU[controlPoints.length]; + } + + public double domainStartV() { + return knotsV[degree]; + } + + public double domainEndV() { + return knotsV[controlPoints[0].length]; + } + + // ------------------------------------------------------------------ + // Subdivision + // ------------------------------------------------------------------ + /** + * Splits this surface in two halves along the u direction, at the given u value. + * Trick: each column of the grid IS a curve in u, + * so {@link DefaultNurbCurve#subdivide(double)} is reused directly, column by column. + * The resulting left/right knot vector is the same for every column + * (it depends only on u and on the original knot vector), so it can be collected once. + * + * @param u parameter at which to split the surface + * @return the two halves, in parameter order + */ + public DefaultNurbSurface[] subdivideU(final double u) { + final int rows = controlPoints.length, cols = controlPoints[0].length; + Vector<?>[][] leftCtrl = null, rightCtrl = null; + double[][] leftW = null, rightW = null; + double[] leftKnots = null, rightKnots = null; + + for (int j = 0; j < cols; j++) { + final Vector<?>[] colCtrl = new Vector<?>[rows]; + final double[] colW = new double[rows]; + for (int i = 0; i < rows; i++) { + colCtrl[i] = controlPoints[i][j]; + colW[i] = weights[i][j]; + } + + final DataPoints colDp = new ArrayDataPoints( + NDArrays.of(List.of(colCtrl), + colCtrl[0].getSampleSystem(), + colCtrl[0].getDataType())); + final DefaultNurbCurve[] halves = new DefaultNurbCurve(colDp, colW, knotsU, degree).subdivide(u); + + if (leftCtrl == null) { + leftCtrl = new Vector<?>[halves[0].nbCtrlPts][cols]; + leftW = new double[halves[0].nbCtrlPts][cols]; + rightCtrl = new Vector<?>[halves[1].nbCtrlPts][cols]; + rightW = new double[halves[1].nbCtrlPts][cols]; + leftKnots = halves[0].knots; + rightKnots = halves[1].knots; + } + for (int i = 0; i < halves[0].nbCtrlPts; i++) { + leftCtrl[i][j] = halves[0].controlPoints[i]; + leftW[i][j] = halves[0].weights[i]; + } + for (int i = 0; i < halves[1].nbCtrlPts; i++) { + rightCtrl[i][j] = halves[1].controlPoints[i]; + rightW[i][j] = halves[1].weights[i]; + } + } + + return new DefaultNurbSurface[]{ + new DefaultNurbSurface(leftCtrl, leftW, leftKnots, knotsV, degree), + new DefaultNurbSurface(rightCtrl, rightW, rightKnots, knotsV, degree) + }; + } + + /** + * Same as {@link #subdivideU} but along the v direction (each ROW is a curve in v). + */ + public DefaultNurbSurface[] subdivideV(final double v) { + final int rows = controlPoints.length; + Vector<?>[][] leftCtrl = null, rightCtrl = null; + double[][] leftW = null, rightW = null; + double[] leftKnots = null, rightKnots = null; + + for (int i = 0; i < rows; i++) { + + final DataPoints dp = new ArrayDataPoints( + NDArrays.of(List.of(controlPoints[i]), + controlPoints[i][0].getSampleSystem(), + controlPoints[i][0].getDataType())); + + final DefaultNurbCurve[] halves = new DefaultNurbCurve(dp, weights[i], knotsV, degree).subdivide(v); + + if (leftCtrl == null) { + leftCtrl = new Vector<?>[rows][halves[0].nbCtrlPts]; + leftW = new double[rows][halves[0].controlPoints.length]; + rightCtrl = new Vector<?>[rows][halves[1].nbCtrlPts]; + rightW = new double[rows][halves[1].controlPoints.length]; + leftKnots = halves[0].knots; + rightKnots = halves[1].knots; + } + leftCtrl[i] = halves[0].controlPoints; + leftW[i] = halves[0].weights; + rightCtrl[i] = halves[1].controlPoints; + rightW[i] = halves[1].weights; + } + + return new DefaultNurbSurface[]{ + new DefaultNurbSurface(leftCtrl, leftW, knotsU, leftKnots, degree), + new DefaultNurbSurface(rightCtrl, rightW, knotsU, rightKnots, degree) + }; + } + + + /** + * Bounding box of the control points, which contains the surface (convex hull property). + * + * @return the bounding box of this surface + */ + public BBox getEnvelope() { + final BBox bbox = new BBox(controlPoints[0][0], controlPoints[0][0]); + for (final Vector<?>[] row : controlPoints) { + for (final Vector<?> p : row) { + bbox.add(p); + } + } + return bbox; + } + + @Override + public int getDegree() { + return degree; + } + + @Override + public double[] getKnots() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public KnotType getKnotSpec() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public BSplineSurfaceForm getSurfaceForm() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public boolean isPolynomial() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public double getArea() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public int getRows() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public int getColumns() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public List<DirectPosition> getControlPoints() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public DataPoints getDataPoints() { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public Curve getHorizontalCurve(double v) { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public Curve getVerticalCurve(double u) { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public DirectPosition getSurface(double u, double v) { + throw new UnsupportedOperationException("Not supported yet."); + } + + @Override + public Identifier getName() { + throw new UnsupportedOperationException("Not supported yet."); + } + +} diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java index 3aa91e0f7e..08d21abbb8 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/Intersection.java @@ -29,6 +29,8 @@ import org.apache.sis.geometries.Point; import org.apache.sis.geometries.curve.LineString; import org.apache.sis.geometries.curve.MultiLineString; import org.apache.sis.geometries.internal.shared.DefaultDataPoints; +import org.apache.sis.geometries.internal.shared.DefaultNurbCurve; +import org.apache.sis.geometries.internal.shared.DefaultNurbSurface; import org.apache.sis.geometries.mesh.MeshPrimitive; import org.apache.sis.geometries.mesh.MeshPrimitiveVisitor; import org.apache.sis.geometries.mesh.MultiMeshPrimitive; @@ -36,8 +38,14 @@ import org.apache.sis.geometries.surface.PreparedTIN; import org.apache.sis.geometries.surface.Triangle; import org.apache.sis.maths.Array; import org.apache.sis.maths.Cursor; +import static org.apache.sis.maths.Maths.clamp; +import org.apache.sis.maths.Matrices; import org.apache.sis.maths.NDArrays; +import org.apache.sis.maths.ReadOnly; import org.apache.sis.maths.Tuple; +import org.apache.sis.maths.Vector; +import org.apache.sis.maths.Vectors; +import org.opengis.geometry.Envelope; import org.opengis.referencing.operation.TransformException; @@ -50,6 +58,30 @@ public final class Intersection { private Intersection(){} + static Vector<?> subtract(final ReadOnly.Vector<?> a, final ReadOnly.Tuple<?> b) { + return a.copy().subtract(b); + } + + static Vector<?> scale(final Tuple<?> a, final double s) { + return Vectors.castOrWrap(a).copy().scale(s); + } + + /** + * Tests whether two bounding boxes overlap, within the given tolerance. + * Boxes separated by less than the tolerance are still considered as overlapping, + * which {@link Envelope#intersects Envelope.intersects(…)} cannot express. + */ + private static boolean boxesOverlap(final Envelope box1, final Envelope box2, final double tol) { + for (int i = 0, n = box1.getDimension(); i < n; i++) { + if (box1.getMaximum(i) + tol < box2.getMinimum(i) + || box2.getMaximum(i) + tol < box1.getMinimum(i)) { + return false; + } + } + return true; + } + + /** * Inherit attributes from TIN. * The POSITION attribute is ignored. @@ -124,7 +156,7 @@ public final class Intersection { final Cursor cursor = positions.cursor(); final List<Integer> values = new ArrayList<>(); while (cursor.next()) { - Tuple position = cursor.samples(); + final Tuple<?> position = cursor.samples(); if (evaluator.evaluate(position).isPresent()) { values.add(Math.toIntExact(cursor.coordinate())); } @@ -155,8 +187,8 @@ public final class Intersection { for (int i = 0, n = lines.getNumGeometries(); i < n; i++) { final LineString line = lines.getGeometryN(i); final Array segment = line.getDataPoints().getAttributeArray(AttributesType.ATT_POSITION); - final Tuple s1 = segment.get(0); - final Tuple s2 = segment.get(1); + final Tuple<?> s1 = segment.get(0); + final Tuple<?> s2 = segment.get(1); try (Stream<Triangle> stream = tin.getPatches(line.getEnvelope())) { final Iterator<Triangle> iterator = stream.iterator(); @@ -164,9 +196,9 @@ public final class Intersection { while (iterator.hasNext()) { final Triangle triangle = iterator.next(); final Array corners = triangle.getExteriorRing().getDataPoints().getAttributeArray(AttributesType.ATT_POSITION); - final Tuple c0 = corners.get(0); - final Tuple c1 = corners.get(1); - final Tuple c2 = corners.get(2); + final Tuple<?> c0 = corners.get(0); + final Tuple<?> c1 = corners.get(1); + final Tuple<?> c2 = corners.get(2); final List<Tuple> clip = SutherlandHodgman.clip(Arrays.asList(s1,s2,s1), Arrays.asList(c0,c1,c2,c0)); if (clip.size() >= 2) { //inherit attributes @@ -176,8 +208,8 @@ public final class Intersection { final double y2 = c1.get(1); final double x3 = c2.get(0); final double y3 = c2.get(1); - final Tuple p1 = clip.get(0); - final Tuple p2 = clip.get(1); + final Tuple<?> p1 = clip.get(0); + final Tuple<?> p2 = clip.get(1); final double[] bary1 = Triangle.getBarycentricValue2D(x1, y1, x2, y2, x3, y3, p1.get(0), p1.get(1), 0.0, false); final double[] bary2 = Triangle.getBarycentricValue2D(x1, y1, x2, y2, x3, y3, p2.get(0), p2.get(1), 0.0, false); final Point point1 = triangle.interpolate(bary1); @@ -195,4 +227,318 @@ public final class Intersection { return intersection; } + /** + * Intersection operations on BSpline curves and surfaces. + */ + public static interface BSpline { + /** + * Computes the intersection points between this curve and the given one (both must have the same dimension). + * Exploits the convex hull property (approximated here by a bounding box, which is simpler and good enough for + * pruning) in order to subdivide recursively down to nearly straight segments, then refines each candidate with a + * Gauss-Newton method (numerical derivatives). + * + * @param first the first curve + * @param second the curve to intersect with + * @param tol tolerance on both the size of the subdivision intervals and the final accepted distance between the + * two curves + * @return the intersection points, possibly empty + */ + public static List<CurveIntersectionPoint> intersect(final DefaultNurbCurve first, final DefaultNurbCurve second, final double tol) { + final List<CurveIntersectionPoint> results = new ArrayList<>(); + intersectRecursive(first, second, tol, results); + return results; + } + + private static void intersectRecursive(final DefaultNurbCurve c1, final DefaultNurbCurve c2, final double tol, final List<CurveIntersectionPoint> results) { + if (!boxesOverlap(c1.getEnvelope(), c2.getEnvelope(), tol)) { + return; + } + + final double span1 = c1.domainEnd() - c1.domainStart(); + final double span2 = c2.domainEnd() - c2.domainStart(); + + if (span1 <= tol && span2 <= tol) { + refineAndCollect(c1, c2, tol, results); + return; + } + + if (span1 >= span2) { + final double mid = (c1.domainStart() + c1.domainEnd()) / 2; + final DefaultNurbCurve[] halves = c1.subdivide(mid); + intersectRecursive(halves[0], c2, tol, results); + intersectRecursive(halves[1], c2, tol, results); + } else { + final double mid = (c2.domainStart() + c2.domainEnd()) / 2; + final DefaultNurbCurve[] halves = c2.subdivide(mid); + intersectRecursive(c1, halves[0], tol, results); + intersectRecursive(c1, halves[1], tol, results); + } + } + + /** + * Initial guess (nearly linear segments) then Gauss-Newton refinement. + */ + private static void refineAndCollect(final DefaultNurbCurve c1, final DefaultNurbCurve c2, final double tol, final List<CurveIntersectionPoint> results) { + final double u1a = c1.domainStart(), u1b = c1.domainEnd(); + final double u2a = c2.domainStart(), u2b = c2.domainEnd(); + + // Closest points between the two chords, as a starting guess. + final double[] start1 = c1.evaluate(u1a).toArrayDouble(); + final double[] end1 = c1.evaluate(u1b).toArrayDouble(); + final double[] start2 = c2.evaluate(u2a).toArrayDouble(); + final double[] end2 = c2.evaluate(u2b).toArrayDouble(); + final double[] ratio = new double[2]; + Distance.distanceSquare(start1, end1, new double[start1.length], + start2, end2, new double[start2.length], ratio, 1e-14); + double u1 = u1a + ratio[0] * (u1b - u1a); + double u2 = u2a + ratio[1] * (u2b - u2a); + + for (int iter = 0; iter < 10; iter++) { + final Vector<?> c1u = c1.evaluate(u1), c2u = c2.evaluate(u2); + final Vector<?> f = subtract(c1u, c2u); + + final Vector<?> d1 = numericDerivative(c1, u1); + final Vector<?> d2 = scale(numericDerivative(c2, u2), -1); + + final double a = d1.dot(d1), b = d1.dot(d2), cc = d2.dot(d2); + final double rhs0 = -d1.dot(f), rhs1 = -d2.dot(f); + final double denom = a * cc - b * b; + if (Math.abs(denom) < 1e-14) { + break; + } + + final double du1 = (rhs0 * cc - rhs1 * b) / denom; + final double du2 = (a * rhs1 - b * rhs0) / denom; + + u1 = clamp(u1 + du1, u1a, u1b); + u2 = clamp(u2 + du2, u2a, u2b); + } + + final Vector<?> p1 = c1.evaluate(u1), p2 = c2.evaluate(u2); + final double dist = subtract(p1, p2).length(); + if (dist > Math.max(tol, 1e-6)) { + return; + } + + final double u1f = u1, u2f = u2; + final boolean duplicate = results.stream().anyMatch(r + -> Math.abs(r.u1 - u1f) < 1e-4 && Math.abs(r.u2 - u2f) < 1e-4); + if (!duplicate) { + results.add(new CurveIntersectionPoint(u1, u2, p1)); + } + } + + /** + * Derivative estimated by centered finite difference, staying inside the curve domain. + */ + private static Vector<?> numericDerivative(final DefaultNurbCurve curve, final double u) { + final double lo = curve.domainStart(), hi = curve.domainEnd(); + final double h = Math.max((hi - lo) * 1e-4, 1e-8); + final double uPlus = Math.min(u + h, hi), uMinus = Math.max(u - h, lo); + final double denom = uPlus - uMinus; + if (denom == 0) { + return scale(curve.evaluate(u), 0); + } + return scale(subtract(curve.evaluate(uPlus), curve.evaluate(uMinus)), 1.0 / denom); + } + + + // ------------------------------------------------------------------ + // Intersection (bounding box + recursive quadtree subdivision + // + regularized Gauss-Newton) + // ------------------------------------------------------------------ + /** + * Computes a sampling of the intersection curve between this surface and the given one. Same principle as for + * curves (bounding box + recursive subdivision + Gauss-Newton), with two notable differences: + * + * 1) The parameter space is now 4 dimensional (u1,v1,u2,v2), so the subdivision is a quadtree: at each step, the + * direction (u or v) of the surface having the largest interval is split, among the 4 candidates. 2) The system to + * solve is under-determined (3 position equations for 4 unknowns), so unlike the curve-curve case there is no + * isolated solution but a whole curve of solutions. The Gauss-Newton normal equations are therefore regularized + * (Tikhonov ridge) so as to still get a point that "falls back" on that curve near the starting point — the result + * is thus a set of points sampling the intersection curve, not isolated points as for two curves. + * + * @param first the surface to intersect with + * @param second the surface to intersect with this one + * @param tol tolerance on the size of the subdivision patches and on the final accepted distance between the two + * surfaces + * @return points sampling the intersection curve, possibly empty + */ + public static List<SurfaceIntersectionPoint> intersect(final DefaultNurbSurface first, final DefaultNurbSurface second, final double tol) { + final List<SurfaceIntersectionPoint> results = new ArrayList<>(); + intersectRecursive(first, second, tol, results, 0); + return results; + } + + private static void intersectRecursive(final DefaultNurbSurface s1, final DefaultNurbSurface s2, final double tol, + final List<SurfaceIntersectionPoint> results, final int depth) { + if (!boxesOverlap(s1.getEnvelope(), s2.getEnvelope(), tol)) { + return; + } + + final double spanU1 = s1.domainEndU() - s1.domainStartU(); + final double spanV1 = s1.domainEndV() - s1.domainStartV(); + final double spanU2 = s2.domainEndU() - s2.domainStartU(); + final double spanV2 = s2.domainEndV() - s2.domainStartV(); + final double maxSpan = Math.max(Math.max(spanU1, spanV1), Math.max(spanU2, spanV2)); + + if (maxSpan <= tol || depth > 40) { + refineAndCollect(s1, s2, tol, results); + return; + } + + if (maxSpan == spanU1) { + final double mid = (s1.domainStartU() + s1.domainEndU()) / 2; + final DefaultNurbSurface[] halves = s1.subdivideU(mid); + intersectRecursive(halves[0], s2, tol, results, depth + 1); + intersectRecursive(halves[1], s2, tol, results, depth + 1); + } else if (maxSpan == spanV1) { + final double mid = (s1.domainStartV() + s1.domainEndV()) / 2; + final DefaultNurbSurface[] halves = s1.subdivideV(mid); + intersectRecursive(halves[0], s2, tol, results, depth + 1); + intersectRecursive(halves[1], s2, tol, results, depth + 1); + } else if (maxSpan == spanU2) { + final double mid = (s2.domainStartU() + s2.domainEndU()) / 2; + final DefaultNurbSurface[] halves = s2.subdivideU(mid); + intersectRecursive(s1, halves[0], tol, results, depth + 1); + intersectRecursive(s1, halves[1], tol, results, depth + 1); + } else { + final double mid = (s2.domainStartV() + s2.domainEndV()) / 2; + final DefaultNurbSurface[] halves = s2.subdivideV(mid); + intersectRecursive(s1, halves[0], tol, results, depth + 1); + intersectRecursive(s1, halves[1], tol, results, depth + 1); + } + } + + /** + * Initial guess (center of the patches) then regularized Gauss-Newton refinement (4 unknowns, 3 equations). + */ + private static void refineAndCollect(final DefaultNurbSurface s1, final DefaultNurbSurface s2, final double tol, + final List<SurfaceIntersectionPoint> results) { + final double u1a = s1.domainStartU(), u1b = s1.domainEndU(); + final double v1a = s1.domainStartV(), v1b = s1.domainEndV(); + final double u2a = s2.domainStartU(), u2b = s2.domainEndU(); + final double v2a = s2.domainStartV(), v2b = s2.domainEndV(); + + double u1 = (u1a + u1b) / 2, v1 = (v1a + v1b) / 2; + double u2 = (u2a + u2b) / 2, v2 = (v2a + v2b) / 2; + + final double lambda = 1e-8; // Tikhonov regularization (under-determined system) + for (int iter = 0; iter < 15; iter++) { + final Vector<?> f = subtract(s1.evaluate(u1, v1), s2.evaluate(u2, v2)); + + final Vector<?> dU1 = numericDerivativeU(s1, u1, v1); + final Vector<?> dV1 = numericDerivativeV(s1, u1, v1); + final Vector<?> dU2 = scale(numericDerivativeU(s2, u2, v2), -1); + final Vector<?> dV2 = scale(numericDerivativeV(s2, u2, v2), -1); + final Vector<?>[] cols = {dU1, dV1, dU2, dV2}; + + final double[][] JtJ = new double[4][4]; + final double[] Jtf = new double[4]; + for (int a = 0; a < 4; a++) { + Jtf[a] = -cols[a].dot(f); + for (int b = 0; b < 4; b++) { + JtJ[a][b] = cols[a].dot(cols[b]) + (a == b ? lambda : 0); + } + } + + final double[] delta = Matrices.solve(JtJ, Jtf); + if (delta == null) { + break; + } + + u1 = clamp(u1 + delta[0], u1a, u1b); + v1 = clamp(v1 + delta[1], v1a, v1b); + u2 = clamp(u2 + delta[2], u2a, u2b); + v2 = clamp(v2 + delta[3], v2a, v2b); + } + + final Vector<?> p1 = s1.evaluate(u1, v1), p2 = s2.evaluate(u2, v2); + final double dist = subtract(p1, p2).length(); + if (dist > Math.max(tol, 1e-6)) { + return; + } + + final double mergeTol = Math.max(tol * 5, 1e-4); + final boolean duplicate = results.stream().anyMatch(r + -> subtract(r.point, p1).length() < mergeTol); + if (!duplicate) { + results.add(new SurfaceIntersectionPoint(u1, v1, u2, v2, p1)); + } + } + + /** + * Partial derivative along u, estimated by centered finite difference inside the surface domain. + */ + private static Vector<?> numericDerivativeU(final DefaultNurbSurface s, final double u, final double v) { + final double lo = s.domainStartU(), hi = s.domainEndU(); + final double h = Math.max((hi - lo) * 1e-4, 1e-8); + final double uPlus = Math.min(u + h, hi), uMinus = Math.max(u - h, lo); + final double denom = uPlus - uMinus; + if (denom == 0) { + return scale(s.evaluate(u, v), 0); + } + return scale(subtract(s.evaluate(uPlus, v), s.evaluate(uMinus, v)), 1.0 / denom); + } + + /** + * Partial derivative along v, estimated by centered finite difference inside the surface domain. + */ + private static Vector<?> numericDerivativeV(final DefaultNurbSurface s, final double u, final double v) { + final double lo = s.domainStartV(), hi = s.domainEndV(); + final double h = Math.max((hi - lo) * 1e-4, 1e-8); + final double vPlus = Math.min(v + h, hi), vMinus = Math.max(v - h, lo); + final double denom = vPlus - vMinus; + if (denom == 0) { + return scale(s.evaluate(u, v), 0); + } + return scale(subtract(s.evaluate(u, vPlus), s.evaluate(u, vMinus)), 1.0 / denom); + } + + /** + * An intersection point which has been found: its two parameters and the point in space. + */ + public static final class CurveIntersectionPoint { + + public final double u1, u2; + public final Vector<?> point; + + public CurveIntersectionPoint(final double u1, final double u2, final Vector<?> point) { + this.u1 = u1; + this.u2 = u2; + this.point = point; + } + + @Override + public String toString() { + return String.format("u1=%.6f, u2=%.6f -> %s", u1, u2, point); + } + } + + /** + * An intersection point between two surfaces: its four parameters and the point in space. Unlike the curve-curve + * intersection (which yields isolated points), the intersection of two surfaces is generally a 3D CURVE (4 + * unknowns, 3 equations => 1 remaining degree of freedom). This class therefore represents a single sample of that + * curve, not an isolated point in the strict sense. + */ + public static final class SurfaceIntersectionPoint { + + public final double u1, v1, u2, v2; + public final Vector<?> point; + + public SurfaceIntersectionPoint(final double u1, final double v1, final double u2, final double v2, final Vector<?> point) { + this.u1 = u1; + this.v1 = v1; + this.u2 = u2; + this.v2 = v2; + this.point = point; + } + + @Override + public String toString() { + return String.format("u1=%.4f, v1=%.4f, u2=%.4f, v2=%.4f -> %s", u1, v1, u2, v2, point); + } + } + } + } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/SutherlandHodgman.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/SutherlandHodgman.java index b74e8f8c1a..a9e0c1b5c7 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/SutherlandHodgman.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/operation/SutherlandHodgman.java @@ -43,11 +43,11 @@ public final class SutherlandHodgman { * @return Sequence of tuple for the result polygon */ public static List<Tuple> clip(List<Tuple> subject, List<Tuple> clip){ - final List<Tuple> outputList = new ArrayList(subject); + final List<Tuple> outputList = new ArrayList<>(subject); for (int i = 0, n = clip.size() - 1; i < n; i++){ - final Tuple clipEdgeStart = clip.get(i); - final Tuple clipEdgeEnd = clip.get(i + 1); + final Tuple<?> clipEdgeStart = clip.get(i); + final Tuple<?> clipEdgeEnd = clip.get(i + 1); final List<Tuple> inputList = new ArrayList<>(outputList); if (inputList.isEmpty()) break; @@ -55,7 +55,7 @@ public final class SutherlandHodgman { Tuple start = inputList.get(inputList.size() - 1); for (int k = 0, kn = inputList.size(); k < kn; k++){ - final Tuple end = (Tuple) inputList.get(k); + final Tuple<?> end = inputList.get(k); if (isInside(clipEdgeStart, clipEdgeEnd, end)){ if (!isInside(clipEdgeStart, clipEdgeEnd, start)){ @@ -72,11 +72,11 @@ public final class SutherlandHodgman { return outputList; } - private static boolean isInside(Tuple edgeStart, Tuple edgeEnd, Tuple point){ + private static boolean isInside(Tuple<?> edgeStart, Tuple<?> edgeEnd, Tuple<?> point){ return Maths.lineSide(edgeStart, edgeEnd, point) > 0; } - private static Tuple computeIntersection(Tuple start1, Tuple end1, Tuple start2, Tuple end2){ + private static Tuple computeIntersection(Tuple<?> start1, Tuple<?> end1, Tuple<?> start2, Tuple<?> end2){ final double[] buffer1 = new double[start2.getDimension()]; final double[] buffer2 = new double[start2.getDimension()]; final double[] ratio = new double[2]; diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/Sphere.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/Sphere.java index cc9aa65f9f..01e14f3a24 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/Sphere.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/Sphere.java @@ -20,6 +20,7 @@ import java.util.List; import org.apache.sis.geometries.AttributesType; import org.apache.sis.geometries.BBox; import org.apache.sis.geometries.Curve; +import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.Geometries; import org.apache.sis.geometries.GeometryType; import org.apache.sis.geometries.internal.shared.AbstractGeometry; @@ -173,7 +174,7 @@ public final class Sphere extends AbstractGeometry implements ParametricCurveSur } @Override - public List<DirectPosition> getDataPoints() { + public DataPoints getDataPoints() { throw new UnsupportedOperationException("Not supported yet."); } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/BSplineSurface.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/BSplineSurface.java index c612ca8d62..bef96530e0 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/BSplineSurface.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/BSplineSurface.java @@ -16,10 +16,12 @@ */ package org.apache.sis.geometries.surface; +import org.apache.sis.geometries.AttributesType; import org.apache.sis.geometries.GeometryType; import org.apache.sis.geometries.curve.KnotType; import static org.opengis.annotation.Specification.ISO_19107; import org.opengis.annotation.UML; +import org.opengis.referencing.crs.CoordinateReferenceSystem; /** @@ -46,9 +48,32 @@ public interface BSplineSurface extends ParametricCurveSurface { boolean isPolynomial(); @Override - GeometryType getHorizontalCurveType(); + default CoordinateReferenceSystem getCoordinateReferenceSystem() { + return getDataPoints().getCoordinateReferenceSystem(); + } @Override - GeometryType getVerticalCurveType(); + default void setCoordinateReferenceSystem(CoordinateReferenceSystem crs) throws IllegalArgumentException { + getDataPoints().setCoordinateReferenceSystem(crs); + } + @Override + default AttributesType getAttributesType() { + return getDataPoints().getAttributesType(); + } + + @Override + default boolean isEmpty() { + return getDataPoints().isEmpty(); + } + + @Override + default GeometryType getHorizontalCurveType() { + return GeometryType.SPLINECURVE; + } + + @Override + default GeometryType getVerticalCurveType() { + return GeometryType.SPLINECURVE; + } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/NURB.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/NurbSurface.java similarity index 74% rename from incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/NURB.java rename to incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/NurbSurface.java index 9d7ccf9bf2..ab439e94ba 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/curve/NURB.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/NurbSurface.java @@ -14,17 +14,21 @@ * See the License for the specific language governing permissions and * limitations under the License. */ -package org.apache.sis.geometries.curve; - -import static org.opengis.annotation.Specification.ISO_19107; -import org.opengis.annotation.UML; - +package org.apache.sis.geometries.surface; /** + * TODO : missing in ISO:19107 ? need to recheck this one. + * * * @author Johann Sorel (Geomatys) */ -@UML(identifier="NURB", specification=ISO_19107) // section 7.13.8 -public interface NURB extends BSplineCurve { +public interface NurbSurface extends BSplineSurface { + + public static final String TYPE = "NURBSSURFACE"; + + @Override + default String getGeometryType() { + return TYPE; + } } diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/ParametricCurveSurface.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/ParametricCurveSurface.java index d4a6b4bd6e..2afecf0ed3 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/ParametricCurveSurface.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/surface/ParametricCurveSurface.java @@ -18,6 +18,7 @@ package org.apache.sis.geometries.surface; import java.util.List; import org.apache.sis.geometries.Curve; +import org.apache.sis.geometries.DataPoints; import org.apache.sis.geometries.GeometryType; import org.apache.sis.geometries.Surface; import static org.opengis.annotation.Specification.ISO_19107; @@ -45,7 +46,7 @@ public interface ParametricCurveSurface extends Surface, ReferenceSystem { @UML(identifier="dataPoints", specification=ISO_19107) // section 8.3.2.5 @Override - List<DirectPosition> getDataPoints(); + DataPoints getDataPoints(); @UML(identifier="horizontalCurveType", specification=ISO_19107) // section 8.3.2.2, 8.3.2.7 GeometryType getHorizontalCurveType(); diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Array.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Array.java index add97bd378..9fdbfc6a90 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Array.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Array.java @@ -16,6 +16,7 @@ */ package org.apache.sis.maths; +import java.util.function.Function; import java.util.stream.Stream; import java.util.stream.StreamSupport; import org.opengis.referencing.crs.CoordinateReferenceSystem; @@ -285,6 +286,22 @@ public interface Array extends NDArray { return array; } + /** + * Extract all tuples as a Vector array. + */ + default Vector<?>[] toArray() { + return toArray(0, getLength()); + } + + /** + * Extract all tuples as a Vector array. + */ + default Vector<?>[] toArray(long offset, long nbTuple) { + return stream(false).skip(offset).limit(nbTuple) + .map((Tuple<?> t) -> Vectors.create(t.getSampleSystem(), t.getDataType()).set(t)) + .toArray(Vector<?>[]::new); + } + /** * Apply given transformation to all tuples. * diff --git a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/NDArrays.java b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/NDArrays.java index da108e0303..94fc4b370d 100644 --- a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/NDArrays.java +++ b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/NDArrays.java @@ -43,6 +43,10 @@ public final class NDArrays { return of(vectors, SampleSystem.ofSize(dimension), dataType); } + public static <T extends ReadOnly.Tuple<?>> Array of(T[] vectors, int dimension, DataType dataType) { + return of(List.of(vectors), SampleSystem.ofSize(dimension), dataType); + } + public static Array of(List<? extends ReadOnly.Tuple<?>> vectors, SampleSystem type, DataType dataType) { final int dimension = type.getSize(); final Array array;
