From commits-return-6416-apmail-sis-commits-archive=sis.apache.org@sis.apache.org Mon Oct 5 20:24:55 2015 Return-Path: X-Original-To: apmail-sis-commits-archive@www.apache.org Delivered-To: apmail-sis-commits-archive@www.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id DA5D218DAA for ; Mon, 5 Oct 2015 20:24:55 +0000 (UTC) Received: (qmail 92440 invoked by uid 500); 5 Oct 2015 20:24:55 -0000 Delivered-To: apmail-sis-commits-archive@sis.apache.org Received: (qmail 92412 invoked by uid 500); 5 Oct 2015 20:24:55 -0000 Mailing-List: contact commits-help@sis.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: sis-dev@sis.apache.org Delivered-To: mailing list commits@sis.apache.org Received: (qmail 92403 invoked by uid 99); 5 Oct 2015 20:24:55 -0000 Received: from Unknown (HELO spamd2-us-west.apache.org) (209.188.14.142) by apache.org (qpsmtpd/0.29) with ESMTP; Mon, 05 Oct 2015 20:24:55 +0000 Received: from localhost (localhost [127.0.0.1]) by spamd2-us-west.apache.org (ASF Mail Server at spamd2-us-west.apache.org) with ESMTP id 56E1C1A0AC9 for ; Mon, 5 Oct 2015 20:24:55 +0000 (UTC) X-Virus-Scanned: Debian amavisd-new at spamd2-us-west.apache.org X-Spam-Flag: NO X-Spam-Score: 1.79 X-Spam-Level: * X-Spam-Status: No, score=1.79 tagged_above=-999 required=6.31 tests=[KAM_ASCII_DIVIDERS=0.8, KAM_LAZY_DOMAIN_SECURITY=1, T_RP_MATCHES_RCVD=-0.01] autolearn=disabled Received: from mx1-us-east.apache.org ([10.40.0.8]) by localhost (spamd2-us-west.apache.org [10.40.0.9]) (amavisd-new, port 10024) with ESMTP id AqNYTKziWAHb for ; Mon, 5 Oct 2015 20:24:52 +0000 (UTC) Received: from mailrelay1-us-west.apache.org (mailrelay1-us-west.apache.org [209.188.14.139]) by mx1-us-east.apache.org (ASF Mail Server at mx1-us-east.apache.org) with ESMTP id BC5A842C69 for ; Mon, 5 Oct 2015 20:24:51 +0000 (UTC) Received: from svn01-us-west.apache.org (svn.apache.org [10.41.0.6]) by mailrelay1-us-west.apache.org (ASF Mail Server at mailrelay1-us-west.apache.org) with ESMTP id 3DEE3E045B for ; Mon, 5 Oct 2015 20:24:51 +0000 (UTC) Received: from svn01-us-west.apache.org (localhost [127.0.0.1]) by svn01-us-west.apache.org (ASF Mail Server at svn01-us-west.apache.org) with ESMTP id 3B1CF3A03C8 for ; Mon, 5 Oct 2015 20:24:51 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 8bit Subject: svn commit: r1706916 - in /sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection: ConformalProjection.java NormalizedProjection.java TransverseMercator.java Date: Mon, 05 Oct 2015 20:24:51 -0000 To: commits@sis.apache.org From: desruisseaux@apache.org X-Mailer: svnmailer-1.0.9 Message-Id: <20151005202451.3B1CF3A03C8@svn01-us-west.apache.org> Author: desruisseaux Date: Mon Oct 5 20:24:50 2015 New Revision: 1706916 URL: http://svn.apache.org/viewvc?rev=1706916&view=rev Log: Tune the application of trigonometric identities in map projections. Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java 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/ConformalProjection.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java?rev=1706916&r1=1706915&r2=1706916&view=diff ============================================================================== --- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java [UTF-8] (original) +++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java [UTF-8] Mon Oct 5 20:24:50 2015 @@ -49,7 +49,7 @@ import static java.lang.Math.*; * * @author Martin Desruisseaux (Geomatys) * @since 0.6 - * @version 0.6 + * @version 0.7 * @module */ abstract class ConformalProjection extends NormalizedProjection { @@ -72,22 +72,6 @@ abstract class ConformalProjection exten static final double EXCENTRICITY_THRESHOLD = 0.16; /** - * Whether to use the original formulas a published by EPSG, or their form modified using trigonometric identities. - * The modified form uses trigonometric identifies for reducing the amount of calls to the {@link Math#sin(double)} - * method. The identities used are: - * - *
    - *
  • sin(2⋅x) = 2⋅sin(x)⋅cos(x)
    • - *
  • sin(3⋅x) = (3 - 4⋅sin²(x))⋅sin(x)
    • - *
  • sin(4⋅x) = (4 - 8⋅sin²(x))⋅sin(x)⋅cos(x)
    • - *
- * - * 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). - */ - private static final boolean ORIGINAL_FORMULA = false; - - /** * Coefficients in the series expansion used by {@link #φ(double)}. * *

Consider those fields as final. They are not only of the purpose of {@link #readObject(ObjectInputStream)}.

@@ -129,7 +113,7 @@ abstract class ConformalProjection exten c4χ = 811/ 11520.* e8 + 29/240.* e6 + 7/48.* e4; c6χ = 81/ 1120.* e8 + 7/120.* e6; c8χ = 4279/161280.* e8; - if (!ORIGINAL_FORMULA) { + if (ALLOW_TRIGONOMETRIC_IDENTITIES) { c4χ *= 2; c6χ *= 4; c8χ *= 8; @@ -196,7 +180,7 @@ abstract class ConformalProjection exten * 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. */ - if (ORIGINAL_FORMULA) { + if (!ALLOW_TRIGONOMETRIC_IDENTITIES) { φ += c8χ * sin(8*φ) + c6χ * sin(6*φ) + c4χ * sin(4*φ) @@ -204,15 +188,12 @@ abstract class ConformalProjection exten } else { /* * Same formula than above, be 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 some machines. + * calls to Math.sin/cos instead than 5. The performance gain is twice faster on tested machine. */ - final double sin2φ = sin(2*φ); - final double sin_cos2φ = cos(2*φ) * sin2φ; - final double sin_sin2φ = sin2φ * sin2φ; - φ += c8χ * (0.50 - sin_sin2φ)*sin_cos2φ // ÷8 compared to original formula - + c6χ * (0.75 - sin_sin2φ)*sin2φ // ÷4 compared to original formula - + c4χ * ( sin_cos2φ) // ÷2 compared to original formula - + c2χ * sin2φ; + final double sin_2φ = sin(2*φ); + final double sin2 = sin_2φ * sin_2φ; + φ += ((c4χ + c8χ * (0.50 - sin2)) * cos(2*φ) + + (c2χ + c6χ * (0.75 - sin2))) * sin_2φ; } /* * Note: a previous version checked if the value of the smallest term c8χ⋅sin(8φ) was smaller than Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java?rev=1706916&r1=1706915&r2=1706916&view=diff ============================================================================== --- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java [UTF-8] (original) +++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java [UTF-8] Mon Oct 5 20:24:50 2015 @@ -116,7 +116,7 @@ import java.util.Objects; * @author Rueben Schulz (UBC) * @author Rémi Maréchal (Geomatys) * @since 0.6 - * @version 0.6 + * @version 0.7 * @module * * @see ContextualParameters @@ -129,6 +129,40 @@ public abstract class NormalizedProjecti private static final long serialVersionUID = 1969740225939106310L; /** + * Whether to use the original formulas a published by EPSG, or their form modified using trigonometric identities. + * The modified form uses trigonometric identifies for reducing the amount of calls to the {@link Math#sin(double)} + * and similar methods. Some identities used are: + * + *
    + *
  • sin(2θ) = 2⋅sinθ⋅cosθ
    • + *
  • cos(2θ) = cos²θ - sin²θ
    • + *
  • sin(3θ) = (3 - 4⋅sin²θ)⋅sinθ
    • + *
  • cos(3θ) = (4⋅cos³θ) - 3⋅cosθ
    • + *
  • sin(4θ) = (4 - 8⋅sin²θ)⋅sinθ⋅cosθ
    • + *
  • cos(4θ) = (8⋅cos⁴θ) - (8⋅cos²θ) + 1
    • + *
+ * + * Hyperbolic formulas: + * + *
    + *
  • sinh(2θ) = 2⋅sinhθ⋅coshθ
    • + *
  • cosh(2θ) = cosh²θ + sinh²θ = 2⋅cosh²θ - 1 = 1 + 2⋅sinh²θ
    • + *
  • sinh(3θ) = (3 + 4⋅sinh²θ)⋅sinhθ
    • + *
  • cosh(3θ) = ((4⋅cosh²θ) - 3)⋅coshθ
    • + *
  • sinh(4θ) = (1 + 2⋅sinh²θ)⋅4.sinhθ⋅coshθ + * = 4.cosh(2θ).sinhθ⋅coshθ
    • + *
  • cosh(4θ) = (8⋅cosh⁴θ) - (8⋅cosh²θ) + 1 + * = 8⋅cosh²(θ) ⋅ (cosh²θ - 1) + 1 + * = 8⋅cosh²(θ) ⋅ sinh²(θ) + 1 + * = 2⋅sinh²(2θ) + 1
    • + *
+ * + * 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). + */ + static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true; + + /** * Maximum difference allowed when comparing longitudes or latitudes in radians. * The current value takes the system-wide angular tolerance value (equivalent to * about 1 cm on Earth) converted to radians. 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=1706916&r1=1706915&r2=1706916&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 20:24:50 2015 @@ -17,6 +17,8 @@ package org.apache.sis.referencing.operation.projection; import java.util.EnumMap; +import java.io.IOException; +import java.io.ObjectInputStream; import org.opengis.util.FactoryException; import org.opengis.parameter.ParameterDescriptor; import org.opengis.referencing.operation.MathTransform; @@ -57,7 +59,7 @@ import static org.apache.sis.math.MathFu * @author Martin Desruisseaux (Geomatys) * @author Rémi Maréchal (Geomatys) * @since 0.6 - * @version 0.6 + * @version 0.7 * @module * * @see Mercator @@ -67,46 +69,15 @@ public class TransverseMercator extends /** * For cross-version compatibility. */ - private static final long serialVersionUID = -4717976245811852528L; - - /** - * Whether to use the original formulas a published by EPSG, or their form modified using trigonometric identities. - * The modified form uses trigonometric identifies for reducing the amount of calls to the {@link Math#sin(double)} - * and similar method. The identities used are: - * - *
    - *
  • sin(2θ) = 2⋅sinθ⋅cosθ
    • - *
  • cos(2θ) = cos²θ - sin²θ
    • - *
  • sin(3θ) = (3 - 4⋅sin²θ)⋅sinθ
    • - *
  • cos(3θ) = (4⋅cos³θ) - 3⋅cosθ
    • - *
  • sin(4θ) = (4 - 8⋅sin²θ)⋅sinθ⋅cosθ
    • - *
  • cos(4θ) = (8⋅cos⁴θ) - (8⋅cos²θ) + 1
    • - *
- * - * Hyperbolic formulas: - *
    - *
  • sinh(2θ) = 2⋅sinhθ⋅coshθ
    • - *
  • cosh(2θ) = cosh²θ + sinh²θ = 2cosh²θ - 1 = 1 + 2sinh²θ
    • - *
  • sinh(3θ) = (3 + 4⋅sinh²θ)⋅sinhθ
    • - *
  • cosh(3θ) = ((4⋅cosh²θ) - 3)coshθ
    • - *
  • sinh(4θ) = (1 + 2⋅sinh²θ)⋅4.sinhθ⋅coshθ - * = 4.cosh(2θ).sinhθ⋅coshθ
    • - *
  • cosh(4θ) = (8⋅cosh⁴θ) - (8⋅cosh²θ) + 1 - * = 8.cosh²θ(cosh²θ - 1) + 1 - * = 8.cosh²(θ).sinh²(θ) + 1 - * = 2.sinh²(2θ) + 1
    • - *
- * - * 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). - */ - private static final boolean ORIGINAL_FORMULA = false; + private static final long serialVersionUID = 8167193470189463354L; /** * Internal coefficients for computation, depending only on value of excentricity. * Defined in §1.3.5.1 of IOGP Publication 373-7-2 – Geomatics Guidance Note number 7, part 2 – April 2015. + * + *

Consider those fields as final. They are not only of the purpose of {@link #readObject(ObjectInputStream)}.

*/ - private final double h1, h2, h3, h4, ih1, ih2, ih3, ih4; + private transient double h1, h2, h3, h4, ih1, ih2, ih3, ih4; /** * Creates a Transverse Mercator projection from the given parameters. @@ -159,7 +130,6 @@ public class TransverseMercator extends * 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. */ - final double n; final DoubleDouble B; { // For keeping the 't' variable locale. /* @@ -168,35 +138,16 @@ public class TransverseMercator extends */ final DoubleDouble t = initializer.axisLengthRatio(); // t = b/a t.ratio_1m_1p(); // t = (1 - t) / (1 + t) - n = t.doubleValue(); + initialize(t.doubleValue()); /* * Compute B = (1 + n²/4 + n⁴/64) / (1 + n) */ B = new DoubleDouble(t); // B = n B.square(); B.series(1, 0.25, 1./64); // B = (1 + n²/4 + n⁴/64) - t.add(1,0); + t.add(1, 0); B.divide(t); // B = (1 + n²/4 + n⁴/64) / (1 + n) } - final double n2 = n * n; - final double n3 = n2 * n; - final double n4 = n2 * n2; - /* - * Coefficients for direct projection. - * Add the smallest values first in order to reduce rounding errors. - */ - h1 = ( 41. / 180)*n4 + ( 5. / 16)*n3 + (-2. / 3)*n2 + n/2; - h2 = ( 557. / 1440)*n4 + (-3. / 5)*n3 + (13. / 48)*n2; - h3 = ( -103. / 140)*n4 + (61. / 240)*n3; - h4 = (49561. / 161280)*n4; - /* - * Coefficients for inverse projection. - * Add the smallest values first in order to reduce rounding errors. - */ - ih1 = ( -1. / 360)*n4 + (37. / 96)*n3 + (-2. / 3)*n2 + n/2; - ih2 = (-437. / 1440)*n4 + ( 1. / 15)*n3 + ( 1. / 48)*n2; - ih3 = ( -37. / 840)*n4 + (17. / 480)*n3; - ih4 = (4397. / 161280)*n4; /* * Compute M₀ = B⋅(ξ₁ + ξ₂ + ξ₃ + ξ₄) and negate in anticipation for what will be needed * in the denormalization matrix. We opportunistically use double-double arithmetic but @@ -237,6 +188,63 @@ public class TransverseMercator extends } /** + * Restores transient fields after deserialization. + */ + private void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException { + in.defaultReadObject(); + /* + * 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(1, 0); + t.subtract(excentricitySquared, 0); + t.sqrt(); + t.ratio_1m_1p(); + initialize(t.doubleValue()); + } + + /** + * Computes the transient fields after construction or deserialization. + * 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 excentricity as: + * + *
b/a = √(1 - ℯ²)
+ * + * @param n The value of {@code f / (2-f)} where {@code f} is the flattening factor. + */ + private void initialize(final double n) { + final double n2 = n * n; + final double n3 = n2 * n; + final double n4 = n2 * n2; + /* + * Coefficients for direct projection. + * Add the smallest values first in order to reduce rounding errors. + */ + h1 = ( 41. / 180)*n4 + ( 5. / 16)*n3 + (-2. / 3)*n2 + n/2; + h2 = ( 557. / 1440)*n4 + (-3. / 5)*n3 + (13. / 48)*n2; + h3 = ( -103. / 140)*n4 + (61. / 240)*n3; + h4 = (49561. / 161280)*n4; + /* + * Coefficients for inverse projection. + * Add the smallest values first in order to reduce rounding errors. + */ + ih1 = ( -1. / 360)*n4 + (37. / 96)*n3 + (-2. / 3)*n2 + n/2; + ih2 = (-437. / 1440)*n4 + ( 1. / 15)*n3 + ( 1. / 48)*n2; + ih3 = ( -37. / 840)*n4 + (17. / 480)*n3; + ih4 = (4397. / 161280)*n4; + } + + /** * Creates a new projection initialized to the same parameters than the given one. */ TransverseMercator(final TransverseMercator other) { @@ -313,7 +321,7 @@ public class TransverseMercator extends */ 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; - if (ORIGINAL_FORMULA) { + if (!ALLOW_TRIGONOMETRIC_IDENTITIES) { sin_2ξ0 = sin(2*ξ0); cos_2ξ0 = cos(2*ξ0); sin_4ξ0 = sin(4*ξ0); @@ -342,7 +350,7 @@ public class TransverseMercator extends */ 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; - if (ORIGINAL_FORMULA) { + if (!ALLOW_TRIGONOMETRIC_IDENTITIES) { sinh_2η0 = sinh(2*η0); cosh_2η0 = cosh(2*η0); sinh_4η0 = sinh(4*η0); @@ -474,7 +482,7 @@ public class TransverseMercator extends cos_2ξ, cos_4ξ, cos_6ξ, cos_8ξ, sinh_2η, sinh_4η, sinh_6η, sinh_8η, cosh_2η, cosh_4η, cosh_6η, cosh_8η; - if (ORIGINAL_FORMULA) { + if (!ALLOW_TRIGONOMETRIC_IDENTITIES) { sin_2ξ = sin(2*ξ); cos_2ξ = cos(2*ξ); sin_4ξ = sin(4*ξ);