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: Apply to ConformalProjection the same Snyder 3-35 formulas (trigonometric identities) than we did for MeridianArcBased. Include the coefficients in the serialization, since recomputing values after deserialization may not produce the exact same results than the values computed at construction time. Since we apply a Snyder formula now, those fields shouls be more stable than before.
Date Sun, 28 Apr 2019 22:22:33 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 851bd10  Apply to ConformalProjection the same Snyder 3-35 formulas (trigonometric
identities) than we did for MeridianArcBased. Include the coefficients in the serialization,
since recomputing values after deserialization may not produce the exact same results than
the values computed at construction time. Since we apply a Snyder formula now, those fields
shouls be more stable than before.
851bd10 is described below

commit 851bd10fafb601a75c4f4b0c33cebf557c3bcbbb
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Apr 29 00:19:14 2019 +0200

    Apply to ConformalProjection the same Snyder 3-35 formulas (trigonometric identities)
than we did for MeridianArcBased.
    Include the coefficients in the serialization, since recomputing values after deserialization
may not produce the exact
    same results than the values computed at construction time. Since we apply a Snyder formula
now, those fields shouls be
    more stable than before.
---
 .../operation/projection/ConformalProjection.java  | 220 ++++++---------------
 .../operation/projection/EqualAreaProjection.java  |   4 -
 .../projection/LambertConicConformal.java          |   1 -
 .../referencing/operation/projection/Mercator.java |   1 -
 .../operation/projection/MeridianArcBased.java     |   4 +-
 .../operation/projection/ObliqueMercator.java      |   1 -
 .../operation/projection/PolarStereographic.java   |   1 -
 .../operation/projection/TransverseMercator.java   | 182 ++++++++++-------
 .../projection/ConformalProjectionTest.java        |   2 +-
 .../projection/MercatorMethodComparison.java       |   2 +-
 .../sis/referencing/operation/projection/NoOp.java |   1 -
 11 files changed, 175 insertions(+), 244 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
index 3c598bd..35f0151 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
@@ -16,8 +16,6 @@
  */
 package org.apache.sis.referencing.operation.projection;
 
-import java.io.IOException;
-import java.io.ObjectInputStream;
 import org.apache.sis.internal.referencing.Resources;
 
 import static java.lang.Math.*;
@@ -26,9 +24,12 @@ import static java.lang.Math.*;
 /**
  * Base class of {@link LambertConicConformal}, {@link Mercator} and {@link PolarStereographic}
projections.
  * All those projections have in common the property of being <cite>conformal</cite>,
i.e. they preserve
- * angles locally. However we do not put this base class in public API because we do not
(yet) guarantee
- * than all conformal projections will extend this base class.
- * Note that no projection can be both conformal and equal-area.
+ * angles locally. However we do not put this base class in public API because not all conformal
projections
+ * will extend this base class. For example {@link TransverseMercator}, despite being conformal,
uses very
+ * different formulas.
+ *
+ * <p>Note that no projection can be both conformal and equal-area. This restriction
is implemented in class
+ * hierarchy with {@link ConformalProjection} and {@link EqualAreaProjection} being two distinct
classes.</p>
  *
  * <p>This base class can been seen as a generalization of <cite>Lambert Conic
Conformal</cite> projection,
  * which includes some other projections like Mercator and Polar Stereographic as special
cases.
@@ -54,7 +55,7 @@ import static java.lang.Math.*;
  * parallel is a pole, the polar Stereographic results.” (Snyder, page 105)</div>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.7
+ * @version 1.0
  * @since   0.6
  * @module
  */
@@ -62,66 +63,7 @@ abstract class ConformalProjection extends NormalizedProjection {
     /**
      * For cross-version compatibility.
      */
-    private static final long serialVersionUID = 458860570536642265L;
-
-    /**
-     * {@code false} for using the original formulas as published by EPSG, or {@code true}
for using formulas
-     * modified using trigonometric identities. The use of trigonometric identities is for
reducing the amount
-     * of calls to the {@link Math#sin(double)} and similar methods. Some identities used
are:
-     *
-     * <ul>
-     *   <li>sin(2θ) = 2⋅sinθ⋅cosθ</li>
-     *   <li>cos(2θ) = cos²θ - sin²θ</li>
-     *   <li>sin(3θ) = (3 - 4⋅sin²θ)⋅sinθ</li>
-     *   <li>cos(3θ) = (4⋅cos³θ) - 3⋅cosθ</li>
-     *   <li>sin(4θ) = (4 - 8⋅sin²θ)⋅sinθ⋅cosθ</li>
-     *   <li>cos(4θ) = (8⋅cos⁴θ) - (8⋅cos²θ) + 1</li>
-     * </ul>
-     *
-     * Hyperbolic formulas (used in Transverse Mercator projection):
-     *
-     * <ul>
-     *   <li>sinh(2θ) = 2⋅sinhθ⋅coshθ</li>
-     *   <li>cosh(2θ) = cosh²θ + sinh²θ   =   2⋅cosh²θ - 1   =   1 + 2⋅sinh²θ</li>
-     *   <li>sinh(3θ) = (3 + 4⋅sinh²θ)⋅sinhθ</li>
-     *   <li>cosh(3θ) = ((4⋅cosh²θ) - 3)⋅coshθ</li>
-     *   <li>sinh(4θ) = (1 + 2⋅sinh²θ)⋅4.sinhθ⋅coshθ
-     *                = 4.cosh(2θ).sinhθ⋅coshθ</li>
-     *   <li>cosh(4θ) = (8⋅cosh⁴θ) - (8⋅cosh²θ) + 1
-     *                = 8⋅cosh²(θ) ⋅ (cosh²θ - 1) + 1
-     *                = 8⋅cosh²(θ) ⋅ sinh²(θ) + 1
-     *                = 2⋅sinh²(2θ) + 1</li>
-     * </ul>
-     *
-     * Snyder 3-34 to 3-39 uses above identities for deriving the following equivalence:
-     *
-     * <pre>
-     *     If:     f(φ) = A⋅sin(2φ) + B⋅sin(4φ) + C⋅sin(6φ) + D⋅sin(8φ)
-     *     Then:   f(φ) = sin(2φ)⋅(A′ + cos(2φ)⋅(B′ + cos(2φ)⋅(C′ + D′⋅cos(2φ))))
-     *     Where:  A′ = A - C
-     *             B′ = 2B - 4D
-     *             C′ = 4C
-     *             D′ = 8D
-     * </pre>
-     *
-     * Similar, but with cosine instead than sin and the addition of a constant:
-     *
-     * <pre>
-     *     If:     f(φ) = A + B⋅cos(2φ) + C⋅cos(4φ) + D⋅cos(6φ) + E⋅cos(8φ)
-     *     Then:   f(φ) = A′ + cos(2φ)⋅(B′ + cos(2φ)⋅(C′ + cos(2φ)⋅(D′
+ E′⋅cos(2φ))))
-     *     Where:  A′ = A - C + E
-     *             B′ = B - 3D
-     *             C′ = 2C - 8E
-     *             D′ = 4D
-     *             E′ = 8E
-     * </pre>
-     *
-     * Note that since this boolean is static final, the compiler should exclude the code
in the branch that is never
-     * executed (no need to comment-out that code).
-     *
-     * @see #identityEquals(double, double)
-     */
-    static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true;
+    private static final long serialVersionUID = -292755890184138414L;
 
     /**
      * The threshold value of {@link #eccentricity} at which we consider the accuracy of
the
@@ -139,88 +81,85 @@ abstract class ConformalProjection extends NormalizedProjection {
     /**
      * {@code true} if the {@link #eccentricity} value is greater than or equals to {@link
#ECCENTRICITY_THRESHOLD},
      * in which case the {@link #φ(double)} method will need to use an iterative method.
-     *
-     * <p><strong>Consider this field as final!</strong>
-     * It is not final only for the purpose of {@link #readObject(ObjectInputStream)}.
-     * This field is not serialized because its value may depend on the version of this
-     * {@code ConformalProjection} class.</p>
      */
-    private transient boolean useIterations;
+    private final boolean useIterations;
 
     /**
-     * Coefficients in the series expansion of the inverse projection,
-     * depending only on {@linkplain #eccentricity eccentricity} value.
-     * The series expansion is of the following form, where f(θ) is typically sin(θ):
+     * Coefficients in the series expansion of the inverse projection, depending only on
{@linkplain #eccentricity
+     * eccentricity} value. The series expansion is published under the following form, where
χ is the conformal latitude:
      *
-     *     <blockquote>ci₂⋅f(2θ) + ci₄⋅f(4θ) + ci₆⋅f(6θ) + ci₈⋅f(8θ)</blockquote>
+     *     <blockquote>c₂⋅sin(2χ) + c₄⋅sin(4χ) + c₆⋅sin(6χ) + c₈⋅sin(8χ)</blockquote>
      *
-     * This {@code ConformalProjection} class uses those coefficients in {@link #φ(double)}.
-     * However some subclasses may compute those coefficients differently and use them in
a
-     * different series expansion (but for the same purpose).
+     * However we rewrite above series expansion for taking advantage of trigonometric identities.
+     * The equation become (with different <var>c</var> values):
      *
-     * <p><strong>Consider those fields as final!</strong> They are not
final only for sub-class
-     * constructors convenience and for the purpose of {@link #readObject(ObjectInputStream)}.</p>
+     *     <blockquote>sin(2χ)⋅(c₂ + cos(2χ)⋅(c₄ + cos(2χ)⋅(c₆ + cos(2χ)⋅c₈)))</blockquote>
      *
-     * @see #computeCoefficients()
+     * <div class="note"><b>Serialization note:</b>
+     * we do not strictly need to serialize those fields since they could be computed after
deserialization.
+     * Bu we serialize them anyway in order to simplify a little bit this class (it allows
us to keep those
+     * fields final) and because values computed after deserialization could be slightly
different than the
+     * ones computed after construction if a future version of the constructor uses the double-double
values
+     * provided by {@link Initializer}.</div>
      */
-    transient double ci2, ci4, ci6, ci8;
+    private final double c2χ, c4χ, c6χ, c8χ;
 
     /**
      * Creates a new normalized projection from the parameters computed by the given initializer.
-     *
-     * <p>It is sub-classes responsibility to invoke {@code super.computeCoefficients()}
in their
-     * constructor when ready, or to compute the coefficients themselves.</p>
+     * This constructor computes the coefficients in the series expansion from the {@link
#eccentricitySquared} value.
      *
      * @param  initializer  the initializer for computing map projection internal parameters.
      */
     ConformalProjection(final Initializer initializer) {
         super(initializer);
-    }
-
-    /**
-     * Computes the coefficients in the series expansions from the {@link #eccentricitySquared}
value.
-     * This method shall be invoked after {@code ConformalProjection} construction or deserialization.
-     */
-    void computeCoefficients() {
         useIterations = (eccentricity >= ECCENTRICITY_THRESHOLD);
         final double e2 = eccentricitySquared;
         final double e4 = e2 * e2;
         final double e6 = e2 * e4;
         final double e8 = e4 * e4;
         /*
-         * For each line below, add the smallest values first in order to reduce rounding
errors.
-         * The smallest values are the one using the eccentricity raised to the highest power.
-         */
-        ci2  =    13/   360. * e8  +   1/ 12. * e6  +  5/24. * e4  +  e2/2;
-        ci4  =   811/ 11520. * e8  +  29/240. * e6  +  7/48. * e4;
-        ci6  =    81/  1120. * e8  +   7/120. * e6;
-        ci8  =  4279/161280. * e8;
-        /*
-         * 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 'c' field.
+         * Published equation is of the following form:
+         *
+         *     φ  =  χ + A⋅sin(2χ) + B⋅sin(4χ) + C⋅sin(6χ) + D⋅sin(8χ)    
                 (Snyder 3-34)
+         *
+         * We rewrite as:
+         *
+         *     φ  =  χ + sin(2χ)⋅(A′ + cos(2χ)⋅(B′ + cos(2χ)⋅(C′ + cos(2χ)⋅D′)))
           (Snyder 3-35)
+         *     A′ =  A - C
+         *     B′ =  2B - 4D
+         *     C′ =  4C
+         *     D′ =  8D
+         *
+         * Published coefficients are:
+         *
+         *     A  =    13/360   ⋅ℯ⁸  +   1/12 ⋅ℯ⁶  +  5/24⋅ℯ⁴  +  ℯ²/2
+         *     B  =   811/11520 ⋅ℯ⁸  +  29/240⋅ℯ⁶  +  7/48⋅ℯ⁴
+         *     C  =    81/1120  ⋅ℯ⁸  +   7/120⋅ℯ⁶
+         *     D  =  4279/161280⋅ℯ⁸
+         *
+         * We replace A B C D by A′ B′ C′ D′ using Snyder 3-35. The result are coefficients
below.
+         * For each line, we add the smallest values first in order to reduce rounding errors.
+         * The smallest values are the ones using the eccentricity raised to the highest
power.
          */
-        if (ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            // Multiplication by powers of 2 does not bring any additional rounding error.
-            ci4 *= 2;
-            ci6 *= 4;
-            ci8 *= 8;
-        }
+        c2χ =  -73./ 2016 * e8  +   1./40  * e6  +  5./24 * e4  +  1./2 * e2;
+        c4χ =  233./ 6720 * e8  +  29./120 * e6  +  7./24 * e4;
+        c6χ =   81./  280 * e8  +   7./30  * e6;
+        c8χ = 4279./20160 * e8;
     }
 
     /**
      * Creates a new projection initialized to the values of the given one. This constructor
may be invoked after
-     * we determined that the default implementation can be replaced by an other one, for
example using spherical
+     * we determined that the default implementation can be replaced by another one, for
example using spherical
      * formulas instead than the ellipsoidal ones. This constructor allows to transfer all
parameters to the new
      * instance without recomputing them.
      */
     ConformalProjection(final ConformalProjection other) {
         super(other);
         useIterations = other.useIterations;
-        ci2 = other.ci2;
-        ci4 = other.ci4;
-        ci6 = other.ci6;
-        ci8 = other.ci8;
+        c2χ = other.c2χ;
+        c4χ = other.c4χ;
+        c6χ = other.c6χ;
+        c8χ = other.c8χ;
     }
 
     /**
@@ -228,12 +167,6 @@ abstract class ConformalProjection extends NormalizedProjection {
      * This formula is also part of other projections, since Mercator can be considered as
a special case of
      * Lambert Conic Conformal for instance.
      *
-     * <div class="note"><b>Warning:</b>
-     * this method is valid only if the series expansion coefficients where computed by the
-     * {@link #computeCoefficients()} method defined in this class. This is not the case
of
-     * {@link TransverseMercator} for instance. Note however that even in the later case,
-     * this method can still be seen as a special case of {@code TransverseMercator} formulas.</div>
-     *
      * <p>This function is <em>almost</em> the converse of the {@link #expOfNorthing(double,
double)} function.
      * In a Mercator inverse projection, the value of the {@code expOfSouthing} argument
is {@code exp(-y)}.</p>
      *
@@ -271,24 +204,15 @@ abstract class ConformalProjection extends NormalizedProjection {
         /*
          * Add a correction for the flattened shape of the Earth. The correction can be represented
by an
          * infinite series. Here, we apply only the first 4 terms. Those terms are given
by §1.3.3 in the
-         * EPSG guidance note. Note that we add those terms in reverse order, beginning with
the smallest
-         * values, for reducing rounding errors due to IEEE 754 arithmetic.
+         * EPSG guidance note and modified by application of Snyder 3-35. We add those terms
in reverse order,
+         * beginning with the smallest values, for reducing rounding errors due to IEEE 754
arithmetic.
+         *
+         * For the original formula as published by EPSG, see MercatorMethodComparison.bySeriesExpansion()
+         * in the test suite.
          */
-        if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            φ += ci8 * sin(8*φ)
-               + ci6 * sin(6*φ)
-               + ci4 * sin(4*φ)
-               + ci2 * sin(2*φ);
-        } else {
-            /*
-             * Same formula than above, but rewriten using trigonometric identities in order
to have only two
-             * calls to Math.sin/cos instead than 5. The performance gain is twice faster
on tested machine.
-             */
-            final double sin_2φ = sin(2*φ);
-            final double sin2 = sin_2φ * sin_2φ;
-            φ += ((ci4 + ci8 * (0.50 - sin2)) * cos(2*φ)
-                + (ci2 + ci6 * (0.75 - sin2))) * sin_2φ;
-        }
+        final double sin_2φ = sin(2*φ);
+        final double cos_2φ = cos(2*φ);
+        φ += sin_2φ * (c2χ + cos_2φ * (c4χ + cos_2φ * (c6χ + cos_2φ * c8χ)));
         /*
          * Note: a previous version checked if the value of the smallest term c8χ⋅sin(8φ)
was smaller than
          * the iteration tolerance. But this was not reliable enough. We use now a hard coded
threshold
@@ -402,24 +326,4 @@ abstract class ConformalProjection extends NormalizedProjection {
     final double dy_dφ(final double sinφ, final double cosφ) {
         return (1 / cosφ)  -  eccentricitySquared * cosφ / (1 - eccentricitySquared * (sinφ*sinφ));
     }
-
-    /**
-     * Restores transient fields after deserialization.
-     */
-    private void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException
{
-        in.defaultReadObject();
-        computeCoefficients();
-    }
-
-    /**
-     * Verifies if a trigonometric identity produced the expected value. This method is used
in assertions only,
-     * for values close to the [-1 … +1] range. The tolerance threshold is approximately
1.5E-12 (note that it
-     * still about 7000 time greater than {@code Math.ulp(1.0)}).
-     *
-     * @see #ALLOW_TRIGONOMETRIC_IDENTITIES
-     */
-    static boolean identityEquals(final double actual, final double expected) {
-        return !(abs(actual - expected) >                               // Use !(a >
b) instead of (a <= b) in order to tolerate NaN.
-                (ANGULAR_TOLERANCE / 1000) * max(1, abs(expected)));    // Increase tolerance
for values outside the [-1 … +1] range.
-    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
index 40a22d4..26c6579 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
@@ -222,10 +222,6 @@ abstract class EqualAreaProjection extends NormalizedProjection {
             final double t4β = 0.5 - sin2_β;                                        //
= sin(4β) / ( 4⋅sin(2β))
             final double t8β = (cos2_β - sin2_β)*(cos2_β*cos2_β - cos2_β + 1./8); 
 // = sin(8β) / (32⋅sin(2β))
 
-            assert ConformalProjection.identityEquals(t2β, sin(2*β) / ( 2      ));
-            assert ConformalProjection.identityEquals(t4β, sin(4*β) / ( 8 * t2β));
-            assert ConformalProjection.identityEquals(t8β, sin(8*β) / (64 * t2β));
-
             φ = (ci8*t8β  +  ci4*t4β  +  ci2) * t2β  +  β;
         }
         /*
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
index 08898ac..8796382 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
@@ -210,7 +210,6 @@ public class LambertConicConformal extends ConformalProjection {
     @Workaround(library="JDK", version="1.7")
     private LambertConicConformal(final Initializer initializer) {
         super(initializer);
-        super.computeCoefficients();
         double φ0 = initializer.getAndStore(((initializer.variant & 1) != 0) ?  // Odd
'type' are SP1, even 'type' are SP2.
                 LambertConformal1SP.LATITUDE_OF_ORIGIN : LambertConformal2SP.LATITUDE_OF_FALSE_ORIGIN);
         /*
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
index c1b672f..3a98b54 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
@@ -203,7 +203,6 @@ public class Mercator extends ConformalProjection {
     @Workaround(library="JDK", version="1.7")
     private Mercator(final Initializer initializer) {
         super(initializer);
-        super.computeCoefficients();
         this.variant = initializer.variant;
         /*
          * The "Longitude of natural origin" parameter is found in all Mercator projections
and is mandatory.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
index db20419..81a8010 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
@@ -136,7 +136,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
          *    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
+         *       by:       f(φ) = sin(2φ)⋅(A′ + cos(2φ)⋅(B′ + cos(2φ)⋅(C′
+ cos(2φ)⋅D′)))                   Snyder 3-35
          *       with:       A′ = A - C
          *                   B′ = 2B - 4D
          *                   C′ = 4C
@@ -144,7 +144,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
          *
          *    3) Replacement of  sin2φ  by  2⋅sinφ⋅cosφ  and  cos2φ  by  1 - 2sin²φ:
          *
-         *                 f(φ) = sinφ⋅cosφ⋅(C₁ + sin²φ⋅(C₂ + sin²φ⋅(C₃
+ C₄⋅sin²φ)))
+         *                 f(φ) = sinφ⋅cosφ⋅(C₁ + sin²φ⋅(C₂ + sin²φ⋅(C₃
+ sin²φ⋅C₄)))
          *       with:      C₁ = 2⋅(A′ + B′ + C′ + D′)
          *                  C₂ = 2⋅(−2B′ − 4C′ − 6D′)
          *                  C₃ = 2⋅(4C′ + 12D′)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
index 9f20ecc..3131561 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueMercator.java
@@ -147,7 +147,6 @@ public class ObliqueMercator extends ConformalProjection {
      */
     private ObliqueMercator(final Initializer initializer) {
         super(initializer);
-        super.computeCoefficients();
         final double λc = toRadians(initializer.getAndStore(LONGITUDE_OF_CENTRE));
         final double φc = toRadians(initializer.getAndStore(LATITUDE_OF_CENTRE));
         final double sinφ  = sin(φc);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
index f722e40..eeb987e 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
@@ -140,7 +140,6 @@ public class PolarStereographic extends ConformalProjection {
     @Workaround(library="JDK", version="1.7")
     private PolarStereographic(final Initializer initializer) {
         super(initializer);
-        super.computeCoefficients();
         final byte variant = initializer.variant;
         /*
          * "Standard parallel" and "Latitude of origin" should be mutually exclusive,
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
index 866f993..669b1d3 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
@@ -59,7 +59,7 @@ import static org.apache.sis.internal.referencing.provider.TransverseMercator.*;
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @author  Rémi Maréchal (Geomatys)
- * @version 0.8
+ * @version 1.0
  *
  * @see Mercator
  * @see ObliqueMercator
@@ -67,11 +67,59 @@ import static org.apache.sis.internal.referencing.provider.TransverseMercator.*;
  * @since 0.6
  * @module
  */
-public class TransverseMercator extends ConformalProjection {
+public class TransverseMercator extends NormalizedProjection {
     /**
      * For cross-version compatibility.
      */
-    private static final long serialVersionUID = 8167193470189463354L;
+    private static final long serialVersionUID = -627685138188387835L;
+
+    /**
+     * {@code false} for using the original formulas as published by EPSG, or {@code true}
for using formulas
+     * modified using trigonometric identities. The use of trigonometric identities is for
reducing the amount
+     * of calls to the {@link Math#sin(double)} and similar methods. Some identities used
are:
+     *
+     * <ul>
+     *   <li>sin(2θ) = 2⋅sinθ⋅cosθ</li>
+     *   <li>cos(2θ) = cos²θ - sin²θ</li>
+     *   <li>sin(3θ) = (3 - 4⋅sin²θ)⋅sinθ</li>
+     *   <li>cos(3θ) = (4⋅cos³θ) - 3⋅cosθ</li>
+     *   <li>sin(4θ) = (4 - 8⋅sin²θ)⋅sinθ⋅cosθ</li>
+     *   <li>cos(4θ) = (8⋅cos⁴θ) - (8⋅cos²θ) + 1</li>
+     * </ul>
+     *
+     * Hyperbolic formulas:
+     *
+     * <ul>
+     *   <li>sinh(2θ) = 2⋅sinhθ⋅coshθ</li>
+     *   <li>cosh(2θ) = cosh²θ + sinh²θ   =   2⋅cosh²θ - 1   =   1 + 2⋅sinh²θ</li>
+     *   <li>sinh(3θ) = (3 + 4⋅sinh²θ)⋅sinhθ</li>
+     *   <li>cosh(3θ) = ((4⋅cosh²θ) - 3)⋅coshθ</li>
+     *   <li>sinh(4θ) = (1 + 2⋅sinh²θ)⋅4.sinhθ⋅coshθ
+     *                = 4.cosh(2θ).sinhθ⋅coshθ</li>
+     *   <li>cosh(4θ) = (8⋅cosh⁴θ) - (8⋅cosh²θ) + 1
+     *                = 8⋅cosh²(θ) ⋅ (cosh²θ - 1) + 1
+     *                = 8⋅cosh²(θ) ⋅ sinh²(θ) + 1
+     *                = 2⋅sinh²(2θ) + 1</li>
+     * </ul>
+     *
+     * Note that since this boolean is static final, the compiler should exclude the code
in the branch that is never
+     * executed (no need to comment-out that code).
+     *
+     * @see #identityEquals(double, double)
+     */
+    private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true;
+
+    /**
+     * Verifies if a trigonometric identity produced the expected value. This method is used
in assertions only,
+     * for values close to the [-1 … +1] range. The tolerance threshold is approximately
1.5E-12 (note that it
+     * still about 7000 time greater than {@code Math.ulp(1.0)}).
+     *
+     * @see #ALLOW_TRIGONOMETRIC_IDENTITIES
+     */
+    private static boolean identityEquals(final double actual, final double expected) {
+        return !(abs(actual - expected) >                               // Use !(a >
b) instead of (a <= b) in order to tolerate NaN.
+                (ANGULAR_TOLERANCE / 1000) * max(1, abs(expected)));    // Increase tolerance
for values outside the [-1 … +1] range.
+    }
 
     /**
      * Coefficients in the series expansion of the forward projection,
@@ -83,10 +131,20 @@ public class TransverseMercator extends ConformalProjection {
      * Those coefficients are named h₁, h₂, h₃ and h₄ in §1.3.5.1 of
      * IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015.
      *
-     * <p><strong>Consider those fields as final!</strong>
-     * They are not final only for the purpose of {@link #computeCoefficients()}.</p>
+     * <div class="note"><b>Serialization note:</b>
+     * we do not strictly need to serialize those fields since they could be computed after
deserialization.
+     * Bu we serialize them anyway in order to simplify a little bit this class (it allows
us to keep those
+     * fields final) and because values computed after deserialization could be slightly
different than the
+     * ones computed after construction since the constructor uses the double-double values
provided by
+     * {@link Initializer}.</div>
+     */
+    private final double cf2, cf4, cf6, cf8;
+
+    /**
+     * Coefficients in the series expansion of the inverse projection,
+     * depending only on {@linkplain #eccentricity eccentricity} value.
      */
-    private transient double cf2, cf4, cf6, cf8;
+    private final double ci2, ci4, ci6, ci8;
 
     /**
      * Creates a Transverse Mercator projection from the given parameters.
@@ -136,15 +194,46 @@ public class TransverseMercator extends ConformalProjection {
          * Opportunistically use double-double arithmetic for computation of B since we will
store
          * it in the denormalization matrix, and there is no sine/cosine functions involved
here.
          */
+        double cf4, cf6, cf8, ci4, ci6, ci8;
         final DoubleDouble B;
-        {   // For keeping the 't' variable locale.
+        {   // For keeping the `t` and `n` variables locale.
             /*
-             * EPSG gives:      n  =  f / (2-f)
-             * We rewrite as:   n  =  (1 - b/a) / (1 + b/a)
+             * The n parameter is defined by the EPSG guide as:
+             *
+             *     n = f / (2-f)
+             *
+             * Where f is the flattening factor (1 - b/a). This equation can be rewritten
as:
+             *
+             *     n = (1 - b/a) / (1 + b/a)
+             *
+             * As much as possible, b/a should be computed from the map projection parameters.
+             * However if those parameters are not available anymore, then they can be computed
+             * from the eccentricity as:
+             *
+             *     b/a = √(1 - ℯ²)
              */
             final DoubleDouble t = initializer.axisLengthRatio();   // t  =  b/a
             t.ratio_1m_1p();                                        // t  =  (1 - t) / (1
+ t)
-            computeCoefficients(t.doubleValue());
+            final double n  = t.doubleValue();                      // n = f / (2-f)
+            final double n2 = n  * n;
+            final double n3 = n2 * n;
+            final double n4 = n2 * n2;
+            /*
+             * Coefficients for the forward projections.
+             * Add the smallest values first in order to reduce rounding errors.
+             */
+            cf2 = (   41. /    180)*n4  +  ( 5. /  16)*n3  +  (-2. /  3)*n2  +  n/2;
+            cf4 = (  557. /   1440)*n4  +  (-3. /   5)*n3  +  (13. / 48)*n2;
+            cf6 = ( -103. /    140)*n4  +  (61. / 240)*n3;
+            cf8 = (49561. / 161280)*n4;
+            /*
+             * Coefficients for the inverse projections.
+             * Add the smallest values first in order to reduce rounding errors.
+             */
+            ci2 = (  -1. /    360)*n4  +  (37. /  96)*n3  +  (-2. /  3)*n2  +  n/2;
+            ci4 = (-437. /   1440)*n4  +  ( 1. /  15)*n3  +  ( 1. / 48)*n2;
+            ci6 = ( -37. /    840)*n4  +  (17. / 480)*n3;
+            ci8 = (4397. / 161280)*n4;
             /*
              * Compute B  =  (1 + n²/4 + n⁴/64) / (1 + n)
              */
@@ -203,69 +292,12 @@ public class TransverseMercator extends ConformalProjection {
             cf6 *= 16;   ci6 *= 16;
             cf8 *= 64;   ci8 *= 64;
         }
-    }
-
-    /**
-     * Automatically invoked after deserialization for restoring transient fields.
-     */
-    @Override
-    final void computeCoefficients() {
-        /*
-         * Double-double precision is necessary for computing 'n', otherwise we have rounding
errors
-         * in the 3 last digits. Note that we still have sometime a 1 ULP difference compared
to the
-         * 'n' value at serialization time.
-         */
-        final DoubleDouble t = new DoubleDouble(1d);
-        t.subtract(eccentricitySquared);
-        t.sqrt();
-        t.ratio_1m_1p();
-        computeCoefficients(t.doubleValue());
-        if (ALLOW_TRIGONOMETRIC_IDENTITIES) {
-            // Same scaling than in the constructor.
-            cf4 *=  4;   ci4 *=  4;
-            cf6 *= 16;   ci6 *= 16;
-            cf8 *= 64;   ci8 *= 64;
-        }
-    }
-
-    /**
-     * Computes the transient fields after construction or deserialization.
-     * The <var>n</var> parameter is defined by the EPSG guide as:
-     *
-     *     <blockquote>n = f / (2-f)</blockquote>
-     *
-     * Where <var>f</var> is the flattening factor (1 - b/a). This equation can
be rewritten as:
-     *
-     *     <blockquote>n = (1 - b/a) / (1 + b/a)</blockquote>
-     *
-     * As much as possible, b/a should be computed from the map projection parameters.
-     * However if those parameters are not available anymore, then they can be computed
-     * from the eccentricity as:
-     *
-     *     <blockquote>b/a = √(1 - ℯ²)</blockquote>
-     *
-     * @param  n  the value of {@code f / (2-f)} where {@code f} is the flattening factor.
-     */
-    private void computeCoefficients(final double n) {
-        final double n2 = n  * n;
-        final double n3 = n2 * n;
-        final double n4 = n2 * n2;
-        /*
-         * Coefficients for the forward projections.
-         * Add the smallest values first in order to reduce rounding errors.
-         */
-        cf2 = (   41. /    180)*n4  +  ( 5. /  16)*n3  +  (-2. /  3)*n2  +  n/2;
-        cf4 = (  557. /   1440)*n4  +  (-3. /   5)*n3  +  (13. / 48)*n2;
-        cf6 = ( -103. /    140)*n4  +  (61. / 240)*n3;
-        cf8 = (49561. / 161280)*n4;
-        /*
-         * Coefficients for the inverse projections.
-         * Add the smallest values first in order to reduce rounding errors.
-         */
-        ci2 = (  -1. /    360)*n4  +  (37. /  96)*n3  +  (-2. /  3)*n2  +  n/2;
-        ci4 = (-437. /   1440)*n4  +  ( 1. /  15)*n3  +  ( 1. / 48)*n2;
-        ci6 = ( -37. /    840)*n4  +  (17. / 480)*n3;
-        ci8 = (4397. / 161280)*n4;
+        this.cf4 = cf4;
+        this.cf6 = cf6;
+        this.cf8 = cf8;
+        this.ci4 = ci4;
+        this.ci6 = ci6;
+        this.ci8 = ci8;
     }
 
     /**
@@ -277,6 +309,10 @@ public class TransverseMercator extends ConformalProjection {
         cf4 = other.cf4;
         cf6 = other.cf6;
         cf8 = other.cf8;
+        ci2 = other.ci2;
+        ci4 = other.ci4;
+        ci6 = other.ci6;
+        ci8 = other.ci8;
     }
 
     /**
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
index 010a944..d024bef 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
@@ -36,7 +36,7 @@ import static org.apache.sis.referencing.operation.projection.NormalizedProjecti
  * Tests the {@link ConformalProjection} class.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.6
+ * @version 1.0
  * @since   0.6
  * @module
  */
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
index d02ba59..dd31382 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
@@ -77,7 +77,7 @@ public final class MercatorMethodComparison {   // No 'strictfp' keyword
here si
     private final double c2χ, c4χ, c6χ, c8χ;
 
     /**
-     * Creates a new instance for the excentricty of the WGS84 ellipsoid, which is approximately
0.08181919084262157.
+     * Creates a new instance for the eccentricity of the WGS84 ellipsoid, which is approximately
0.08181919084262157.
      * Reminder: the eccentricity of a sphere is 0.
      */
     public MercatorMethodComparison() {
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java
index ce95d0d..7db509c 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/NoOp.java
@@ -74,7 +74,6 @@ final strictfp class NoOp extends ConformalProjection {
         super(new Initializer(new DefaultOperationMethod(
                 Collections.singletonMap(DefaultOperationMethod.NAME_KEY, parameters.getDescriptor().getName()),
                 2, 2, parameters.getDescriptor()), parameters, Collections.emptyMap(), (byte)
0));
-        super.computeCoefficients();
     }
 
     /**


Mime
View raw message