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
The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
new 0b2d6fba0c feat(Geometry): add Vector.slerp, support Tuple arrays in
ArrayFactory
0b2d6fba0c is described below
commit 0b2d6fba0cd2912d7c4f137da76ce8b44dd1d653
Author: jsorel <[email protected]>
AuthorDate: Tue Sep 8 17:08:31 2026 +0200
feat(Geometry): add Vector.slerp, support Tuple arrays in ArrayFactory
---
.../org.apache.sis.geometry/main/module-info.java | 1 +
.../solid/polyhedron/AbstractPolyhedron.java | 5 +-
.../geometries/solid/polyhedron/Dodecahedron.java | 2 +
.../geometries/solid/polyhedron/Hexahedron.java | 2 +
.../geometries/solid/polyhedron/Icosahedron.java | 2 +
.../geometries/solid/polyhedron/Octahedron.java | 2 +
.../solid/polyhedron/RhombicTriacontahedron.java | 2 +
.../geometries/solid/polyhedron/Tetrahedron.java | 2 +
.../solid/polyhedron/TruncatedIcosahedron.java | 2 +
.../sis/geometries/spherical/GreatCircleArc.java | 40 +++++++
.../geometries/spherical/SphericalTriangle.java | 72 ++++++++++---
.../main/org/apache/sis/maths/ArrayFactory.java | 56 ++++++++--
.../main/org/apache/sis/maths/Vector.java | 22 ++++
.../main/org/apache/sis/maths/Vectors.java | 101 ++++++++++++++++++
.../geometries/spherical/GreatCircleArcTest.java | 50 +++++++++
.../spherical/SphericalTriangleTest.java | 68 ++++++++++++
.../test/org/apache/sis/maths/VectorsTest.java | 115 +++++++++++++++++++++
17 files changed, 518 insertions(+), 26 deletions(-)
diff --git a/incubator/src/org.apache.sis.geometry/main/module-info.java
b/incubator/src/org.apache.sis.geometry/main/module-info.java
index 53e2b9255f..db680f0f4c 100644
--- a/incubator/src/org.apache.sis.geometry/main/module-info.java
+++ b/incubator/src/org.apache.sis.geometry/main/module-info.java
@@ -37,6 +37,7 @@ module org.apache.sis.geometry {
exports org.apache.sis.geometries.surface;
exports org.apache.sis.geometries.solid;
exports org.apache.sis.geometries.solid.polyhedron;
+ exports org.apache.sis.geometries.spherical;
exports org.apache.sis.maths;
exports org.apache.sis.scene;
exports org.apache.sis.scene.light;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/AbstractPolyhedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/AbstractPolyhedron.java
index 0ac8a35f27..d9d399ad88 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/AbstractPolyhedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/AbstractPolyhedron.java
@@ -28,6 +28,7 @@ import org.apache.sis.geometries.surface.Polygon;
import org.apache.sis.maths.Array;
import org.apache.sis.maths.NDArrays;
import org.apache.sis.maths.ReadOnly;
+import org.apache.sis.maths.SampleSystem;
import org.apache.sis.maths.Vector3D;
import org.opengis.geometry.Envelope;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
@@ -137,7 +138,7 @@ public abstract sealed class AbstractPolyhedron extends
AbstractGeometry impleme
* @return unit direction vector
*/
protected static ReadOnly.Vector<?> fromLatLon(double latRad, double
lonRad) {
- return new Vector3D.Double().setFromLatLon(latRad, lonRad);
+ return new
Vector3D.Double(SampleSystem.cartesian(3)).setFromLatLon(latRad, lonRad);
}
/**
@@ -150,7 +151,7 @@ public abstract sealed class AbstractPolyhedron extends
AbstractGeometry impleme
* @return direction vector
*/
protected static ReadOnly.Vector<?> fromLatLon(double latRad, double
lonRad, double radius) {
- return new Vector3D.Double().setFromLatLon(latRad,
lonRad).scale(radius);
+ return new
Vector3D.Double(SampleSystem.cartesian(3)).setFromLatLon(latRad,
lonRad).scale(radius);
}
/**
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Dodecahedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Dodecahedron.java
index a79ffe04bd..508602b5d4 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Dodecahedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Dodecahedron.java
@@ -73,6 +73,8 @@ public final class Dodecahedron extends AbstractPolyhedron{
{15, 12, 11, 19, 14, 15}
};
+ public Dodecahedron(){}
+
@Override
public int getFaceCount() {
return FACES.length;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Hexahedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Hexahedron.java
index c89f800214..050614cf16 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Hexahedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Hexahedron.java
@@ -51,6 +51,8 @@ public final class Hexahedron extends AbstractPolyhedron{
{6, 4, 5, 7, 6}
};
+ public Hexahedron() {}
+
@Override
public int getFaceCount() {
return FACES.length;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Icosahedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Icosahedron.java
index 83df0b730c..5c82036213 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Icosahedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Icosahedron.java
@@ -73,6 +73,8 @@ public final class Icosahedron extends AbstractPolyhedron{
{11, 5, 7, 11}
};
+ public Icosahedron() {}
+
@Override
public int getFaceCount() {
return FACES.length;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Octahedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Octahedron.java
index af98108e87..dfcb67227a 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Octahedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Octahedron.java
@@ -51,6 +51,8 @@ public final class Octahedron extends AbstractPolyhedron{
{3, 1, 5, 3}
};
+ public Octahedron() {}
+
@Override
public int getFaceCount() {
return FACES.length;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/RhombicTriacontahedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/RhombicTriacontahedron.java
index ea05f96ff3..62c7b7698c 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/RhombicTriacontahedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/RhombicTriacontahedron.java
@@ -110,6 +110,8 @@ public final class RhombicTriacontahedron extends
AbstractPolyhedron{
{ 9, 24, 11, 27, 9}
};
+ public RhombicTriacontahedron() {}
+
@Override
public int getFaceCount() {
return FACES.length;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Tetrahedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Tetrahedron.java
index d6b88ee4c6..57d6d8f10f 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Tetrahedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/Tetrahedron.java
@@ -45,6 +45,8 @@ public final class Tetrahedron extends AbstractPolyhedron{
{2, 1, 3, 2}
};
+ public Tetrahedron(){}
+
@Override
public int getFaceCount() {
return FACES.length;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/TruncatedIcosahedron.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/TruncatedIcosahedron.java
index b05ec4b539..31283a05bd 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/TruncatedIcosahedron.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/solid/polyhedron/TruncatedIcosahedron.java
@@ -159,6 +159,8 @@ public final class TruncatedIcosahedron extends
AbstractPolyhedron{
{47, 46, 42, 43, 54, 55, 47}
};
+ public TruncatedIcosahedron() {}
+
@Override
public int getFaceCount() {
return FACES.length;
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/GreatCircleArc.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/GreatCircleArc.java
index 0d553ed186..67075662a5 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/GreatCircleArc.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/GreatCircleArc.java
@@ -77,4 +77,44 @@ public final class GreatCircleArc {
return sphere.getRadius() * Math.acos(cosAngle);
}
+ /**
+ * Get the point at the given fraction of the arc, interpolated along the
great circle.
+ * <p>
+ * A fraction of 0 returns {@linkplain #getA() A}, a fraction of 1 returns
+ * {@linkplain #getB() B} and 0.5 returns the arc middle. Fractions outside
+ * the [0 .. 1] range are clamped.
+ *
+ * @param fraction position on the arc, in range [0 .. 1]
+ * @return point at given fraction, as a unit direction vector from the
sphere center
+ * @throws IllegalArgumentException if the two arc ends are antipodal,
since they do
+ * not define a unique great circle
+ */
+ public ReadOnly.Vector<?> pointAt(double fraction) {
+ return interpolate(vecA, vecB, fraction);
+ }
+
+ /**
+ * Spherical linear interpolation between two unit direction vectors.
+ * <p>
+ * The result is a unit vector on the great circle passing through both
vectors, at the
+ * given fraction of the angle between them. Coincident vectors
interpolate to
+ * themselves. Antipodal vectors are rejected : every great circle passes
through both
+ * of them, so there is no arc to interpolate along.
+ *
+ * @param vecA first point, as a unit direction vector from the sphere
center
+ * @param vecB second point, as a unit direction vector from the sphere
center
+ * @param fraction position on the arc, in range [0 .. 1]
+ * @return point at given fraction, as a unit direction vector from the
sphere center
+ * @throws IllegalArgumentException if the two vectors are antipodal
+ * @see org.apache.sis.maths.Vector#slerp(ReadOnly.Tuple, double)
+ */
+ public static ReadOnly.Vector<?> interpolate(ReadOnly.Vector<?> vecA,
ReadOnly.Vector<?> vecB, double fraction) {
+ /*
+ * Slerp already keeps the length of unit vectors, but normalizing
removes the
+ * rounding it leaves behind. Worth the cost because subdividing a
cell feeds the
+ * result back in, so the drift would otherwise pile up over the
levels.
+ */
+ return vecA.copy().slerp(vecB, fraction).normalize();
+ }
+
}
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/SphericalTriangle.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/SphericalTriangle.java
index 19d696db85..981c965bdf 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/SphericalTriangle.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/geometries/spherical/SphericalTriangle.java
@@ -19,6 +19,7 @@ package org.apache.sis.geometries.spherical;
import org.apache.sis.geometries.solid.Sphere;
import org.apache.sis.maths.Maths;
import org.apache.sis.maths.ReadOnly;
+import org.apache.sis.util.ArgumentChecks;
/**
@@ -163,22 +164,67 @@ public final class SphericalTriangle {
* @return regular 4 triangle subdivision
*/
public SphericalTriangle[] quadSubdivide() {
- final ReadOnly.Vector<?> vecAB = middle(vecA, vecB);
- final ReadOnly.Vector<?> vecBC = middle(vecB, vecC);
- final ReadOnly.Vector<?> vecCA = middle(vecC, vecA);
- return new SphericalTriangle[]{
- new SphericalTriangle(sphere, vecA, vecAB, vecCA),
- new SphericalTriangle(sphere, vecAB, vecB, vecBC),
- new SphericalTriangle(sphere, vecCA, vecBC, vecC),
- new SphericalTriangle(sphere, vecAB, vecBC, vecCA)
- };
+ return subdivide(2);
}
/**
- * Get middle unit vector between two vectors.
- */
- private static ReadOnly.Vector<?> middle(ReadOnly.Vector<?> p,
ReadOnly.Vector<?> q) {
- return p.copy().add(q).normalize();
+ * Subdivide triangle in {@code n*n} triangles perfectly overlapping this
triangle.
+ * <p>
+ * Each edge is divided in {@code n} equal arcs, and the triangle is cut
along the
+ * lines joining those divisions, following the same construction as
+ * {@link #quadSubdivide() } generalized to any ratio. All returned
triangles
+ * keep the CCW order of this triangle.
+ * <p>
+ * The returned array lists the {@code n*(n+1)/2} triangles pointing the
same way
+ * as this triangle first, ordered by rows starting at corner A and,
within a row,
+ * from the A→B edge toward the A→C edge. The {@code n*(n-1)/2} triangles
pointing
+ * the opposite way follow, in the same row order. For {@code n = 2} this
yields
+ * corner A, corner B, corner C then the center triangle.
+ *
+ * @param n number of divisions of each edge, must be 1 or more
+ * @return regular {@code n*n} triangle subdivision
+ * @throws IllegalArgumentException if {@code n} is less than 1
+ */
+ public SphericalTriangle[] subdivide(final int n) {
+ ArgumentChecks.ensureStrictlyPositive("n", n);
+ if (n == 1) {
+ return new SphericalTriangle[]{this};
+ }
+ /*
+ * Build the triangular lattice of vertices. Row r, for r in [0 .. n],
holds r+1
+ * vertices spread on the arc going from the A→B edge to the A→C edge
at fraction
+ * r/n. Row 0 degenerates to corner A, row n is the B→C edge.
+ */
+ final ReadOnly.Vector<?>[][] lattice = new ReadOnly.Vector<?>[n+1][];
+ lattice[0] = new ReadOnly.Vector<?>[]{vecA};
+ for (int r = 1; r <= n; r++) {
+ //the last row ends on the original corners, keep them as they are
+ final ReadOnly.Vector<?> left = (r == n) ? vecB :
GreatCircleArc.interpolate(vecA, vecB, r / (double) n);
+ final ReadOnly.Vector<?> right = (r == n) ? vecC :
GreatCircleArc.interpolate(vecA, vecC, r / (double) n);
+ final ReadOnly.Vector<?>[] row = new ReadOnly.Vector<?>[r+1];
+ row[0] = left;
+ row[r] = right;
+ for (int k = 1; k < r; k++) {
+ row[k] = GreatCircleArc.interpolate(left, right, k / (double)
r);
+ }
+ lattice[r] = row;
+ }
+
+ final SphericalTriangle[] result = new SphericalTriangle[n*n];
+ int i = 0;
+ //triangles pointing the same way as this triangle
+ for (int r = 0; r < n; r++) {
+ for (int k = 0; k <= r; k++) {
+ result[i++] = new SphericalTriangle(sphere, lattice[r][k],
lattice[r+1][k], lattice[r+1][k+1]);
+ }
+ }
+ //triangles pointing the opposite way
+ for (int r = 1; r < n; r++) {
+ for (int k = 0; k < r; k++) {
+ result[i++] = new SphericalTriangle(sphere, lattice[r][k],
lattice[r+1][k+1], lattice[r][k+1]);
+ }
+ }
+ return result;
}
}
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/ArrayFactory.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/ArrayFactory.java
index cee1b99a82..5d18fce398 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/ArrayFactory.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/ArrayFactory.java
@@ -17,6 +17,7 @@
package org.apache.sis.maths;
import java.lang.foreign.SegmentAllocator;
+import java.util.Arrays;
import java.util.Collection;
import java.util.Iterator;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
@@ -173,6 +174,15 @@ public interface ArrayFactory {
throw new IllegalArgumentException("Values iterable is
not made of Tuple or primitive array");
}
}
+ } else if (values != null && values.getClass().isArray()) {
+ final Class<?> componentType =
values.getClass().getComponentType();
+ if (ReadOnly.Tuple.class.isAssignableFrom(componentType)) {
+ final int size = java.lang.reflect.Array.getLength(values);
+ if (size != 0) {
+ final ReadOnly.Tuple tuple = (ReadOnly.Tuple)
java.lang.reflect.Array.get(values, 0);
+ return tuple.getSampleSystem();
+ }
+ }
}
return SampleSystem.ofSize(1);
}
@@ -187,8 +197,13 @@ public interface ArrayFactory {
return nd.getShape();
} else if (values != null && values.getClass().isArray()) {
final int size = java.lang.reflect.Array.getLength(values);
- if ((size % nbSample) != 0) throw new
IllegalArgumentException("Values size : " + size + "is not a multiple of sample
system size : " + nbSample);
- return new long[]{size / nbSample};
+ final Class<?> componentType =
values.getClass().getComponentType();
+ if (ReadOnly.Tuple.class.isAssignableFrom(componentType)) {
+ return new long[]{size};
+ } else {
+ if ((size % nbSample) != 0) throw new
IllegalArgumentException("Values size : " + size + " is not a multiple of
sample system size : " + nbSample);
+ return new long[]{size / nbSample};
+ }
} else if (values instanceof Collection<?> col) {
return new long[]{col.size()};
}
@@ -202,7 +217,15 @@ public interface ArrayFactory {
return nd.getDataType();
} else if (values != null && values.getClass().isArray()) {
final Class<?> componentType =
values.getClass().getComponentType();
- return DataType.forPrimitiveType(componentType, false);
+ if (ReadOnly.Tuple.class.isAssignableFrom(componentType)) {
+ final int size = java.lang.reflect.Array.getLength(values);
+ if (size != 0) {
+ final ReadOnly.Tuple tuple = (ReadOnly.Tuple)
java.lang.reflect.Array.get(values, 0);
+ return tuple.getDataType();
+ }
+ } else {
+ return DataType.forPrimitiveType(componentType, false);
+ }
} else if (values instanceof Collection<?> col) {
final Iterator<?> ite = col.iterator();
if (ite.hasNext()) {
@@ -224,16 +247,27 @@ public interface ArrayFactory {
final int nbDim = target.getSampleSystem().getSize();
if (values != null) {
+ final Class<? extends Object> valueClass = values.getClass();
if (values instanceof Array array) {
target.set(0, array, 0, array.getLength());
- } else if (values != null && values.getClass().isArray()) {
- int idx = 0;
- final Cursor cursor = target.cursor();
- while (cursor.next()) {
- final Tuple tuple = cursor.samples();
- for (int i = 0; i < nbDim; i++) {
- tuple.set(i,
java.lang.reflect.Array.getDouble(values, idx));
- idx++;
+ } else if (values != null && valueClass.isArray()) {
+ Class<?> componentType = valueClass.getComponentType();
+ if (ReadOnly.Tuple.class.isAssignableFrom(componentType)) {
+ int idx = 0;
+ final Cursor cursor = target.cursor();
+ while (cursor.next()) {
+ final Tuple tuple = cursor.samples();
+ tuple.set((ReadOnly.Tuple)
java.lang.reflect.Array.get(values, idx));
+ }
+ } else {
+ int idx = 0;
+ final Cursor cursor = target.cursor();
+ while (cursor.next()) {
+ final Tuple tuple = cursor.samples();
+ for (int i = 0; i < nbDim; i++) {
+ tuple.set(i,
java.lang.reflect.Array.getDouble(values, idx));
+ idx++;
+ }
}
}
} else if (values instanceof Collection<?> col) {
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vector.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vector.java
index ce4ded62ea..8f467e7acd 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vector.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vector.java
@@ -102,6 +102,28 @@ public interface Vector<T extends Vector<T>> extends
Tuple<T>, ReadOnly.Vector<T
return (T) this;
}
+ /**
+ * Spherical linear interpolation from this vector to the other vector.
+ * <p>
+ * Where {@link #lerp(ReadOnly.Tuple, double) } moves along the straight
line joining
+ * the two vectors, this moves along the arc joining them, sweeping the
angle between
+ * them proportionally to the ratio. Both vectors having the same length,
the result
+ * keeps it too, which makes this the interpolation to use for directions.
When their
+ * lengths differ, the angle is still swept evenly and the length is
interpolated
+ * linearly.
+ *
+ * @param other vector to interpolate toward
+ * @param ratio interpolation factor, zero for this vector and one for the
other vector
+ * @return this vector
+ * @throws IllegalArgumentException if both vectors point in opposite
directions, in
+ * which case they define no unique arc to interpolate along
+ * @see Vectors#slerp(double[], double[], double, double[])
+ */
+ default T slerp(ReadOnly.Tuple<?> other, double ratio) {
+ set( Vectors.slerp(toArrayDouble(), other.toArrayDouble(), ratio));
+ return (T) this;
+ }
+
@Override
public T copy();
diff --git
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vectors.java
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vectors.java
index 7edfdb8dfd..02f73f1fec 100644
---
a/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vectors.java
+++
b/incubator/src/org.apache.sis.geometry/main/org/apache/sis/maths/Vectors.java
@@ -762,6 +762,19 @@ public final class Vectors {
return lerp(start, end, ratio, null);
}
+ /**
+ * Spherically interpolates between given vectors.
+ *
+ * @param start start vector (return value for ratio == 0.)
+ * @param end end vector (return value for ratio == 1.)
+ * @param ratio : 0 is close to start vector, 1 is on end vector
+ * @return the interpolated vector.
+ * @throws IllegalArgumentException if both vectors point in opposite
directions
+ */
+ public static double[] slerp(final double[] start, final double[] end,
final double ratio) {
+ return slerp(start, end, ratio, null);
+ }
+
//
/////////////////////////////////////////////////////////////////////////
// OPERATIONS WITH BUFFER
//////////////////////////////////////////////////
//
/////////////////////////////////////////////////////////////////////////
@@ -1192,6 +1205,94 @@ public final class Vectors {
return buffer;
}
+ /**
+ * Spherically interpolates between given vectors, result is stored in
buffer.
+ * If buffer is null, a new vector is created.
+ * Vectors must have the same size.
+ * <p>
+ * Where {@link #lerp(double[], double[], double, double[]) } moves along
the straight
+ * line joining the two vectors, this moves along the arc joining them, so
the angle
+ * swept is proportional to the ratio. Both vectors keeping the same
length, the result
+ * keeps it too, which makes this the interpolation to use for directions.
+ * <p>
+ * The direction follows the usual
+ * <var>sin((1-ratio)·Ω)/sin(Ω)·start + sin(ratio·Ω)/sin(Ω)·end</var>,
where Ω is the
+ * angle between the two vectors, and the length is interpolated linearly.
For vectors
+ * of equal length that is exactly the formula above ; for vectors of
different lengths
+ * the angle is still swept evenly, which the plain formula would not do
because it
+ * would let the longer vector pull the direction toward itself.
+ * <p>
+ * When the two vectors are nearly parallel the arc is indistinguishable
from the
+ * straight line, and dividing by <var>sin(Ω)</var> would lose precision,
so a linear
+ * interpolation is used instead. A vector of length zero has no
direction, so that
+ * case is linear too.
+ *
+ * @param start start vector (return value for ratio == 0.)
+ * @param end end vector (return value for ratio == 1.)
+ * @param ratio : 0 is close to start vector, 1 is on end vector
+ * @param buffer must have same size as start and end vectors or be null.
+ * @return the interpolated vector, buffer if not null
+ * @throws IllegalArgumentException if both vectors point in opposite
directions, in
+ * which case they define no unique arc to interpolate along
+ */
+ public static double[] slerp(final double[] start, final double[] end,
final double ratio, double[] buffer) {
+ if( start.length != end.length ) {
+ throw new IllegalArgumentException("Both vectors must have same
length.");
+ }
+ if( buffer == null ){
+ buffer = new double[start.length];
+ } else if( start.length != buffer.length ) {
+ throw new IllegalArgumentException("Buffer must have same
length as start and end vectors.");
+ }
+
+ if (ratio <= 0) {
+ System.arraycopy(start, 0, buffer, 0, start.length);
+ return buffer;
+ } else if (ratio >= 1) {
+ System.arraycopy(end, 0, buffer, 0, end.length);
+ return buffer;
+ }
+
+ /*
+ * Angle between the two directions. Dividing by the lengths keeps the
angle right
+ * for vectors which are not unit vectors. A vector of length zero has
no direction,
+ * which the parallel case below handles as a linear interpolation.
+ */
+ final double lengthStart = length(start);
+ final double lengthEnd = length(end);
+ final double lengths = lengthStart * lengthEnd;
+ final double cosAngle = (lengths > 0) ? Maths.clamp(dot(start, end) /
lengths, -1, 1) : 1;
+ final double angle = Math.acos(cosAngle);
+ final double sinAngle = Math.sin(angle);
+
+ if (sinAngle < 1E-9) {
+ if (cosAngle < 0) {
+ throw new IllegalArgumentException("Cannot interpolate between
vectors pointing in "
+ + "opposite directions, they define no unique arc.");
+ }
+ //nearly parallel vectors, the arc is the straight line joining
them
+ for(int i=0;i<start.length;i++){
+ buffer[i] = (1-ratio)*start[i] + ratio * end[i];
+ }
+ return buffer;
+ }
+
+ /*
+ * The ratios below sweep the angle evenly for vectors of length one,
so they are
+ * applied to the directions rather than to the vectors, and the
result is brought
+ * back to the interpolated length. Doing it on the vectors directly
would let the
+ * longer one pull the direction toward itself.
+ */
+ final double ratioStart = Math.sin(angle * (1-ratio)) / sinAngle /
lengthStart;
+ final double ratioEnd = Math.sin(angle * ratio ) / sinAngle /
lengthEnd;
+ final double interpolatedLength = (1-ratio)*lengthStart + ratio *
lengthEnd;
+
+ for(int i=0;i<start.length;i++){
+ buffer[i] = interpolatedLength * (ratioStart*start[i] + ratioEnd *
end[i]);
+ }
+ return buffer;
+ }
+
//
/////////////////////////////////////////////////////////////////////////
// OPERATIONS WITH MULTIPLE ELEMENTS AT THE SAME TIME
//////////////////////
//
/////////////////////////////////////////////////////////////////////////
diff --git
a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/GreatCircleArcTest.java
b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/GreatCircleArcTest.java
index 74457759a5..d742aa3356 100644
---
a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/GreatCircleArcTest.java
+++
b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/GreatCircleArcTest.java
@@ -59,4 +59,54 @@ public class GreatCircleArcTest {
assertEquals(0, emptyArc.getLength(), TOLERANCE);
}
+ /**
+ * Interpolation must land on the arc ends, on the arc middle, and space
the points
+ * evenly in angle along the way.
+ */
+ @Test
+ public void pointAtTest() {
+ final GreatCircleArc arc = new GreatCircleArc(new Sphere(3), a, b);
+
+ assertArrayEquals(a.toArrayDouble(), arc.pointAt(0).toArrayDouble(),
TOLERANCE);
+ assertArrayEquals(b.toArrayDouble(), arc.pointAt(1).toArrayDouble(),
TOLERANCE);
+ //fractions outside the range are clamped
+ assertArrayEquals(a.toArrayDouble(),
arc.pointAt(-0.5).toArrayDouble(), TOLERANCE);
+ assertArrayEquals(b.toArrayDouble(), arc.pointAt(1.5).toArrayDouble(),
TOLERANCE);
+
+ //the middle of a quarter circle is at 45 degrees
+ final double h = Math.sqrt(0.5);
+ assertArrayEquals(new double[] {h, h, 0},
arc.pointAt(0.5).toArrayDouble(), TOLERANCE);
+
+ //a third of the way is at 30 degrees, which a straight line
interpolation would miss
+ final double third = Math.PI / 6;
+ assertArrayEquals(new double[] {Math.cos(third), Math.sin(third), 0},
+ arc.pointAt(1.0 / 3).toArrayDouble(), TOLERANCE);
+
+ //every interpolated point must be a unit vector at the expected angle
from the start
+ for (int i = 0; i <= 10; i++) {
+ final double fraction = i / 10.0;
+ final ReadOnly.Vector<?> point = arc.pointAt(fraction);
+ assertEquals(1, point.length(), TOLERANCE, "Interpolated points
must be unit vectors");
+ assertEquals(fraction * Math.PI / 2, Math.acos(point.dot(a)),
TOLERANCE,
+ "Points must be evenly spaced in angle");
+ }
+ }
+
+ /**
+ * Coincident points interpolate to themselves, and antipodal points must
be rejected
+ * rather than silently produce NaN : they do not define a unique great
circle.
+ */
+ @Test
+ public void pointAtDegenerateTest() {
+ final GreatCircleArc coincident = new GreatCircleArc(new Sphere(3), a,
a);
+ assertArrayEquals(a.toArrayDouble(),
coincident.pointAt(0.0).toArrayDouble(), TOLERANCE);
+ assertArrayEquals(a.toArrayDouble(),
coincident.pointAt(0.5).toArrayDouble(), TOLERANCE);
+ assertArrayEquals(a.toArrayDouble(),
coincident.pointAt(1.0).toArrayDouble(), TOLERANCE);
+
+ final GreatCircleArc antipodal = new GreatCircleArc(new Sphere(3), a,
new Vector3D.Double(-1, 0, 0));
+ //the ends themselves are still well defined
+ assertArrayEquals(a.toArrayDouble(),
antipodal.pointAt(0.0).toArrayDouble(), TOLERANCE);
+ assertThrows(IllegalArgumentException.class, () ->
antipodal.pointAt(0.5));
+ assertThrows(IllegalArgumentException.class, () ->
antipodal.pointAt(0.25));
+ }
}
diff --git
a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/SphericalTriangleTest.java
b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/SphericalTriangleTest.java
index ac001962a8..e80fc887b1 100644
---
a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/SphericalTriangleTest.java
+++
b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/geometries/spherical/SphericalTriangleTest.java
@@ -142,4 +142,72 @@ public class SphericalTriangleTest {
return p.copy().add(q).normalize().toArrayDouble();
}
+ /**
+ * A subdivision in one is the triangle itself.
+ */
+ @Test
+ public void subdivideByOneTest() {
+ final SphericalTriangle triangle = new SphericalTriangle(new
Sphere(3), a, b, c);
+ final SphericalTriangle[] children = triangle.subdivide(1);
+ assertEquals(1, children.length);
+ assertSame(triangle, children[0]);
+ }
+
+ /**
+ * Children of a subdivision must exactly tile their parent, so their
areas must add up
+ * to the parent area, they must all keep the CCW order and their centers
must all fall
+ * inside the parent.
+ */
+ @Test
+ public void subdivideTilesTheParentTest() {
+ final SphericalTriangle triangle = new SphericalTriangle(new
Sphere(3), a, b, c);
+ for (final int n : new int[] {2, 3, 4, 5}) {
+ final SphericalTriangle[] children = triangle.subdivide(n);
+ assertEquals(n*n, children.length, "Subdivision by " + n);
+
+ double sum = 0;
+ for (SphericalTriangle child : children) {
+ sum += child.getArea();
+ /*
+ * The spherical excess is computed from the distances between
the corners,
+ * so it stays positive whatever the order. The triple product
is what tells
+ * the order apart, and it must stay positive for a CCW
triangle.
+ */
+ assertTrue(child.getA().dot(child.getB().cross(child.getC()))
> 0,
+ "Child of a subdivision by " + n + " must keep the CCW
order");
+ assertTrue(triangle.contains(child.getCentroidVector()),
+ "Child center must fall inside the parent");
+ }
+ assertEquals(triangle.getArea(), sum, 1e-9, "Children of a
subdivision by " + n + " must tile it");
+ }
+ }
+
+ /**
+ * A subdivision must reuse the parent corners as they are, and place the
new vertices
+ * on the parent edges.
+ */
+ @Test
+ public void subdivideKeepsParentCornersTest() {
+ final SphericalTriangle triangle = new SphericalTriangle(new
Sphere(3), a, b, c);
+ final SphericalTriangle[] children = triangle.subdivide(3);
+ assertEquals(9, children.length);
+ /*
+ * The 6 triangles pointing the same way as the parent come first, by
rows starting
+ * at corner A, so corner A is in the first one, corner B in the first
of the last
+ * row and corner C in the last of that row.
+ */
+ assertSame(a, children[0].getA());
+ assertSame(b, children[3].getB());
+ assertSame(c, children[5].getC());
+ }
+
+ /**
+ * Subdividing must reject a ratio below one.
+ */
+ @Test
+ public void subdivideInvalidTest() {
+ final SphericalTriangle triangle = new SphericalTriangle(new
Sphere(3), a, b, c);
+ assertThrows(IllegalArgumentException.class, () ->
triangle.subdivide(0));
+ assertThrows(IllegalArgumentException.class, () ->
triangle.subdivide(-1));
+ }
}
diff --git
a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/VectorsTest.java
b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/VectorsTest.java
index e94d32d8b6..5c09bbc92d 100644
---
a/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/VectorsTest.java
+++
b/incubator/src/org.apache.sis.geometry/test/org/apache/sis/maths/VectorsTest.java
@@ -124,4 +124,119 @@ public class VectorsTest {
assertArrayEquals(new double[] {32767, 32767, 32767},
coord2.toArrayDouble(), DELTA);
}
+
+ private static final double SLERP_TOLERANCE = 1e-12;
+
+ /**
+ * Spherical interpolation must sweep the angle evenly, and must keep the
length of
+ * the vectors, which is what sets it apart from a linear interpolation.
+ */
+ @Test
+ public void slerpTest() {
+ final double[] start = {1, 0, 0};
+ final double[] end = {0, 1, 0};
+
+ //the ends are returned as they are, and ratios outside the range are
clamped
+ assertArrayEquals(start, Vectors.slerp(start, end, 0),
SLERP_TOLERANCE);
+ assertArrayEquals(end, Vectors.slerp(start, end, 1),
SLERP_TOLERANCE);
+ assertArrayEquals(start, Vectors.slerp(start, end, -1),
SLERP_TOLERANCE);
+ assertArrayEquals(end, Vectors.slerp(start, end, 2),
SLERP_TOLERANCE);
+
+ //the middle of a quarter circle is at 45 degrees, not at the middle
of the chord
+ final double h = Math.sqrt(0.5);
+ assertArrayEquals(new double[] {h, h, 0}, Vectors.slerp(start, end,
0.5), SLERP_TOLERANCE);
+ //where a linear interpolation would fall short, inside the circle
+ assertEquals(Math.sqrt(0.5), Vectors.length(Vectors.lerp(start, end,
0.5)), SLERP_TOLERANCE);
+
+ //a third of the way is a third of the angle, so 30 degrees
+ assertArrayEquals(new double[] {Math.cos(Math.PI/6),
Math.sin(Math.PI/6), 0},
+ Vectors.slerp(start, end, 1.0/3), SLERP_TOLERANCE);
+
+ //the length is kept, and the angle grows proportionally to the ratio
+ for (int i = 0; i <= 10; i++) {
+ final double ratio = i / 10.0;
+ final double[] point = Vectors.slerp(start, end, ratio);
+ assertEquals(1, Vectors.length(point), SLERP_TOLERANCE, "Length
must be kept");
+ assertEquals(ratio * Math.PI/2, Math.acos(Vectors.dot(point,
start)), SLERP_TOLERANCE,
+ "Angle must grow proportionally to the ratio");
+ }
+ }
+
+ /**
+ * Spherical interpolation must work in any dimension, and on vectors
which are not
+ * unit vectors.
+ */
+ @Test
+ public void slerpOtherLengthsAndDimensionsTest() {
+ //2D vectors of length 3
+ final double[] start2D = {3, 0};
+ final double[] end2D = {0, 3};
+ final double[] mid2D = Vectors.slerp(start2D, end2D, 0.5);
+ assertEquals(3, Vectors.length(mid2D), SLERP_TOLERANCE);
+ assertArrayEquals(new double[] {3*Math.sqrt(0.5), 3*Math.sqrt(0.5)},
mid2D, SLERP_TOLERANCE);
+
+ //4D vectors
+ final double[] start4D = {1, 0, 0, 0};
+ final double[] end4D = {0, 0, 0, 1};
+ final double[] mid4D = Vectors.slerp(start4D, end4D, 0.5);
+ assertEquals(1, Vectors.length(mid4D), SLERP_TOLERANCE);
+ assertArrayEquals(new double[] {Math.sqrt(0.5), 0, 0, Math.sqrt(0.5)},
mid4D, SLERP_TOLERANCE);
+
+ //vectors of different lengths still sweep the angle evenly
+ final double[] shortV = {1, 0, 0};
+ final double[] longV = {0, 4, 0};
+ final double[] mid = Vectors.slerp(shortV, longV, 0.5);
+ assertEquals(Math.PI/4, Math.acos(Vectors.dot(mid, shortV) /
Vectors.length(mid)), SLERP_TOLERANCE);
+ }
+
+ /**
+ * Interpolating between aligned vectors : parallel vectors have no arc to
sweep so a
+ * linear interpolation is used, and opposite vectors define no unique arc
at all so
+ * they must be rejected rather than silently produce NaN.
+ */
+ @Test
+ public void slerpAlignedTest() {
+ final double[] a = {1, 0, 0};
+
+ //parallel vectors of the same length
+ assertArrayEquals(a, Vectors.slerp(a, a.clone(), 0.5),
SLERP_TOLERANCE);
+ //parallel vectors of different lengths interpolate linearly
+ assertArrayEquals(new double[] {2, 0, 0}, Vectors.slerp(a, new
double[] {3, 0, 0}, 0.5), SLERP_TOLERANCE);
+ //a vector of length zero has no direction, the interpolation stays
linear
+ assertArrayEquals(new double[] {0.5, 0, 0}, Vectors.slerp(a, new
double[] {0, 0, 0}, 0.5), SLERP_TOLERANCE);
+
+ //opposite vectors
+ final double[] opposite = {-1, 0, 0};
+ assertThrows(IllegalArgumentException.class, () -> Vectors.slerp(a,
opposite, 0.5));
+ //but the ends themselves are still well defined
+ assertArrayEquals(a, Vectors.slerp(a, opposite, 0), SLERP_TOLERANCE);
+ assertArrayEquals(opposite, Vectors.slerp(a, opposite, 1),
SLERP_TOLERANCE);
+ }
+
+ /**
+ * Vectors of different sizes cannot be interpolated, and a buffer must
match.
+ */
+ @Test
+ public void slerpInvalidTest() {
+ final double[] v3 = {1, 0, 0};
+ final double[] v2 = {0, 1};
+ assertThrows(IllegalArgumentException.class, () -> Vectors.slerp(v3,
v2, 0.5));
+ assertThrows(IllegalArgumentException.class, () -> Vectors.slerp(v3,
v3.clone(), 0.5, new double[2]));
+ }
+
+ /**
+ * The interpolation must also be reachable from the vector itself,
mutating it.
+ */
+ @Test
+ public void slerpOnVectorTest() {
+ final Vector3D.Double start = new Vector3D.Double(1, 0, 0);
+ final Vector3D.Double end = new Vector3D.Double(0, 1, 0);
+ final Vector3D.Double result = start.slerp(end, 0.5);
+ //the vector is mutated in place and returned
+ assertSame(start, result);
+ final double h = Math.sqrt(0.5);
+ assertArrayEquals(new double[] {h, h, 0}, result.toArrayDouble(),
SLERP_TOLERANCE);
+ //the other vector is left alone
+ assertArrayEquals(new double[] {0, 1, 0}, end.toArrayDouble(),
SLERP_TOLERANCE);
+ }
}
