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: + * + *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*ξ);