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.meridianArc(…).
Date Tue, 23 Apr 2019 17:45:59 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 8a492bc  Implement MeridionalDistance.meridianArc(…).
8a492bc is described below

commit 8a492bc47a272e7ccefbe690aa44dd9377e9149d
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Tue Apr 23 19:45:00 2019 +0200

    Implement MeridionalDistance.meridianArc(…).
---
 .../projection/MeridionalDistanceBased.java        | 102 ++++++++++++++++++++-
 .../projection/MeridionalDistanceTest.java         |   9 +-
 2 files changed, 105 insertions(+), 6 deletions(-)

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 0d2b2ed..192a097 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
@@ -62,6 +62,16 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
     private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true;
 
     /**
+     * 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.
+     * The values depends on the form of equation implemented by {@code meridianArc(…)}.
We do not use the
+     * form published commonly found in publication 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.
+     */
+    private transient double c0, c1, c2, c3, c4;
+
+    /**
      * Creates a new normalized projection from the parameters computed by the given initializer.
      */
     MeridionalDistanceBased(final Initializer initializer) {
@@ -74,6 +84,11 @@ 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;
     }
 
     /**
@@ -89,7 +104,90 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
      * This method shall be invoked after {@code MeridionalDistanceBased} construction or
deserialization.
      */
     private void computeCoefficients() {
-        final double e2 = eccentricitySquared;
-        // TODO
+        final double e2  = eccentricitySquared;
+        final double e4  = e2*e2;
+        final double e6  = e2*e4;
+        final double e8  = e4*e4;
+        final double e10 = e2*e8;
+        /*
+         * 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φ
+         *                                         + …]
+         *
+         * But we do not use this equation as-is. Instead, we start from a series providing
two more terms (e⁸ and e¹⁰)
+         * 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
+         *   Conversion in the Gauss-Krüger Projection. Bulletin of the Geospatial Information
Authority of Japan, Vol.59.
+         *
+         * 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²).
+         *
+         *    2) Application of trigonometric identities:
+         *
+         *       replace:  f(φ) = A⋅sin(2φ) + B⋅sin(4φ) + C⋅sin(6φ) + D⋅sin(8φ)
                            Snyder 3-34
+         *       by:       f(φ) = sin(2φ)⋅(A′ + cos(2φ)⋅(B′ + cos(2φ)⋅(C′
+ D′⋅cos(2φ))))                   Snyder 3-35
+         *       with:       A′ = A - C
+         *                   B′ = 2B - 4D
+         *                   C′ = 4C
+         *                   D′ = 8D
+         *
+         *    3) Replacement of  sin2φ  by  2⋅sinφ⋅cosφ  and  cos2φ  by  1 - 2sin²φ:
+         *
+         *                 f(φ) = sinφ⋅cosφ⋅(C₁ + sin²φ⋅(C₂ + sin²φ⋅(C₃
+ C₄⋅sin²φ)))
+         *       with:      C₁ = 2⋅(A′ + B′ + C′ + D′)
+         *                  C₂ = 2⋅(−2B′ − 4C′ − 6D′)
+         *                  C₃ = 2⋅(4C′ + 12D′)
+         *                  C₄ = 2⋅(-8D′)
+         *
+         * They are the coefficients computed below (the C₀ coefficient is not affected
by trigonometric identities).
+         * We add the terms from smallest values to largest values for accuracy reasons.
Most denominators are power
+         * of 2, which result in exact representations in IEEE 754 double type. Derivation
from Kawase formula can be
+         * viewed at: https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods
+         */
+        c0 =  -441./0x10000 * e10  +  -175./0x4000  * e8  +    -5./0x100   * e6  +   -3./0x40
* e4  +  -1./4 * e2  +  1;
+        c1 =  1575./0x20000 * e10  +   175./0x4000  * e8  +     5./0x100   * e6  +    3./0x40
* e4  +  -3./4 * e2;
+        c2 = -2625./0x8000  * e10  +   175./0x6000  * e8  +  5120./0x60000 * e6  +  -15./0x20
* e4;
+        c3 =   735./0x800   * e10  +  2240./0x60000 * e8  +   -35./0x60    * e6;
+        c4 = -2205./0x1000  * e10  +  -315./0x400   * e8;
+    }
+
+    /**
+     * Calculates the meridian distance (M). This is the distance along the central meridian
from the equator to {@code φ}.
+     * 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  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.
+     */
+    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.
+    }
+
+    /**
+     * Gets the derivative of this {@link #meridianArc(double, double, double)} method.
+     *
+     * @return the derivative at the specified latitude.
+     */
+    final double dM_dφ(final double sinφ2, final double cosφ2) {
+        return c0 +
+               c1 * (sinφ2 -   cosφ2) + sinφ2*(
+               c2 * (sinφ2 - 3*cosφ2) + sinφ2*(
+               c3 * (sinφ2 - 5*cosφ2) + sinφ2*
+               c4 * (    1 - 7*cosφ2)));
     }
 }
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 35ee297..2db8320 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
@@ -78,7 +78,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} after making sure that
+     * implemented by {@link MeridionalDistanceBased#meridianArc(double, double, double)}
after making sure that
      * this value is in agreement with {@link #reference(double)}.
      *
      * <p>References:</p>
@@ -108,7 +108,7 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
     }
 
     /**
-     * Compares {@link MeridionalDistanceBased} calculation with formulas taken as references.
+     * Compares {@link MeridionalDistanceBased#meridianArc(double, double, double)} with
formulas taken as references.
      */
     @Test
     public void compareWithReference() {
@@ -118,8 +118,9 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
             final double φ = random.nextDouble() * PI - PI/2;
             final double reference = reference(φ);
             final double accurate  = referenceMoreAccurate(φ);
-            assertEquals("Accurate formula disagree with reference.", reference, accurate,
2E-10);
-            // TODO: invoke MeridionalDistanceBased method.
+            final double actual    = projection.meridianArc(φ, sin(φ), cos(φ));
+            assertEquals("Accurate formula disagrees with reference.", reference, accurate,
2E-10);
+            assertEquals("Implementation disagrees with reference.",   accurate,  actual,
  1E-13);
         }
     }
 }


Mime
View raw message