sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1706925 - /sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
Date Mon, 05 Oct 2015 21:41:13 GMT
Author: desruisseaux
Date: Mon Oct  5 21:41:13 2015
New Revision: 1706925

URL: http://svn.apache.org/viewvc?rev=1706925&view=rev
Log:
Factor out some more constants (only when rewriting the equations with trigonometric identities)
using the same technic than the one we used for Lambert Conic Conformal. Actually the saving
of
a few multiplications is probably unnoticeable, but the main intend is rather to see some
more
symmetry emerging from the formulas, which is often a good sign in map projection implementations.

Modified:
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java?rev=1706925&r1=1706924&r2=1706925&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
[UTF-8] Mon Oct  5 21:41:13 2015
@@ -185,6 +185,17 @@ public class TransverseMercator extends
         final MatrixSIS denormalize = context.getMatrix(false);
         denormalize.convertBefore(0, B, null);
         denormalize.convertBefore(1, B, M0);
+        /*
+         * When rewriting equations using trigonometric identities, some constants appear.
+         * For example sin(2θ) = 2⋅sinθ⋅cosθ, so we can factor out the 2 constant
into the
+         * corresponding 'h' field.  Note: this factorization can only be performed after
+         * the constructor finished to compute other constants.
+         */
+        if (ALLOW_TRIGONOMETRIC_IDENTITIES) {
+            h2 *=  4;   ih2 *=  4;
+            h3 *= 16;   ih3 *= 16;
+            h4 *= 64;   ih4 *= 64;
+        }
     }
 
     /**
@@ -202,6 +213,11 @@ public class TransverseMercator extends
         t.sqrt();
         t.ratio_1m_1p();
         initialize(t.doubleValue());
+        if (ALLOW_TRIGONOMETRIC_IDENTITIES) {
+            h2 *=  4;   ih2 *=  4;
+            h3 *= 16;   ih3 *= 16;
+            h4 *= 64;   ih4 *= 64;
+        }
     }
 
     /**
@@ -295,13 +311,11 @@ public class TransverseMercator extends
     {
         final double λ     = srcPts[srcOff];
         final double φ     = srcPts[srcOff + 1];
-
-        final double ℯsinφ = excentricity * sin(φ);
-        final double Q     = asinh(tan(φ)) - atanh(ℯsinφ) * excentricity;
         final double sinλ  = sin(λ);
+        final double ℯsinφ = sin(φ) * excentricity;
+        final double Q     = asinh(tan(φ)) - atanh(ℯsinφ) * excentricity;
         final double coshQ = cosh(Q);
         final double η0    = atanh(sinλ / coshQ);
-
         /*
          * Original formula: η0 = atanh(sin(λ) * cos(β)) where
          * cos(β) = cos(atan(sinh(Q)))
@@ -314,16 +328,15 @@ public class TransverseMercator extends
          */
         final double coshη0 = cosh(η0);
         final double ξ0     = asin(tanh(Q) * coshη0);
-
         /*
          * Compute sin(2⋅ξ₀), sin(4⋅ξ₀), sin(6⋅ξ₀), sin(8⋅ξ₀) and same
for cos, but using the following
          * trigonometric identities in order to reduce the number of calls to Math.sin and
cos methods.
          */
-        final double sin_2ξ0, sin_4ξ0, sin_6ξ0, sin_8ξ0,
-                     cos_2ξ0, cos_4ξ0, cos_6ξ0, cos_8ξ0;
+        final double sin_2ξ0 = sin(2*ξ0);
+        final double cos_2ξ0 = cos(2*ξ0);
+        final double sin_4ξ0, sin_6ξ0, sin_8ξ0,
+                     cos_4ξ0, cos_6ξ0, cos_8ξ0;
         if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            sin_2ξ0 = sin(2*ξ0);
-            cos_2ξ0 = cos(2*ξ0);
             sin_4ξ0 = sin(4*ξ0);
             cos_4ξ0 = cos(4*ξ0);
             sin_6ξ0 = sin(6*ξ0);
@@ -331,28 +344,25 @@ public class TransverseMercator extends
             sin_8ξ0 = sin(8*ξ0);
             cos_8ξ0 = cos(8*ξ0);
         } else {
-            sin_2ξ0 = sin(2*ξ0);                              // sin(2⋅ξ₀);
-            cos_2ξ0 = cos(2*ξ0);                              // cos(2⋅ξ₀)
             final double sin2 = sin_2ξ0 * sin_2ξ0;
             final double cos2 = cos_2ξ0 * cos_2ξ0;
-            sin_4ξ0 = 2 * sin_2ξ0 * cos_2ξ0;                  // sin(4⋅ξ₀)
-            cos_4ξ0 = cos2 - sin2;                            // cos(4⋅ξ₀)
-            sin_6ξ0 = (3 - 4*sin2) * sin_2ξ0;                 // sin(6⋅ξ₀)
-            cos_6ξ0 = (4*cos2 - 3) * cos_2ξ0;                 // cos(6⋅ξ₀)
-            sin_8ξ0 = 4*cos_4ξ0 * (sin_2ξ0 * cos_2ξ0);        // sin(8⋅ξ₀)
-            cos_8ξ0 = 1 - 2*sin_4ξ0*sin_4ξ0;                  // cos(8⋅ξ₀)
+            sin_4ξ0 = sin_2ξ0 * cos_2ξ0;                      // sin(4⋅ξ₀) ÷ 2
+            cos_4ξ0 = (cos2 - sin2) * 0.5;                    // cos(4⋅ξ₀) ÷ 2
+            sin_6ξ0 = (0.75 - sin2) * sin_2ξ0;                // sin(6⋅ξ₀) ÷ 4
+            cos_6ξ0 = (cos2 - 0.75) * cos_2ξ0;                // cos(6⋅ξ₀) ÷ 4
+            sin_8ξ0 =         cos_4ξ0 * sin_4ξ0;              // sin(8⋅ξ₀) ÷ 8
+            cos_8ξ0 = 0.125 - sin_4ξ0 * sin_4ξ0;              // cos(8⋅ξ₀) ÷ 8
         }
-
         /*
          * Compute sinh(2⋅ξ₀), sinh(4⋅ξ₀), sinh(6⋅ξ₀), sinh(8⋅ξ₀) and
same for cosh, but using the following
          * hyperbolic identities in order to reduce the number of calls to Math.sinh and
cosh methods.
          * Note that the formulas are very similar to the above ones, with only some signs
reversed.
          */
-        final double sinh_2η0, sinh_4η0, sinh_6η0, sinh_8η0,
-                     cosh_2η0, cosh_4η0, cosh_6η0, cosh_8η0;
+        final double sinh_2η0 = sinh(2*η0);
+        final double cosh_2η0 = cosh(2*η0);
+        final double sinh_4η0, sinh_6η0, sinh_8η0,
+                     cosh_4η0, cosh_6η0, cosh_8η0;
         if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            sinh_2η0 = sinh(2*η0);
-            cosh_2η0 = cosh(2*η0);
             sinh_4η0 = sinh(4*η0);
             cosh_4η0 = cosh(4*η0);
             sinh_6η0 = sinh(6*η0);
@@ -360,18 +370,15 @@ public class TransverseMercator extends
             sinh_8η0 = sinh(8*η0);
             cosh_8η0 = cosh(8*η0);
         } else {
-            sinh_2η0 = sinh(2*η0);                              // sinh(2⋅η₀);
-            cosh_2η0 = cosh(2*η0);                              // cosh(2⋅η₀)
             final double sinh2 = sinh_2η0 * sinh_2η0;
             final double cosh2 = cosh_2η0 * cosh_2η0;
-            sinh_4η0 = 2 * sinh_2η0 * cosh_2η0;                 // sinh(4⋅η₀)
-            cosh_4η0 = cosh2 + sinh2;                           // cosh(4⋅η₀)
-            sinh_6η0 = (3 + 4*sinh2) * sinh_2η0;                // sinh(6⋅η₀)
-            cosh_6η0 = (4*cosh2 - 3) * cosh_2η0;                // cosh(6⋅η₀)
-            sinh_8η0 = 4*cosh_4η0 * (sinh_2η0 * cosh_2η0);      // sinh(8⋅η₀)
-            cosh_8η0 = 1 + 2*sinh_4η0*sinh_4η0;                 // cosh(8⋅η₀)
+            sinh_4η0 = sinh_2η0 * cosh_2η0;                     // sinh(4⋅η₀) ÷
2
+            cosh_4η0 = (cosh2 + sinh2) * 0.5;                   // cosh(4⋅η₀) ÷ 2
+            sinh_6η0 = (sinh2 +  0.75) * sinh_2η0;              // sinh(6⋅η₀) ÷ 4
+            cosh_6η0 = (cosh2 -  0.75) * cosh_2η0;              // cosh(6⋅η₀) ÷ 4
+            sinh_8η0 = sinh_4η0 * cosh_4η0;                     // sinh(8⋅η₀) ÷
8
+            cosh_8η0 = sinh_4η0 * sinh_4η0 + 0.125;             // cosh(8⋅η₀) ÷
8
         }
-
         /*
          * Assuming that (λ, φ) ↦ Proj((λ, φ))
          * where Proj is defined by: Proj((λ, φ)) : (η(λ, φ), ξ(λ, φ)).
@@ -396,19 +403,18 @@ public class TransverseMercator extends
             dstPts[dstOff    ] = η;
             dstPts[dstOff + 1] = ξ;
         }
-
         if (!derivate) {
             return null;
         }
 
-        final double cosλ          = cos(λ);                                     //-- λ
-        final double cosφ          = cos(φ);                                     //-- φ
-        final double cosh2Q        = coshQ * coshQ;                              //-- Q
+        final double cosλ          = cos(λ);                                        //--
λ
+        final double cosφ          = cos(φ);                                        //--
φ
+        final double cosh2Q        = coshQ * coshQ;                                 //--
Q
         final double sinhQ         = sinh(Q);
         final double tanhQ         = tanh(Q);
-        final double cosh2Q_sin2λ  = cosh2Q - sinλ * sinλ;                       //--
Qλ
-        final double sinhη0        = sinh(η0);                                   //-- η0
-        final double sqrt1_thQchη0 = sqrt(1 - tanhQ * tanhQ * coshη0 * coshη0);  //--
Qη0
+        final double cosh2Q_sin2λ  = cosh2Q - (sinλ * sinλ);                        //--
Qλ
+        final double sinhη0        = sinh(η0);                                      //--
η0
+        final double sqrt1_thQchη0 = sqrt(1 - (tanhQ * tanhQ) * (coshη0 * coshη0)); //--
Qη0
 
         //-- dQ_dλ = 0;
         final double dQ_dφ  = 1 / cosφ - excentricitySquared * cosφ / (1 - ℯsinφ *
ℯsinφ);
@@ -418,11 +424,10 @@ public class TransverseMercator extends
 
         final double dξ0_dλ = sinhQ * sinhη0 * cosλ / (cosh2Q_sin2λ * sqrt1_thQchη0);
         final double dξ0_dφ = (dQ_dφ * coshη0 / cosh2Q + dη0_dφ * sinhη0 * tanhQ)
/ sqrt1_thQchη0;
-
         /*
          * Assuming that Jac(Proj((λ, φ))) is the Jacobian matrix of Proj((λ, φ)) function.
          *
-         * So derivative Proj((λ, φ)) is defined by:
+         * So the derivative of Proj((λ, φ)) is defined by:
          *                    ┌                              ┐
          *                    │ dη(λ, φ) / dλ, dη(λ, φ) / dφ │
          * Jac              = │                              │
@@ -471,20 +476,22 @@ public class TransverseMercator extends
                                     final double[] dstPts, final int dstOff)
             throws ProjectionException
     {
-        final double η = srcPts[srcOff    ];
+        final double η = srcPts[srcOff];
         final double ξ = srcPts[srcOff + 1];
         /*
          * Following calculation of sin_2ξ, sin_4ξ, etc. is basically a copy-and-paste
of the code in transform(…).
          * Its purpose is the same than for transform(…): reduce the amount of calls to
Math.sin(double) and other
          * methods.
          */
-        final double sin_2ξ,  sin_4ξ,  sin_6ξ,  sin_8ξ,
-                     cos_2ξ,  cos_4ξ,  cos_6ξ,  cos_8ξ,
-                     sinh_2η, sinh_4η, sinh_6η, sinh_8η,
-                     cosh_2η, cosh_4η, cosh_6η, cosh_8η;
+        final double sin_2ξ  = sin (2*ξ);
+        final double cos_2ξ  = cos (2*ξ);
+        final double sinh_2η = sinh(2*η);
+        final double cosh_2η = cosh(2*η);
+        final double sin_4ξ,  sin_6ξ,  sin_8ξ,
+                     cos_4ξ,  cos_6ξ,  cos_8ξ,
+                     sinh_4η, sinh_6η, sinh_8η,
+                     cosh_4η, cosh_6η, cosh_8η;
         if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            sin_2ξ = sin(2*ξ);
-            cos_2ξ = cos(2*ξ);
             sin_4ξ = sin(4*ξ);
             cos_4ξ = cos(4*ξ);
             sin_6ξ = sin(6*ξ);
@@ -492,8 +499,6 @@ public class TransverseMercator extends
             sin_8ξ = sin(8*ξ);
             cos_8ξ = cos(8*ξ);
 
-            sinh_2η = sinh(2*η);
-            cosh_2η = cosh(2*η);
             sinh_4η = sinh(4*η);
             cosh_4η = cosh(4*η);
             sinh_6η = sinh(6*η);
@@ -501,27 +506,23 @@ public class TransverseMercator extends
             sinh_8η = sinh(8*η);
             cosh_8η = cosh(8*η);
         } else {
-            sin_2ξ = sin(2*ξ);
-            cos_2ξ = cos(2*ξ);
             final double sin2 = sin_2ξ * sin_2ξ;
             final double cos2 = cos_2ξ * cos_2ξ;
-            sin_4ξ = 2 * sin_2ξ * cos_2ξ;
-            cos_4ξ = cos2 - sin2;
-            sin_6ξ = (3 - 4*sin2) * sin_2ξ;
-            cos_6ξ = (4*cos2 - 3) * cos_2ξ;
-            sin_8ξ = 4*cos_4ξ * (sin_2ξ * cos_2ξ);
-            cos_8ξ = 1 - 2*sin_4ξ*sin_4ξ;
+            sin_4ξ = sin_2ξ * cos_2ξ;                   // sin(4⋅ξ) ÷ 2
+            cos_4ξ = (cos2 - sin2) * 0.5;               // cos(4⋅ξ) ÷ 2
+            sin_6ξ = (0.75 - sin2) * sin_2ξ;            // sin(6⋅ξ) ÷ 4
+            cos_6ξ = (cos2 - 0.75) * cos_2ξ;            // cos(6⋅ξ) ÷ 4
+            sin_8ξ =         cos_4ξ * sin_4ξ;           // sin(8⋅ξ) ÷ 8
+            cos_8ξ = 0.125 - sin_4ξ * sin_4ξ;           // cos(8⋅ξ) ÷ 8
 
-            sinh_2η = sinh(2*η);
-            cosh_2η = cosh(2*η);
             final double sinh2 = sinh_2η * sinh_2η;
             final double cosh2 = cosh_2η * cosh_2η;
-            sinh_4η = 2 * sinh_2η * cosh_2η;
-            cosh_4η = cosh2 + sinh2;
-            sinh_6η = (3 + 4*sinh2) * sinh_2η;
-            cosh_6η = (4*cosh2 - 3) * cosh_2η;
-            sinh_8η = 4*cosh_4η * (sinh_2η * cosh_2η);
-            cosh_8η = 1 + 2*sinh_4η*sinh_4η;
+            sinh_4η = sinh_2η * cosh_2η;                // sinh(4⋅η₀) ÷ 2
+            cosh_4η = (cosh2 + sinh2) * 0.5;            // cosh(4⋅η₀) ÷ 2
+            sinh_6η = (sinh2 +  0.75) * sinh_2η;        // sinh(6⋅η₀) ÷ 4
+            cosh_6η = (cosh2 -  0.75) * cosh_2η;        // cosh(6⋅η₀) ÷ 4
+            sinh_8η = sinh_4η * cosh_4η;                // sinh(8⋅η₀) ÷ 8
+            cosh_8η = sinh_4η * sinh_4η + 0.125;        // cosh(8⋅η₀) ÷ 8
         }
         /*
          * The actual inverse transform.



Mime
View raw message