You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2022/01/31 11:02:38 UTC
[sis] 02/03: Add transform derivative (Jacobian matrix) for rotated pole. This work completes https://issues.apache.org/jira/browse/SIS-533
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
commit c7a27e45c20ad8fac4d4df9b56c588e96c5475e3
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Sun Jan 30 17:06:58 2022 +0100
Add transform derivative (Jacobian matrix) for rotated pole.
This work completes https://issues.apache.org/jira/browse/SIS-533
---
.../operation/transform/RotatedPole.java | 50 ++++++++++++++++------
.../operation/transform/RotatedPoleTest.java | 34 +++++++++------
2 files changed, 60 insertions(+), 24 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java
index 5026903..1831e2d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java
@@ -39,6 +39,7 @@ import org.apache.sis.internal.util.Numerics;
import org.apache.sis.internal.util.Constants;
import org.apache.sis.metadata.iso.citation.Citations;
import org.apache.sis.referencing.ImmutableIdentifier;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.Debug;
@@ -300,25 +301,50 @@ public class RotatedPole extends AbstractMathTransform2D implements Serializable
* transform pre-concatenated to this transform, simply by subtracting λp from the longitude value.
* This is simpler than performing the rotation in Cartesian coordinates.
*/
- double λ = srcPts[srcOff];
- double φ = srcPts[srcOff+1];
- double z = sin(φ);
- double cosφ = cos(φ);
- double y = sin(λ) * cosφ;
- double x = cos(λ) * cosφ;
+ final double λ = srcPts[srcOff];
+ final double φ = srcPts[srcOff+1];
+ final double z = sin(φ);
+ final double cosφ = cos(φ);
+ final double sinλ = sin(λ);
+ final double cosλ = cos(λ);
+ final double x = cosφ * cosλ;
+ final double y = cosφ * sinλ;
+ final double y2 = y * y;
/*
* Apply the rotation around Y axis (so the y value stay unchanged)
* and convert back to spherical coordinates.
*/
- double xr = cosφp * z - sinφp * x;
- double zr = -cosφp * x - sinφp * z;
- double R = fastHypot(xr, y); // The slower hypot(…) is not needed because values are close to 1.
- dstPts[dstOff] = atan2(y, xr);
- dstPts[dstOff+1] = atan2(zr, R);
+ final double xsinφp = x * sinφp;
+ final double xcosφp = x * cosφp;
+ final double zsinφp = z * sinφp;
+ final double zcosφp = z * cosφp;
+ final double xt = zcosφp - xsinφp;
+ final double zt = -xcosφp - zsinφp;
+ final double r2 = xt*xt + y2; // yt = y in ihis algorithm.
+ final double r = sqrt(r2); // The slower hypot(…) is not needed because values are close to 1.
+ if (dstPts != null) {
+ dstPts[dstOff] = atan2(y, xt);
+ dstPts[dstOff+1] = atan2(zt, r);
+ }
if (!derivate) {
return null;
}
- throw new TransformException(); // TODO
+ /*
+ * We used WxMaxima for a first derivation of formulas below,
+ * then simplified the formulas by hand.
+ *
+ * https://svn.apache.org/repos/asf/sis/analysis/Rotated%20pole.wxmx
+ */
+ final double dxφ = cosλ * zsinφp + cosφ * cosφp;
+ final double dyφ = cosλ * zcosφp - cosφ * sinφp;
+ final double zλ = z * sinλ;
+ final double zr = zt / r;
+ final double rc = r2 + zt*zt;
+ return new Matrix2(
+ (xt*x - y2*sinφp) / r2, // ∂x/∂λ
+ -(xt*zλ + y*dxφ) / r2, // ∂x/∂φ
+ (r*cosφp - zr*(x + sinφp*xt))*y / rc, // ∂y/∂λ
+ (r*dyφ + zr*(y*zλ - dxφ*xt)) / rc); // ∂y/∂φ
}
/**
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/RotatedPoleTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/RotatedPoleTest.java
index 8320148..dba1a4c 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/RotatedPoleTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/RotatedPoleTest.java
@@ -44,6 +44,13 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
}
/**
+ * Creates a new test case.
+ */
+ public RotatedPoleTest() {
+ tolerance = Formulas.ANGULAR_TOLERANCE;
+ }
+
+ /**
* Creates a new transform which should be the inverse of current transform according the
* parameters declared in {@link RotatedPole#context}. Those parameters may be wrong even
* if the coordinates transformed by {@code transform.inverse()} are corrects because the
@@ -82,8 +89,6 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
@Test
public void testRotateSouthPoleOnGreenwich() throws FactoryException, TransformException {
transform = RotatedPole.rotateSouthPole(factory(), -60, 0, 0);
- tolerance = Formulas.ANGULAR_TOLERANCE;
- isDerivativeSupported = false;
final double[] coordinates = { // (λ,φ) coordinates to convert.
0, -51,
20, -51,
@@ -111,8 +116,6 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
@DependsOnMethod("testRotateSouthPoleOnGreenwich")
public void testRotateSouthPoleWithAngle() throws FactoryException, TransformException {
transform = RotatedPole.rotateSouthPole(factory(), -50, 20, 10);
- tolerance = Formulas.ANGULAR_TOLERANCE;
- isDerivativeSupported = false;
final double[] coordinates = { // (λ,φ) coordinates to convert.
20, -51,
80, -44,
@@ -137,8 +140,6 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
@Test
public void testRotateSouthToOppositeHemisphere() throws FactoryException, TransformException {
transform = RotatedPole.rotateSouthPole(factory(), 50, 20, 10);
- tolerance = Formulas.ANGULAR_TOLERANCE;
- isDerivativeSupported = false;
final double[] coordinates = { // (λ,φ) coordinates to convert.
20, 51,
80, 44,
@@ -170,8 +171,6 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
@Test
public void testRotateNorthPoleOnGreenwich() throws FactoryException, TransformException {
transform = RotatedPole.rotateNorthPole(factory(), 60, 0, 0);
- tolerance = Formulas.ANGULAR_TOLERANCE;
- isDerivativeSupported = false;
final double[] coordinates = { // (λ,φ) coordinates to convert.
0, 54,
20, 62,
@@ -205,8 +204,6 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
@DependsOnMethod("testRotateNorthPoleOnGreenwich")
public void testRotateNorthPole() throws FactoryException, TransformException {
transform = RotatedPole.rotateNorthPole(factory(), 70, 40, 0);
- tolerance = Formulas.ANGULAR_TOLERANCE;
- isDerivativeSupported = false;
final double[] coordinates = { // (λ,φ) coordinates to convert.
0, 54,
20, 62,
@@ -231,8 +228,6 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
@Test
public void testRotateNorthToOppositeHemisphere() throws FactoryException, TransformException {
transform = RotatedPole.rotateNorthPole(factory(), -50, 20, 10);
- tolerance = Formulas.ANGULAR_TOLERANCE;
- isDerivativeSupported = false;
final double[] coordinates = { // (λ,φ) coordinates to convert.
20, -51,
80, -44,
@@ -247,4 +242,19 @@ public final strictfp class RotatedPoleTest extends MathTransformTestCase {
inverseNorthPoleTransform();
verifyTransform(expected, coordinates);
}
+
+ /**
+ * Tests derivative.
+ *
+ * @throws FactoryException if the transform can not be created.
+ * @throws TransformException if an error occurred while computing a derivative.
+ */
+ @Test
+ public void testDerivative() throws FactoryException, TransformException {
+ transform = RotatedPole.rotateSouthPole(factory(), -50, 0, 0);
+ derivativeDeltas = new double[] {1E-6, 1E-6};
+ verifyDerivative( 0, -51);
+ verifyDerivative( 20, -58);
+ verifyDerivative(-30, -40);
+ }
}