sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 02/02: Implement Mollweide projection derivative (Jacobian matrix).
Date Sat, 28 Jul 2018 13:26:26 GMT
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 188f5f5d1123fc0189fdb455074d9ed7525c0a8a
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Sat Jul 28 15:25:09 2018 +0200

    Implement Mollweide projection derivative (Jacobian matrix).
---
 .../operation/projection/Mollweide.java            | 52 +++++++++++++++++-----
 .../operation/projection/MollweideTest.java        | 38 ++++++++++++++++
 2 files changed, 78 insertions(+), 12 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
index bf2933a..5f96038 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mollweide.java
@@ -23,7 +23,10 @@ import org.opengis.referencing.operation.OperationMethod;
 import org.apache.sis.util.Workaround;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.referencing.operation.transform.ContextualParameters;
+import org.apache.sis.referencing.operation.matrix.NoninvertibleMatrixException;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
+import org.apache.sis.referencing.operation.matrix.Matrix2;
+import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.internal.referencing.Resources;
 
 import static java.lang.Math.*;
@@ -104,7 +107,7 @@ public class Mollweide extends NormalizedProjection {
         final double λ = srcPts[srcOff];
         final double φ = srcPts[srcOff + 1];
         final double sinφ = sin(φ);
-        double θ = 2 * asin(φ / (PI/2));            // Actually θ′ in Snyder formulas.
+        double θp = 2 * asin(φ / (PI/2));           // θ′ in Snyder formulas.
         /*
          * If sin(φ) is 1 or -1 we are on a pole.
          * Iteration would produce NaN values.
@@ -117,18 +120,26 @@ public class Mollweide extends NormalizedProjection {
                 if (--nbIte < 0) {
                     throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
                 }
-                Δθ = (θ + sin(θ) - πsinφ) / (1 + cos(θ));
-                θ -= Δθ;
+                Δθ = (θp + sin(θp) - πsinφ) / (1 + cos(θp));
+                θp -= Δθ;
             } while (abs(Δθ) > 2*ITERATION_TOLERANCE);          // *2 because θ′
is twice the desired angle.
         }
-        θ *= 0.5;
-        dstPts[dstOff    ] = cos(θ) * λ;
-        dstPts[dstOff + 1] = sin(θ);
-        Matrix matrix = null;
-        if (derivate) {
-            throw new ProjectionException("Derivation not supported");      // TODO
+        final double θ = θp * 0.5;
+        final double x = cos(θ) * λ;
+        final double y = sin(θ);
+        if (dstPts != null) {
+            dstPts[dstOff    ] = x;
+            dstPts[dstOff + 1] = y;
+        }
+        if (!derivate) {
+            return null;
+        }
+        try {
+            // TODO: see https://issues.apache.org/jira/browse/SIS-428
+            return Matrices.inverse(inverseDerivate(x, y, θp, sinφ));
+        } catch (NoninvertibleMatrixException e) {
+            throw new ProjectionException(e);
         }
-        return matrix;
     }
 
     /**
@@ -142,8 +153,8 @@ public class Mollweide extends NormalizedProjection {
         final double x  = srcPts[srcOff];
         final double y  = srcPts[srcOff + 1];
         final double θ  = asin(y);
-        final double tθ = 2 * θ;
-        final double φ  = asin((tθ + sin(tθ)) / PI);
+        final double θp = 2 * θ;
+        final double φ  = asin((θp + sin(θp)) / PI);
         double λ = x / cos(θ);
         if (abs(λ) > 2*SQRT_2*PI) {
             λ = Double.NaN;
@@ -151,4 +162,21 @@ public class Mollweide extends NormalizedProjection {
         dstPts[dstOff]   = λ;
         dstPts[dstOff+1] = φ;
     }
+
+    /**
+     * Computes the inverse of projection derivative.
+     *
+     * @param  θp    {@code 2 * asin(y)}
+     * @param  sinφ  {@code (θp + sin(θp)) / PI}
+     *
+     * @see <a href="https://issues.apache.org/jira/browse/SIS-428">SIS-428</a>
+     */
+    private static Matrix inverseDerivate(final double x, final double y, final double θp,
final double sinφ) {
+        final double cosφ  = sqrt(1 - (sinφ * sinφ));
+        final double ym1   = 1 - y*y;
+        final double dλ_dx = 1 / sqrt(ym1);
+        final double dλ_dy = (x*y) * dλ_dx / ym1;
+        final double dφ_dy = 2*dλ_dx*(1 + cos(θp)) / (PI*cosφ);
+        return new Matrix2(dλ_dx, dλ_dy, 0, dφ_dy);
+    }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MollweideTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MollweideTest.java
index bab466e..325554a 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MollweideTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MollweideTest.java
@@ -126,4 +126,42 @@ public final strictfp class MollweideTest extends MapProjectionTestCase
{
         in[0] = Double.NaN;
         verifyTransform(out, in);
     }
+
+    /**
+     * Tests the inverse derivatives at a few points. This method compares the derivatives
computed by
+     * the projection with an estimation of derivatives computed by the finite differences
method.
+     *
+     * @throws FactoryException if an error occurred while creating the map projection.
+     * @throws TransformException if an error occurred while projecting a point.
+     *
+     * @see <a href="https://issues.apache.org/jira/browse/SIS-428">SIS-428</a>
+     */
+    @Test
+    @DependsOnMethod("testTransform")
+    public void testInverseDerivative() throws FactoryException, TransformException {
+        createProjection(false);
+        transform = transform.inverse();
+        derivativeDeltas = new double[] {100, 100};             // Approximatively 100 metres.
+        tolerance = Formulas.ANGULAR_TOLERANCE;
+        verifyDerivative(  912759.823,  5873471.956);
+        verifyDerivative(-7622861.357, -7774469.608);
+    }
+
+    /**
+     * Tests the derivatives at a few points. This method compares the derivatives computed
by
+     * the projection with an estimation of derivatives computed by the finite differences
method.
+     *
+     * @throws FactoryException if an error occurred while creating the map projection.
+     * @throws TransformException if an error occurred while projecting a point.
+     */
+    @Test
+    @DependsOnMethod("testInverseDerivative")
+    public void testDerivative() throws FactoryException, TransformException {
+        createProjection(false);
+        final double delta = (100.0 / 60) / 1852;               // Approximatively 100 metres.
+        derivativeDeltas = new double[] {delta, delta};
+        tolerance = 1E-6;                                       // More severe than Formulas.LINEAR_TOLERANCE.
+        verifyDerivative(15,  30);
+        verifyDerivative(10, -60);
+    }
 }


Mime
View raw message