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 ModifiedAzimuthalEquidistant.inverseTransform(…). Rewrite some equations using trigonometric identities for reducing the amount of trigonometric operations.
Date Wed, 18 Mar 2020 14:58:21 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 3fdbc4c  Implement ModifiedAzimuthalEquidistant.inverseTransform(…). Rewrite some
equations using trigonometric identities for reducing the amount of trigonometric operations.
3fdbc4c is described below

commit 3fdbc4cc2a8f2d89ae1bbb2745f130a7c1f72c73
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Mar 18 15:57:13 2020 +0100

    Implement ModifiedAzimuthalEquidistant.inverseTransform(…).
    Rewrite some equations using trigonometric identities for reducing the amount of trigonometric
operations.
---
 .../operation/projection/AlbersEqualArea.java      |  4 +-
 .../projection/ModifiedAzimuthalEquidistant.java   | 98 +++++++++++++++-------
 .../operation/projection/ObliqueStereographic.java |  4 +-
 .../operation/projection/Polyconic.java            |  3 +-
 4 files changed, 76 insertions(+), 33 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
index 92f1538..c9aa589 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -255,8 +255,8 @@ public class AlbersEqualArea extends EqualAreaProjection {
         /*
          * End of map projection. Now compute the derivative.
          */
-        final double me = 1 - eccentricitySquared;
-        final double dρ_dφ = -0.5 * nm*dqm_dφ(sinφ, cos(φ)*me) / (me*ρ);
+        final double ome = 1 - eccentricitySquared;
+        final double dρ_dφ = -0.5 * nm*dqm_dφ(sinφ, cos(φ)*ome) / (ome*ρ);
         return new Matrix2(cosθ*ρ, dρ_dφ*sinθ,          // ∂x/∂λ, ∂x/∂φ
                           -sinθ*ρ, dρ_dφ*cosθ);         // ∂y/∂λ, ∂y/∂φ
     }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
index 51e5663..9b6bf32 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
@@ -50,28 +50,42 @@ import static org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEqui
  */
 public class ModifiedAzimuthalEquidistant extends NormalizedProjection {
     /**
+     * For compatibility with different versions during deserialization.
+     */
+    private static final long serialVersionUID = -5394894530514589924L;
+
+    /**
      * Sine and cosine of the latitude of origin φ₀.
      */
     private final double sinφ0, cosφ0;
 
     /**
-     * The radius of curvature (assuming <var>a</var>=1) at the latitude of origin
-     * multiplied by the sine of that latitude.
+     * A term involving radius of curvature ν₀, the latitude of origin φ₀ and the eccentricity.
+     * The semi-major axis length <var>a</var> is omitted since it is handled
outside this class.
      */
     private final double ℯ2_ν0_sinφ0;
 
     /**
-     * The ℯ⋅sin(φ₀)/√(1 − ℯ²) term.
+     * The ℯ⋅sin(φ₀)/√(1 − ℯ²) term, used in direct projection.
      */
     private final double G;
 
     /**
      * The ℯ⋅cos(φ₀)/√(1 − ℯ²) term. This is the <var>H</var> term
in EPSG guidance notes
-     * but without the cos(α) term.
+     * but without the cos(α) term (omitted because α depends on the point to project).
+     *
+     * <p>Note that during inverse projection, EPSG guidance notes has a <var>A</var>
as:
+     * −ℯ²⋅cos²φ₀/(1 − ℯ²)⋅cos²α. We opportunistically use Hp² for that
purpose.</p>
      */
     private final double Hp;
 
     /**
+     * The 3⋅ℯ²⋅sin(φ₀)⋅cos(φ₀)/(1 − ℯ²) term. This is the <var>B</var>
term in EPSG guidance notes
+     * for inverse projection but without the terms that depend on coordinates of transformed
point.
+     */
+    private final double Bp;
+
+    /**
      * Work around for RFE #4093999 in Sun's bug database
      * ("Relax constraint on placement of this()/super() call in constructors").
      */
@@ -106,14 +120,16 @@ public class ModifiedAzimuthalEquidistant extends NormalizedProjection
{
     @Workaround(library="JDK", version="1.8")
     private ModifiedAzimuthalEquidistant(final Initializer initializer) {
         super(initializer);
-        final double φ0 = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN));
-        cosφ0 = cos(φ0);
-        sinφ0 = sin(φ0);
-        final double ν0 = initializer.radiusOfCurvature(sinφ0);
+        final double φ0, ν0, f;
+        φ0          = toRadians(initializer.getAndStore(LATITUDE_OF_ORIGIN));
+        cosφ0       = cos(φ0);
+        sinφ0       = sin(φ0);
+        ν0          = initializer.radiusOfCurvature(sinφ0);
         ℯ2_ν0_sinφ0 = eccentricitySquared * ν0 * sinφ0;
-        final double f = eccentricity / sqrt(1 - eccentricitySquared);
-        G  = f * sinφ0;
-        Hp = f * cosφ0;
+        f           = eccentricity / sqrt(1 - eccentricitySquared);
+        G           = f * sinφ0;
+        Hp          = f * cosφ0;
+        Bp          = 3*eccentricitySquared * (sinφ0*cosφ0) / (1 - eccentricitySquared);
 
         final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
         denormalize.convertBefore(0, ν0, null);
@@ -132,24 +148,32 @@ public class ModifiedAzimuthalEquidistant extends NormalizedProjection
{
                             final double[] dstPts, final int dstOff,
                             final boolean derivate) throws ProjectionException
     {
-        final double λ    = srcPts[srcOff  ];
-        final double φ    = srcPts[srcOff+1];
-        final double cosλ = cos(λ);
-        final double sinλ = sin(λ);
-        final double cosφ = cos(φ);
-        final double sinφ = sin(φ);
-        final double rν   = sqrt(1 - eccentricitySquared*(sinφ*sinφ));
-        final double Ψ    = atan((1 - eccentricitySquared)*tan(φ) + ℯ2_ν0_sinφ0*rν/cosφ);
-        final double α    = atan2(sinλ, cosφ0*tan(Ψ) - sinφ0*cosλ);
-        final double sinα = sin(α);
-        final double cosα = cos(α);
-        final double H    = cosα * Hp;
+        final double λ     = srcPts[srcOff  ];
+        final double φ     = srcPts[srcOff+1];
+        final double cosλ  = cos(λ);
+        final double sinλ  = sin(λ);
+        final double cosφ  = cos(φ);
+        final double sinφ  = sin(φ);
+        final double rν    = sqrt(1 - eccentricitySquared*(sinφ*sinφ));
+        final double tanΨ  = ((1 - eccentricitySquared)*sinφ + ℯ2_ν0_sinφ0*rν) / cosφ;
+        final double rcosΨ = sqrt(1 + tanΨ*tanΨ);
+        final double α     = atan2(sinλ, cosφ0*tanΨ - sinφ0*cosλ);
+        final double sinα  = sin(α);
+        final double cosα  = cos(α);
+        final double H     = cosα * Hp;
+        /*
+         * Equations are:    s  =  asin(cos(φ₀)⋅sin(Ψ) − sin(φ₀)⋅cos(Ψ)) ⋅
signum(cos(α))     for small α
+         *                   s  =  asin(sin(λ)⋅cos(Ψ) / sin(α))                    
          for other α
+         *
+         * Using identity:   sin(atan(x))  =  x / √(1 + x²)
+         * Rewrite as:       sin(Ψ)  =   tan(Ψ) / √(1 + tan²Ψ)
+         */
         double c;
-        if (abs(α) < ANGULAR_TOLERANCE) {
-            c = cosφ0*sin(Ψ) - sinφ0*cos(Ψ);
-            if (abs(α) > PI/2) c = -c;
+        if (abs(sinα) < ANGULAR_TOLERANCE) {
+            c = (cosφ0*tanΨ - sinφ0) / rcosΨ;
+            if (cosα < 0) c = -c;
         } else {
-            c = sinλ*cos(Ψ) / sinα;
+            c = sinλ / (rcosΨ * sinα);
         }
         c = asin(c);                    // After this line this is the `s` value in EPSG
guidance notes.
         final double s2 = c  * c;
@@ -181,7 +205,25 @@ public class ModifiedAzimuthalEquidistant extends NormalizedProjection
{
                                     final double[] dstPts, final int dstOff)
             throws ProjectionException
     {
-        throw new ProjectionException("Inverse transform not yet implemented.");
+        final double x    = srcPts[srcOff  ];
+        final double y    = srcPts[srcOff+1];
+        final double α    = atan2(x, y);                // Actually α′ in EPSG guidance
notes.
+        final double sinα = sin(α);
+        final double cosα = cos(α);
+              double negA = Hp * cosα; negA *= negA;    // mA = −A  compared to EPSG guidance
note.
+        final double B    = Bp * (1 + negA) * cosα;
+        final double D2   = x*x + y*y;
+        final double D    = sqrt(D2);                   // D = c′/ν₀, but division by
ν₀ is already done here.
+        final double J    = D + (negA*(1 -   negA)*(D2*D )/6)
+                              - (   B*(1 - 3*negA)*(D2*D2)/24);
+        final double J2   = J*J;
+        final double K    = 1 + (negA*J2/2) - (B*(J2*J)/6);
+        final double sinJ = sin(J);
+        final double sinΨ = sinφ0*cos(J) + cosφ0*sinJ*cosα;
+        final double cosΨ = sqrt(1 - sinΨ*sinΨ);
+        dstPts[dstOff  ]  = asin(sinα*sinJ / cosΨ);
+        dstPts[dstOff+1]  = atan((1 - eccentricitySquared*sinφ0*K / sinΨ) * (sinΨ/cosΨ)
+                               / (1 - eccentricitySquared));
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
index 64adcd0..d443612 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
@@ -370,11 +370,11 @@ public class ObliqueStereographic extends NormalizedProjection {
         final double ψ = log((1 + sinχ) / ((1 - sinχ)*c)) / (2*n);
         double φ = 2*atan(exp(ψ)) - PI/2;                               // First approximation
         final double he = eccentricity/2;
-        final double me = 1 - eccentricitySquared;
+        final double ome = 1 - eccentricitySquared;
         for (int it=0; it<MAXIMUM_ITERATIONS; it++) {
             final double ℯsinφ = eccentricity * sin(φ);
             final double ψi = log(tan(φ/2 + PI/4) * pow((1 - ℯsinφ) / (1 + ℯsinφ),
he));
-            final double Δφ = (ψ - ψi) * cos(φ) * (1 - ℯsinφ*ℯsinφ) / me;
+            final double Δφ = (ψ - ψi) * cos(φ) * (1 - ℯsinφ*ℯsinφ) / ome;
             φ += Δφ;
             if (!(abs(Δφ) > ITERATION_TOLERANCE)) {     // Use '!' for accepting NaN.
                 dstPts[dstOff  ] = λ;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Polyconic.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Polyconic.java
index d9ef536..b4c4bcd 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Polyconic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Polyconic.java
@@ -227,6 +227,7 @@ public class Polyconic extends MeridianArcBased {
         final double y = srcPts[srcOff+1];
         double φ = y;                           // A = (M₀ + (N-FE)/a)      — Snyder
18-18 with M₀=0, FE=0 and a=1.
         final double B = y*y + x*x;             // B = A² + ((E-FE)²/a²)    — Snyder
18-19 with FE=0 and a=1.
+        final double ome = 1 - eccentricitySquared;
         int i = MAXIMUM_ITERATIONS;
         double dφ;
         do {
@@ -239,7 +240,7 @@ public class Polyconic extends MeridianArcBased {
             final double rν    = sqrt(1 - eccentricitySquared * sinφ2);
             final double C     = rν * sinφ/cosφ;
             final double M     = distance(φ, sinφ, cosφ);
-            final double Mp    = (1 - eccentricitySquared) + sinφ2*(ci2 + sinφ2*(ci4 +
sinφ2*ci6));     // Derived from Snyder 18-17
+            final double Mp    = ome + sinφ2*(ci2 + sinφ2*(ci4 + sinφ2*ci6));     // Derived
from Snyder 18-17
             final double M2B   = M*M + B;
             final double sin2φ = sin(2*φ);
             /*


Mime
View raw message