sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Implement MeridionalDistance.inverse(double).
Date Wed, 24 Apr 2019 18:53:00 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


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new 8e03db4  Implement MeridionalDistance.inverse(double).
8e03db4 is described below

commit 8e03db41eb065580d6e22851610dffd3eda82e72
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Apr 24 20:50:48 2019 +0200

    Implement MeridionalDistance.inverse(double).
---
 .../operation/projection/Initializer.java          |   3 +
 .../projection/MeridionalDistanceBased.java        | 105 ++++++++++++++++-----
 .../projection/MeridionalDistanceTest.java         |  75 ++++++++++-----
 3 files changed, 137 insertions(+), 46 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
index c551e42..e5e8b8b 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
@@ -280,6 +280,9 @@ final class Initializer {
     /**
      * Returns {@code b/a} where {@code a} is the semi-major axis length and {@code b} the
semi-minor axis length.
      * We retrieve this value from the eccentricity with {@code b/a = sqrt(1-ℯ²)}.
+     *
+     * <p><b>Tip:</b> for ℯ₁ = [1 - √(1 - ℯ²)] / [1 + √(1 -
ℯ²)]  (equation Snyder 3-24),
+     * invoke {@link DoubleDouble#ratio_1m_1p()} on the returned value.</p>
      */
     final DoubleDouble axisLengthRatio() {
         final DoubleDouble b = new DoubleDouble(1d);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
index 137b031..961e467 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
@@ -18,6 +18,9 @@ package org.apache.sis.referencing.operation.projection;
 
 import java.io.IOException;
 import java.io.ObjectInputStream;
+import org.apache.sis.internal.util.DoubleDouble;
+
+import static java.lang.Math.*;
 
 
 /**
@@ -63,20 +66,33 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
 
     /**
      * Coefficients for the formula implemented by the {@link #meridianArc(double, double,
double)} method.
-     * Values are computed by {@link #computeCoefficients()} at construction time or after
deserialization.
+     * Values are computed by {@link #computeCoefficients(double)} in constructor or after
deserialization.
      * The values depends on the form of equation implemented by {@code meridianArc(…)}.
We do not use the
      * form commonly found in publications since a few algebraic operations allow to replace
the
      * sin(2φ), sin(4φ), sin(6φ) and sin(8φ) terms by only sin²(φ), which is faster
to calculate.
-     * See {@link #computeCoefficients()} comments in method body for more information.
+     * See {@link #computeCoefficients(double)} comments in method body for more information.
      */
     private transient double c0, c1, c2, c3, c4;
 
     /**
+     * Coefficients used in inverse transform. Same comment than for {@link #c0}… applies.
+     */
+    private transient double ci1, ci2, ci3, ci4;
+
+    /**
+     * Denominator of <cite>rectifying latitude</cite> equation. The rectifying
latitude is computed by
+     * µ = M/(1 – ℯ²/4 – 3ℯ⁴/64 – 5ℯ⁶/256 – …)  (Snyder 7-19 with a=1).
+     */
+    private transient double rld;
+
+    /**
      * Creates a new normalized projection from the parameters computed by the given initializer.
      */
     MeridionalDistanceBased(final Initializer initializer) {
         super(initializer);
-        computeCoefficients();
+        DoubleDouble e1 = initializer.axisLengthRatio();
+        e1.ratio_1m_1p();
+        computeCoefficients(e1.doubleValue());
     }
 
     /**
@@ -84,11 +100,16 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
      */
     MeridionalDistanceBased(final MeridionalDistanceBased other) {
         super(other);
-        c0 = other.c0;
-        c1 = other.c1;
-        c2 = other.c2;
-        c3 = other.c3;
-        c4 = other.c4;
+        c0  = other.c0;
+        c1  = other.c1;
+        c2  = other.c2;
+        c3  = other.c3;
+        c4  = other.c4;
+        ci1 = other.ci1;
+        ci2 = other.ci2;
+        ci3 = other.ci3;
+        ci4 = other.ci4;
+        rld = other.rld;
     }
 
     /**
@@ -96,14 +117,17 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
      */
     private void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException
{
         in.defaultReadObject();
-        computeCoefficients();
+        double ba = sqrt(1 - eccentricitySquared);
+        computeCoefficients((1 - ba) / (1 + ba));
     }
 
     /**
      * Computes the coefficients in the series expansions from the {@link #eccentricitySquared}
value.
      * This method shall be invoked after {@code MeridionalDistanceBased} construction or
deserialization.
+     *
+     * @param  ei  value of  [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)]   (Snyder 3-24).
      */
-    private void computeCoefficients() {
+    private void computeCoefficients(final double ei) {
         final double e2  = eccentricitySquared;
         final double e4  = e2*e2;
         final double e6  = e2*e4;
@@ -112,13 +136,13 @@ abstract class MeridionalDistanceBased extends NormalizedProjection
{
         /*
          * Snyder 3-21 and EPSG guidance notes #7 part 2 (February 2019) give us the following
equation:
          *
-         *   M = a[(1 – e²/4 – 3e⁴/64  –  5e⁶/256  – …)⋅φ
-         *          – (3e²/8 + 3e⁴/32  + 45e⁶/1024 + …)⋅sin2φ
-         *                 + (15e⁴/256 + 45e⁶/1024 + …)⋅sin4φ
-         *                            – (35e⁶/3072 + …)⋅sin6φ
+         *   M = a[(1 – ℯ²/4 – 3ℯ⁴/64  –  5ℯ⁶/256  – …)⋅φ
+         *          – (3ℯ²/8 + 3ℯ⁴/32  + 45ℯ⁶/1024 + …)⋅sin2φ
+         *                 + (15ℯ⁴/256 + 45ℯ⁶/1024 + …)⋅sin4φ
+         *                            – (35ℯ⁶/3072 + …)⋅sin6φ
          *                                         + …]
          *
-         * But we do not use this equation as-is. Instead, we start from a series providing
two more terms (e⁸ and e¹⁰)
+         * But we do not use this equation as-is. Instead, we start from a series providing
two more terms (ℯ⁸ and ℯ¹⁰)
          * from the following source (also available with less terms on https://en.wikipedia.org/wiki/Meridian_arc):
          *
          *   Kawase, Kazushige (2011). A General Formula for Calculating Meridian Arc Length
and its Application to Coordinate
@@ -127,9 +151,9 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
          * That more accurate formula is implemented in MeridionalDistanceTest for comparison
purposes.
          * Then we transform that formula as below:
          *
-         *    1) Multiply by b²/a = (1 - e²). This is done by combining some terms. For
example (1 - e²)⋅(1 + ¾e²) =
-         *       (1 + ¾e²) - (1e² + ¾e⁴) = 1 - ¼e² - ¾e⁴. Note that we get the
first two terms of EPSG formula, which
-         *       already include the multiplication by (1 - e²).
+         *    1) Multiply by b²/a = (1 - ℯ²). This is done by combining some terms. For
example (1 - ℯ²)⋅(1 + ¾ℯ²) =
+         *       (1 + ¾ℯ²) - (1ℯ² + ¾ℯ⁴) = 1 - ¼ℯ² - ¾ℯ⁴. Note that
we get the first two terms of EPSG formula, which
+         *       already include the multiplication by (1 - ℯ²).
          *
          *    2) Application of trigonometric identities:
          *
@@ -159,19 +183,42 @@ abstract class MeridionalDistanceBased extends NormalizedProjection
{
         c3 =   735./0x800   * e10  +  2240./0x60000 * e8  +   -35./0x60    * e6;
         c4 = -2205./0x1000  * e10  +  -315./0x400   * e8;
      // c6 =   693./0x20000 * e10  omitted for now (not yet used).
+        /*
+         * Coefficients for inverse transform derived from Snyder 3-26 and EPSG guidance
notes:
+         *
+         *   φ = µ + (3ℯ₁ /  2 – 27ℯ₁³/ 32 + …)⋅sin2µ
+         *        + (21ℯ₁²/ 16 – 55ℯ₁⁴/ 32 + …)⋅sin4µ
+         *                   + (151ℯ₁³/ 96 + …)⋅sin6µ
+         *                  + (1097ℯ₁⁴/512 – …)⋅sin8µ
+         *                                 + …
+         *
+         *   where  ℯ₁ = [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)]                
         Snyder 3-24
+         *          µ  = M/(1 – ℯ²/4 – 3ℯ⁴/64 – 5ℯ⁶/256 – …)      
                 Snyder 7-19 with a=1
+         *
+         * µ is the rectifying latitude.
+         * Derivation of coefficients used by this class are provided in the above-cited
spreadsheet.
+         */
+        rld = (1 - 1./4*e2 - 3./64*e4 - 5./256*e6);       // Part of Snyder 7-19 for computing
rectifying latitude.
+        final double ei2 = ei*ei;
+        final double ei3 = ei*ei2;
+        final double ei4 = ei2*ei2;
+        ci1 =   657./0x40 * ei4  +    31./4 * ei3  +   21./4 * ei2  +  3./1 * ei;
+        ci2 = -5045./0x20 * ei4  +  -151./3 * ei3  +  -21./2 * ei2;
+        ci3 =  3291./0x08 * ei4  +   151./3 * ei3;
+        ci4 = -1097./0x04 * ei4;
     }
 
     /**
-     * Calculates the meridian distance (M). This is the distance along the central meridian
from the equator to {@code φ}.
+     * Computes the meridional distance (M) from equator to a given latitude φ.
      * Special cases:
      * <ul>
      *   <li>If <var>φ</var> is 0°, then this method returns 0.</li>
      * </ul>
      *
-     * @param  φ     latitude to calculate meridian distance for, in radians.
+     * @param  φ     latitude for which to compute the meridional distance, in radians.
      * @param  sinφ  value of sin(φ).
      * @param  cosφ  value of cos(φ).
-     * @return meridian distance for the given latitude, on an ellipsoid of semi-major axis
of 1.
+     * @return meridional distance for the given latitude on an ellipsoid of semi-major axis
of 1.
      */
     final double meridianArc(final double φ, final double sinφ, final double cosφ) {
         final double sinφcosφ = cosφ * sinφ;
@@ -191,4 +238,20 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
                             +   c1
                             +   c0;
     }
+
+    /**
+     * Computes latitude φ from a meridian distance <var>M</var>.
+     *
+     * @param  distance  meridian distance for which to compute the latitude.
+     * @return the latitude of given meridian distance, in radians.
+     * @throws ProjectionException if the iteration does not converge.
+     */
+    final double inverse(final double distance) throws ProjectionException {
+        double    φ  = distance / rld;            // Rectifying latitude µ used as first
approximation.
+        double sinφ  = sin(φ);
+        double cosφ  = cos(φ);
+        double sinφ2 = sinφ * sinφ;
+        φ += sinφ*cosφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4)));             
     // Snyder 3-26.
+        return φ;
+    }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
index c71cf39..a53ee70 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
@@ -17,6 +17,7 @@
 package org.apache.sis.referencing.operation.projection;
 
 import java.util.Random;
+import org.opengis.referencing.operation.MathTransform1D;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.referencing.operation.transform.AbstractMathTransform1D;
 import org.apache.sis.referencing.operation.DefaultOperationMethod;
@@ -53,8 +54,27 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
     private MeridionalDistanceBased create(final boolean ellipsoidal) {
         final DefaultOperationMethod provider = new org.apache.sis.internal.referencing.provider.Sinusoidal();
         final Sinusoidal projection = new Sinusoidal(provider, parameters(provider, ellipsoidal));
+        derivativeDeltas = new double[] {toRadians(0.01)};
         tolerance = NormalizedProjection.ANGULAR_TOLERANCE;     // = linear tolerance on
a sphere of radius 1.
-        transform = projection;
+        transform = new AbstractMathTransform1D() {
+            @Override public double transform (final double φ) {
+                return projection.meridianArc(φ, sin(φ), cos(φ));
+            }
+            @Override public double derivative(final double φ) {
+                final double sinφ = sin(φ);
+                return projection.dM_dφ(sinφ*sinφ);
+            }
+            @Override public MathTransform1D inverse() {
+                return new AbstractMathTransform1D() {
+                    @Override public double transform (final double M) throws TransformException
{
+                        return projection.inverse(M);
+                    }
+                    @Override public double derivative(final double φ) throws TransformException
{
+                        throw new TransformException("Unsupported");
+                    }
+                };
+            }
+        };
         return projection;
     }
 
@@ -74,8 +94,7 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
      * @param  φ  latitude in radians.
      * @return meridional distance from equator to the given latitude on an ellipsoid with
semi-major axis of 1.
      */
-    private double reference(final double φ) {
-        final MeridionalDistanceBased projection = (MeridionalDistanceBased) transform;
+    private static double reference(final MeridionalDistanceBased projection, final double
φ) {
         final double e2 = projection.eccentricitySquared;
         final double e4 = e2*e2;
         final double e6 = e2*e4;
@@ -92,7 +111,7 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
     /**
      * Series expansion with more terms. We use this formulas as a reference for testing
accuracy of the formula
      * implemented by {@link MeridionalDistanceBased#meridianArc(double, double, double)}
after making sure that
-     * this value is in agreement with {@link #reference(double)}.
+     * this value is in agreement with {@link #reference(MeridionalDistanceBased, double)}.
      *
      * <p>References:</p>
      * <ul>
@@ -103,8 +122,7 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
      *       on Wikipedia.</li>
      * </ul>
      */
-    private double referenceMoreAccurate(final double φ) {
-        final MeridionalDistanceBased projection = (MeridionalDistanceBased) transform;
+    private static double referenceMoreAccurate(final MeridionalDistanceBased projection,
final double φ) {
         final double e2  = projection.eccentricitySquared;
         final double e4  = e2*e2;
         final double e6  = e2*e4;
@@ -129,8 +147,8 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
         final Random random = TestUtilities.createRandomNumberGenerator();
         for (int i=0; i<100; i++) {
             final double φ = random.nextDouble() * PI - PI/2;
-            final double reference = reference(φ);
-            final double accurate  = referenceMoreAccurate(φ);
+            final double reference = reference(projection, φ);
+            final double accurate  = referenceMoreAccurate(projection, φ);
             final double actual    = projection.meridianArc(φ, sin(φ), cos(φ));
             assertEquals("Accurate formula disagrees with reference.", reference, accurate,
2E-10);
             assertEquals("Implementation disagrees with our formula.", accurate,  actual,
  1E-13);
@@ -154,13 +172,29 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
     }
 
     /**
+     * Tests the {@link MeridionalDistanceBased#inverse(double)} method on an ellipsoid.
+     *
+     * @throws TransformException should never happen.
+     */
+    @Test
+    @DependsOnMethod("compareWithReference")
+    public void testInverse() throws TransformException {
+        isDerivativeSupported = false;                                      // For focusing
on inverse transform.
+        create(true);
+        verifyInDomain(100);
+    }
+
+    /**
      * Tests the {@link MeridionalDistanceBased#dM_dφ(double)} method on a sphere.
      *
      * @throws TransformException should never happen.
      */
     @Test
+    @DependsOnMethod("compareWithReference")
     public void testDerivativeOnSphere() throws TransformException {
-        testDerivative(false);
+        isInverseTransformSupported = false;                                // For focusing
on derivative.
+        create(false);
+        verifyInDomain(20);
     }
 
     /**
@@ -171,27 +205,18 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
     @Test
     @DependsOnMethod("testDerivativeOnSphere")
     public void testDerivativeOnEllipsoid() throws TransformException {
-        testDerivative(true);
+        isInverseTransformSupported = false;                                // For focusing
on derivative.
+        create(true);
+        verifyInDomain(100);
     }
 
     /**
-     * Implementation of tests for derivative method.
+     * Verifies transform, inverse transform and derivative in the [-90 … 90]° latitude
range.
      */
-    private void testDerivative(final boolean ellipsoidal) throws TransformException {
-        final MeridionalDistanceBased projection = create(ellipsoidal);
-        isInverseTransformSupported = false;
-        derivativeDeltas = new double[] {toRadians(0.01)};
-        transform = new AbstractMathTransform1D() {
-            @Override public double transform (final double φ) {
-                return projection.meridianArc(φ, sin(φ), cos(φ));
-            }
-            @Override public double derivative(final double φ) {
-                final double sinφ = sin(φ);
-                return projection.dM_dφ(sinφ*sinφ);
-            }
-        };
+    private void verifyInDomain(final int numCoordinates) throws TransformException {
         verifyInDomain(new double[] {toRadians(-89)},
                        new double[] {toRadians(+89)},
-                       new int[] {ellipsoidal ? 100 : 20}, null);
+                       new int[] {numCoordinates},
+                       TestUtilities.createRandomNumberGenerator());
     }
 }


Mime
View raw message