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:
- *
- * <ul>
- * <li>sin(2⋅x) = 2⋅sin(x)⋅cos(x)</li>
- * <li>sin(3⋅x) = (3 - 4⋅sin²(x))⋅sin(x)</li>
- * <li>sin(4⋅x) = (4 - 8⋅sin²(x))⋅sin(x)⋅cos(x)</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).
- */
- private static final boolean ORIGINAL_FORMULA = false;
-
- /**
* Coefficients in the series expansion used by {@link #φ(double)}.
*
* <p>Consider those fields as final. They are not only of the purpose of {@link
#readObject(ObjectInputStream)}.</p>
@@ -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:
+ *
+ * <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).
+ */
+ 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:
- *
- * <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²θ = 2cosh²θ - 1 = 1 + 2sinh²θ</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).
- */
- 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.
+ *
+ * <p>Consider those fields as final. They are not only of the purpose of {@link
#readObject(ObjectInputStream)}.</p>
*/
- 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 <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 excentricity as:
+ *
+ * <blockquote>b/a = √(1 - ℯ²)</blockquote>
+ *
+ * @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*ξ);
|