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: Complete the tests for MeridionalDistance inverse method.
Date Thu, 25 Apr 2019 08:06: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


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new 7e66753  Complete the tests for MeridionalDistance inverse method.
7e66753 is described below

commit 7e667532207af65ef6eebc4c286b755c72c279f4
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Thu Apr 25 09:58:36 2019 +0200

    Complete the tests for MeridionalDistance inverse method.
---
 .../operation/projection/Initializer.java          |   2 +-
 .../projection/MeridionalDistanceBased.java        |  19 +++-
 .../projection/MeridionalDistanceTest.java         | 123 ++++++++++++++-------
 3 files changed, 96 insertions(+), 48 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 e5e8b8b..cb87dd4 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
@@ -281,7 +281,7 @@ 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),
+     * <p><b>Tip:</b> for ℯ₁ = [1 - √(1 - ℯ²)] / [1 + √(1 -
ℯ²)]  (Snyder 3-24),
      * invoke {@link DoubleDouble#ratio_1m_1p()} on the returned value.</p>
      */
     final DoubleDouble axisLengthRatio() {
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 961e467..5024acc 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
@@ -212,7 +212,9 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
      * 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>
+     *   <li>If φ is 0°, then this method returns 0.</li>
+     *   <li>If φ=+π/2, then this method returns a value slightly smaller than +π/2,
depending on the eccentricity.</li>
+     *   <li>If φ=-π/2, then this method returns a value slightly greater than -π/2,
depending on the eccentricity.</li>
      * </ul>
      *
      * @param  φ     latitude for which to compute the meridional distance, in radians.
@@ -221,9 +223,8 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
      * @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φ;
-        final double sinφ2    = sinφ * sinφ;
-        return c0*φ + sinφcosφ*(c1 + sinφ2*(c2 + sinφ2*(c3 + sinφ2*c4)));      // TODO:
use Math.fma with JDK9.
+        final double sinφ2 = sinφ * sinφ;
+        return c0*φ + cosφ*sinφ*(c1 + sinφ2*(c2 + sinφ2*(c3 + sinφ2*c4)));     // TODO:
use Math.fma with JDK9.
     }
 
     /**
@@ -244,14 +245,20 @@ abstract class MeridionalDistanceBased extends NormalizedProjection
{
      *
      * @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 {
+    final double latitude(final double distance) {
         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.
+        /*
+         * We could improve accuracy by continuing from here with Newton's iterative method
+         * (see MeridionalDistanceTest.inverse(…) for implementation).  However that method
+         * requires calls to meridianArc(double, …), which is itself an approximation based
+         * on series expansion. Consequently the accuracy of iterative method can not be
+         * better than `meridianArc(…)` accuracy.
+         */
         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 a53ee70..993bfda 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
@@ -21,6 +21,7 @@ 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;
+import org.apache.sis.internal.referencing.Resources;
 import org.apache.sis.test.TestUtilities;
 import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.DependsOn;
@@ -54,27 +55,7 @@ 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 = 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;
     }
 
@@ -139,6 +120,43 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
     }
 
     /**
+     * Computes the latitude of given meridional distance using Newton's method. This method
is described in
+     * <a href="https://en.wikipedia.org/wiki/Meridian_arc#The_inverse_meridian_problem_for_the_ellipsoid">Wikipedia</a>.
+     * This method depends on the forward method and can not be more accurate than that method.
For testing purposes, we
+     * use the {@linkplain #referenceMoreAccurate must accurate method implemented in this
class} regardless performance.
+     *
+     * @param  m  the meridional distance on ellipsoidal of semi-major axis length of 1.
+     * @return latitude for the given meridional distance, in radians.
+     */
+    private static double inverse(final MeridionalDistanceBased projection, final double
m) throws TransformException {
+        final double ITERATION_TOLERANCE = 1E-13;               // Should be less than `referenceMoreAccurate`
accuracy.
+        final double e2 = projection.eccentricitySquared;
+        int i = NormalizedProjection.MAXIMUM_ITERATIONS;
+        double φ = m;                                           // Starting point.
+        do {
+            /*
+             * Meridian radius of curvature:   M(φ) = a(1 - ℯ²) / (1 - ℯ²sinφ)^(3/2).
+             * We compute the inverse of M.
+             */
+            final double sinφ = sin(φ);
+            double iM = 1 - e2*(sinφ*sinφ);
+            iM *= sqrt(iM);
+            iM /= (1 - e2);
+            /*
+             * From Wikipedia: Newton's method:  φᵢ₊₁ = φᵢ - [m(φᵢ) - m] / M(φᵢ)
+             * In principle the numerator should be the derivative of m(φᵢ), but
+             * the meridian radius of curvature can be used instead.
+             */
+            final double t = (m - referenceMoreAccurate(projection, φ)) * iM;
+            φ += t;
+            if (abs(t) <= ITERATION_TOLERANCE) {
+                return φ;
+            }
+        } while (--i >= 0);
+        throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
+    }
+
+    /**
      * Compares {@link MeridionalDistanceBased#meridianArc(double, double, double)} with
formulas taken as references.
      */
     @Test
@@ -172,48 +190,71 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
     }
 
     /**
-     * Tests the {@link MeridionalDistanceBased#inverse(double)} method on an ellipsoid.
+     * Compares {@link MeridionalDistanceBased#latitude(double)} with iterative approach.
+     * On WGS 84 ellipsoid, the meridional distance <var>M</var> is between ±π/2
× 0.9983243 approximately.
      *
-     * @throws TransformException should never happen.
+     * @throws TransformException if the iterative method does not converge.
      */
     @Test
     @DependsOnMethod("compareWithReference")
-    public void testInverse() throws TransformException {
-        isDerivativeSupported = false;                                      // For focusing
on inverse transform.
-        create(true);
-        verifyInDomain(100);
+    public void compareInverse() throws TransformException {
+        final MeridionalDistanceBased projection = create(true);
+        final Random random = TestUtilities.createRandomNumberGenerator();
+        for (int i=0; i<50; i++) {
+            final double m = random.nextDouble() * (PI * 0.998) - (PI/2) * 0.998;
+            final double reference = inverse(projection, m);
+            final double actual    = projection.latitude(m);
+            assertEquals("Implementation disagrees with reference.", reference, actual, 1E-10);
+        }
     }
 
     /**
-     * Tests the {@link MeridionalDistanceBased#dM_dφ(double)} method on a sphere.
+     * Tests all methods, including {@link MeridionalDistanceBased#dM_dφ(double)}, on a
sphere.
      *
-     * @throws TransformException should never happen.
+     * @throws TransformException if an iteration does not converge.
      */
     @Test
-    @DependsOnMethod("compareWithReference")
-    public void testDerivativeOnSphere() throws TransformException {
-        isInverseTransformSupported = false;                                // For focusing
on derivative.
-        create(false);
-        verifyInDomain(20);
+    @DependsOnMethod("compareInverse")
+    public void testOnSphere() throws TransformException {
+        verifyInDomain(false, 20);
     }
 
     /**
-     * Tests the {@link MeridionalDistanceBased#dM_dφ(double)} method on an ellipsoid.
+     * Tests all methods, including {@link MeridionalDistanceBased#dM_dφ(double)}, on an
ellipsoid.
      *
-     * @throws TransformException should never happen.
+     * @throws TransformException if an iteration does not converge.
      */
     @Test
-    @DependsOnMethod("testDerivativeOnSphere")
-    public void testDerivativeOnEllipsoid() throws TransformException {
-        isInverseTransformSupported = false;                                // For focusing
on derivative.
-        create(true);
-        verifyInDomain(100);
+    @DependsOnMethod("testOnSphere")
+    public void testOnEllipsoid() throws TransformException {
+        verifyInDomain(true, 100);
     }
 
     /**
      * Verifies transform, inverse transform and derivative in the [-90 … 90]° latitude
range.
      */
-    private void verifyInDomain(final int numCoordinates) throws TransformException {
+    private void verifyInDomain(final boolean ellipsoidal, final int numCoordinates) throws
TransformException {
+        final MeridionalDistanceBased projection = create(ellipsoidal);
+        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φ);
+            }
+            @Override public MathTransform1D inverse() {
+                return new AbstractMathTransform1D() {
+                    @Override public double transform (final double M) throws TransformException
{
+                        return projection.latitude(M);
+                    }
+                    @Override public double derivative(final double φ) throws TransformException
{
+                        throw new TransformException("Unsupported");
+                    }
+                };
+            }
+        };
         verifyInDomain(new double[] {toRadians(-89)},
                        new double[] {toRadians(+89)},
                        new int[] {numCoordinates},


Mime
View raw message