sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject svn commit: r1693564 - in /sis/branches/JDK8/core: sis-referencing/src/main/java/org/apache/sis/referencing/datum/ sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ sis-referencing/src/test/java/org/apache/sis/referencing/o...
Date Fri, 31 Jul 2015 10:42:33 GMT
Author: desruisseaux
Date: Fri Jul 31 10:42:33 2015
New Revision: 1693564

URL: http://svn.apache.org/r1693564
Log:
Partial rollback of the use of double-double arithmetic in map projection initialization.
Our usage of double-double arithmetic has proven its value in matrix operations, but has
less value in NormalizedProjection subclasses after the point where we use transcendental
functions (sine, logarithmic, etc.) because we have no double-double versions of those functions.
By reducing double-double arithmetic usage in those cases, we keep the code more readable
and avoid to give a false sense of accuracy.

Modified:
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
    sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/TransverseMercator.java
    sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
    sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/TimeDependentBWP.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -23,6 +23,7 @@ import org.apache.sis.internal.util.Nume
 import org.apache.sis.internal.util.DoubleDouble;
 
 import static org.apache.sis.util.ArgumentChecks.*;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
 import static org.apache.sis.internal.referencing.Formulas.JULIAN_YEAR_LENGTH;
 
 
@@ -159,7 +160,7 @@ public class TimeDependentBWP extends Bu
         if (time != null) {
             final long millis = time.getTime() - timeReference;
             if (millis != 0) { // Returns null for 0 as an optimization.
-                final DoubleDouble period = new DoubleDouble(millis, 0);
+                final DoubleDouble period = verbatim(millis);
                 period.divide(1000 * JULIAN_YEAR_LENGTH, 0);
                 return period;
             }

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -28,14 +28,20 @@ import org.apache.sis.referencing.operat
 
 import static java.lang.Math.*;
 import static org.apache.sis.util.ArgumentChecks.ensureNonNull;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
 
 
 /**
  * Helper class for map projection constructions, providing formulas normally needed only
at construction time.
- * Since map projection constructions should not happen very often, we afford using double-double
arithmetic here.
+ * Since map projection constructions should not happen very often, we afford using some
double-double arithmetic.
  * The main intend is not to provide more accurate coordinate conversions (while it may be
a nice side-effect),
- * but rather to increase the chances that the concatenations of (de)normalization matrices
with the matrices of
- * other transforms give back identity matrices when such result is expected.
+ * but to improve the result of concatenations of (de)normalization matrices with the matrices
of other transforms,
+ * as found in transformation chains.
+ *
+ * <p>As a general rule, we stop storing result with double-double precision after
the point where we need
+ * transcendental functions (sine, logarithm, <i>etc.</i>), since we do not have
double-double versions of
+ * those functions. Digits after the {@code double} part are usually not significant in such
cases, except
+ * in some relatively rare scenarios like 1 ± (a result much smaller than 1).</p>
  *
  * @author  Martin Desruisseaux (Geomatys)
  * @since   0.6
@@ -71,7 +77,13 @@ final class Initializer {
     final byte variant;
 
     /**
-     * Creates a new initializer.
+     * Creates a new initializer. The parameters are described in
+     * {@link NormalizedProjection#NormalizedProjection(OperationMethod, Parameters, Map)}.
+     *
+     * @param method     Description of the map projection parameters.
+     * @param parameters The parameters of the projection to be created.
+     * @param roles Parameters to look for <cite>central meridian</cite>, <cite>scale
factor</cite>,
+     *        <cite>false easting</cite>, <cite>false northing</cite>
and other values.
      */
     Initializer(final OperationMethod method, final Parameters parameters,
             final Map<ParameterRole, ? extends ParameterDescriptor<Double>> roles,
@@ -102,25 +114,23 @@ final class Initializer {
                         - getAndStore(roles.get(ParameterRole.FALSE_SOUTHING));
 
         excentricitySquared = new DoubleDouble();
-        final DoubleDouble k = new DoubleDouble(a);  // The value by which to multiply all
results of normalized projection.
+        DoubleDouble k = new DoubleDouble(a);  // The value by which to multiply all results
of normalized projection.
         if (a != b) {
             /*
-             * Equivalent Java code for the following lines:
-             *
-             *     final double rs = b / a;
-             *     excentricitySquared = 1 - (rs * rs);
+             * ℯ² = 1 - (b/a)²
              *
-             * Test show that double-double arithmetic here makes a difference in the 3 last
digits for WGS84 ellipsoid.
-             * Those 3 digits are not significant since the parameter are not so accurate
(furthermore the 'b' parameter
-             * used below may have been computed from the inverse flattening factor).
+             * Double-double arithmetic here makes a difference in the 3 last digits for
WGS84 ellipsoid.
              */
-            final DoubleDouble rs = new DoubleDouble(b);
-            final double eb = rs.error;
-            rs.divide(k);    // rs = b/a
-            rs.multiply(rs);
-            excentricitySquared.value = 1;
-            excentricitySquared.subtract(rs);
-
+            if (DoubleDouble.DISABLED) {
+                final double rs = b / a;
+                excentricitySquared.value = 1 - (rs * rs);
+            } else {
+                final DoubleDouble rs = new DoubleDouble(b);
+                rs.divide(k);    // rs = b/a
+                rs.multiply(rs);
+                excentricitySquared.value = 1;
+                excentricitySquared.subtract(rs);
+            }
             final ParameterDescriptor<Double> radius = roles.get(ParameterRole.LATITUDE_OF_CONFORMAL_SPHERE_RADIUS);
             if (radius != null) {
                 /*
@@ -139,13 +149,8 @@ final class Initializer {
                  *     final double sinφ = sin(toRadians(parameters.doubleValue(radius)));
                  *     k = b / (1 - excentricitySquared * (sinφ*sinφ));
                  */
-                final DoubleDouble t = new DoubleDouble(sin(toRadians(parameters.doubleValue(radius))),
0);
-                t.multiply(t);
-                t.multiply(excentricitySquared);
-                k.clear();
-                k.value = 1;
-                k.subtract(t);
-                k.inverseDivide(b, eb);
+                k = rν2(sin(toRadians(parameters.doubleValue(radius))));
+                k.inverseDivide(b, 0);
             }
         }
         context.normalizeGeographicInputs(λ0);
@@ -198,28 +203,19 @@ final class Initializer {
         return value;
     }
 
-
-
-
-    //////////////////////////////////////////////////////////////////////////////////////////
-    ////////                                                                          ////////
-    ////////                       FORMULAS FROM EPSG or SNYDER                       ////////
-    ////////                                                                          ////////
-    //////////////////////////////////////////////////////////////////////////////////////////
-
     /**
-     * Computes the reciprocal of the radius of curvature of the ellipsoid perpendicular
to the meridian at latitude φ.
-     * That radius of curvature is:
+     * Computes the square of the reciprocal of the radius of curvature of the ellipsoid
+     * perpendicular to the meridian at latitude φ. That radius of curvature is:
      *
      * <blockquote>ν = 1 / √(1 - ℯ²⋅sin²φ)</blockquote>
      *
-     * This method returns 1/ν.
+     * This method returns 1/ν².
      *
      * <div class="section">Relationship with Snyder</div>
      * This is related to functions (14-15) from Snyder (used for computation of scale factors
      * at the true scale latitude) as below:
      *
-     * <blockquote>m = cosφ / rν</blockquote>
+     * <blockquote>m = cosφ / sqrt(rν²)</blockquote>
      *
      * Special cases:
      * <ul>
@@ -228,24 +224,46 @@ final class Initializer {
      *       (otherwise we get {@link Double#NaN}).</li>
      * </ul>
      *
-     * @param  sinφ The sine of the φ latitude in radians.
-     * @return Reciprocal of the radius of curvature of the ellipsoid perpendicular to the
meridian at latitude φ.
-     */
-    final DoubleDouble rν(final double sinφ) {
-        /*
-         * Equivalent Java code:
-         *
-         *     return sqrt(1 - excentricitySquared * (sinφ*sinφ));
-         */
-        final DoubleDouble t = new DoubleDouble(sinφ, 0);
+     * @param  sinφ The sine of the φ latitude.
+     * @return Reciprocal squared of the radius of curvature of the ellipsoid
+     *         perpendicular to the meridian at latitude φ.
+     */
+    private DoubleDouble rν2(final double sinφ) {
+        if (DoubleDouble.DISABLED) {
+            return verbatim(1 - excentricitySquared.value * (sinφ*sinφ));
+        }
+        final DoubleDouble t = verbatim(sinφ);
         t.multiply(t);
         t.multiply(excentricitySquared);
+
+        // Save result of ℯ²⋅sin²φ
         final double value = t.value;
         final double error = t.error;
+
+        // Compute 1 - above. Since above may be small, this
+        // is where double-double arithmetic has more value.
         t.clear();
         t.value = 1;
         t.subtract(value, error);
-        t.sqrt();
         return t;
     }
+
+    /**
+     * Returns the scale factor at latitude φ. This is computed as:
+     *
+     * <blockquote>cosφ / sqrt(rν2(sinφ))</blockquote>
+     *
+     * The result is returned as a {@code double} because the limited precision of {@code
sinφ} and {@code cosφ}
+     * makes the error term meaningless. We use double-double arithmetic only for intermediate
calculation.
+     *
+     * @param  sinφ The sine of the φ latitude.
+     * @param  cosφ The cosine of the φ latitude.
+     * @return Scale factor at latitude φ.
+     */
+    final double scaleAtφ(final double sinφ, final double cosφ) {
+        final DoubleDouble s = rν2(sinφ);
+        s.sqrt();
+        s.inverseDivide(cosφ, 0);
+        return s.value;
+    }
 }

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertConicConformal.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -40,6 +40,7 @@ import org.apache.sis.util.Workaround;
 import static java.lang.Math.*;
 import static java.lang.Double.*;
 import static org.apache.sis.math.MathFunctions.isPositive;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
 
 
 /**
@@ -96,8 +97,20 @@ public class LambertConicConformal exten
 
     /**
      * Constant for the Belgium 2SP case. This is 29.2985 seconds, given here in radians.
+     *
+     * <div class="note"><b>Tip:</b> how to verify the value:
+     * {@preformat java
+     *     BigDecimal a = new BigDecimal(BELGE_A.value);
+     *     a = a.add     (new BigDecimal(BELGE_A.error));
+     *     a = a.multiply(new BigDecimal("57.29577951308232087679815481410517"));
+     *     a = a.multiply(new BigDecimal(60 * 60));
+     *     System.out.println(a);
+     * }
+     * </div>
      */
-    private static final double BELGE_A = 0.00014204313635987700;
+    static Number belgeA() {
+        return new DoubleDouble(-1.420431363598774E-4, -1.1777378450498224E-20);
+    }
 
     /**
      * Internal coefficients for computation, depending only on values of standards parallels.
@@ -243,62 +256,35 @@ public class LambertConicConformal exten
         φ1 = toRadians(φ1);
         φ2 = toRadians(φ2);
         /*
-         * Computes constants. We do not need to use special formulas for the spherical case
below,
+         * Compute constants. We do not need to use special formulas for the spherical case
below,
          * since rν(sinφ) = 1 and expOfNorthing(φ) = tan(π/4 + φ/2) when the excentricity
is zero.
          * However we need special formulas for φ1 ≈ φ2 in the calculation of n, otherwise
we got
          * a 0/0 indetermination.
-         *
-         * Opportunistically use double-double arithmetic below since this is what we will
store in the
-         * (de)normalization matrices. The extra precision that we get is not necessarily
significant,
-         * but we do that more in an attempt to reduce rounding errors in concatenations
of a sequence
-         * of MathTransforms (through matrix multiplications) than for map projection precisions.
-         * Equivalent Java code for the following double-double arithmetic:
-         *
-         *     final double m1 = cos(φ1) / rν(sinφ1);
          */
         final double sinφ1 = sin(φ1);
-        final DoubleDouble m1 = initializer.rν(sinφ1);
-        m1.inverseDivide(cos(φ1), 0);
+        final double m1 = initializer.scaleAtφ(sinφ1, cos(φ1));
         final double t1 = expOfNorthing(φ1, excentricity*sinφ1);
         /*
-         * Computes n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂), which we rewrite as
ln(m₁/m₂) / ln(t₁/t₂)
+         * Compute n = (ln m₁ – ln m₂) / (ln t₁ – ln t₂), which we rewrite as
ln(m₁/m₂) / ln(t₁/t₂)
          * since division is less at risk of precision lost than subtraction. Note that this
equation
          * tends toward 0/0 if φ₁ ≈ φ₂, which force us to do a special check for
the SP1 case.
-         *
-         * Equivalent Java code for the following double-double arithmetic:
-         *
-         *     final double sinφ2 = sin(φ2);
-         *     final double m2 = cos(φ2) / rν(sinφ2);
-         *     final double t2 = expOfNorthing(φ2, excentricity*sinφ2);
-         *     n = log(m1/m2) / log(t1/t2);    // Tend toward 0/0 if φ1 ≈ φ2.
          */
-        final DoubleDouble F = new DoubleDouble();
         if (abs(φ1 - φ2) >= ANGULAR_TOLERANCE) {  // Should be 'true' for 2SP case.
             final double sinφ2 = sin(φ2);
-            final DoubleDouble m2 = initializer.rν(sinφ2);
-            m2.inverseDivide(cos(φ2), 0);
+            final double m2 = initializer.scaleAtφ(sinφ2, cos(φ2));
             final double t2 = expOfNorthing(φ2, excentricity*sinφ2);
-            m2.inverseDivide(m1);
-            F.value = log(m2.value);
-            F.divide(log(t1/t2), 0);
+            n = log(m1/m2) / log(t1/t2);    // Tend toward 0/0 if φ1 ≈ φ2.
         } else {
-            F.value = -sinφ1;
+            n = -sinφ1;
         }
-        n = F.value;
         /*
-         * Scale factor for longitudes, stored now before we modify the F value.
-         */
-        final DoubleDouble sx = new DoubleDouble(F);
-        if (!isNorth) {
-            sx.negate();
-        }
-        /*
-         * Computes F = m₁/(n⋅t₁ⁿ) from Geomatics Guidance Note number 7.
+         * Compute F = m₁/(n⋅t₁ⁿ) from Geomatics Guidance Note number 7.
          * Following constants will be stored in the denormalization matrix, to be applied
          * after the non-linear formulas implemented by this LambertConicConformal class.
          * Opportunistically use double-double arithmetic since the matrix coefficients will
          * be stored in that format anyway. This makes a change in the 2 or 3 last digits.
          */
+        final DoubleDouble F = verbatim(n);
         F.multiply(pow(t1, n), 0);
         F.inverseDivide(m1);
         if (!isNorth) {
@@ -306,16 +292,16 @@ public class LambertConicConformal exten
         }
         /*
          * Compute the radius of the parallel of latitude of the false origin.
-         * This is related to the "ρ0" term in Snyder. From EPG guide:
+         * This is related to the "ρ₀" term in Snyder. From EPG guide:
          *
          *    r = a⋅F⋅tⁿ     where (in our case) a=1 and t is our 'expOfNorthing' function.
          *
          * EPSG uses this term in the computation of  y = FN + rF – r⋅cos(θ).
          */
-        final DoubleDouble rF = new DoubleDouble();    // Initialized to zero.
+        DoubleDouble rF = null;
         if (φ0 != copySign(PI/2, -n)) {    // For reducing the rounding error documented
in expOfNorthing(+π/2).
-            rF.value = pow(expOfNorthing(φ0, excentricity*sin(φ0)), n);
-            rF.multiply(F);
+            rF = new DoubleDouble(F);
+            rF.multiply(pow(expOfNorthing(φ0, excentricity*sin(φ0)), n), 0);
         }
         /*
          * At this point, all parameters have been processed. Now store
@@ -334,14 +320,19 @@ public class LambertConicConformal exten
          *   - Multiply by the scale factor (done by the super-class constructor).
          *   - Add false easting and false northing (done by the super-class constructor).
          */
-        final MatrixSIS normalize = context.getMatrix(true);
-        normalize.convertAfter(0, sx, (initializer.variant == BELGIUM) ? new DoubleDouble(-BELGE_A,
0) : null);
+        DoubleDouble sλ = verbatim(n);
+        DoubleDouble sφ = null;
         if (isNorth) {
-            normalize.convertAfter(1, new DoubleDouble(-1, 0), null);
+            // Reverse the sign of either longitude or latitude values before map projection.
+            sφ = verbatim(-1);
+        } else {
+            sλ.negate();
         }
+        final MatrixSIS normalize   = context.getMatrix(true);
         final MatrixSIS denormalize = context.getMatrix(false);
-        denormalize.convertBefore(0, F, null);
-        F.negate();
+        normalize  .convertAfter(0, sλ, (initializer.variant == BELGIUM) ? belgeA() : null);
+        normalize  .convertAfter(1, sφ, null);
+        denormalize.convertBefore(0, F, null); F.negate();
         denormalize.convertBefore(1, F, rF);
     }
 

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Mercator.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -38,6 +38,7 @@ import org.apache.sis.util.Workaround;
 import static java.lang.Math.*;
 import static java.lang.Double.*;
 import static org.apache.sis.math.MathFunctions.isPositive;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
 
 
 /**
@@ -228,8 +229,7 @@ public class Mercator extends ConformalP
          * if they really want, since we sometime see such CRS definitions.
          */
         final double φ1 = toRadians(initializer.getAndStore(Mercator2SP.STANDARD_PARALLEL));
-        final DoubleDouble k0 = new DoubleDouble(cos(φ1), 0);
-        k0.divide(initializer.rν(sin(φ1)));
+        final Number k0 = verbatim(initializer.scaleAtφ(sin(φ1), cos(φ1)));
         /*
          * In principle we should rotate the central meridian (λ0) in the normalization
transform, as below:
          *
@@ -251,11 +251,11 @@ public class Mercator extends ConformalP
             denormalize.convertBefore(0, null, offset);
         }
         if (φ0 != 0) {
-            denormalize.convertBefore(1, null, new DoubleDouble(-log(expOfNorthing(φ0, excentricity
* sin(φ0)))));
+            denormalize.convertBefore(1, null, verbatim(-log(expOfNorthing(φ0, excentricity
* sin(φ0)))));
         }
         if (variant == MILLER) {
-              normalize.convertBefore(1, new DoubleDouble(0.80), null);
-            denormalize.convertBefore(1, new DoubleDouble(1.25), null);
+            normalize  .convertBefore(1, 0.80, null);
+            denormalize.convertBefore(1, 1.25, null);
         }
         /*
          * At this point we are done, but we add here a little bit a maniac precision hunting.
@@ -281,9 +281,9 @@ public class Mercator extends ConformalP
          * those remaning lines of code.
          */
         if (φ0 == 0 && isPositive(φ1 != 0 ? φ1 : φ0)) {
-            final DoubleDouble revert = new DoubleDouble(-1, 0);
-              normalize.convertBefore(1, revert, null);
-            denormalize.convertBefore(1, revert, null);  // Must be before false easting/northing.
+            final Number reverseSign = verbatim(-1);
+            normalize  .convertBefore(1, reverseSign, null);
+            denormalize.convertBefore(1, reverseSign, null);  // Must be before false easting/northing.
         }
     }
 

Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/PolarStereographic.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -29,7 +29,6 @@ import org.apache.sis.internal.referenci
 import org.apache.sis.internal.referencing.provider.PolarStereographicB;
 import org.apache.sis.internal.referencing.provider.PolarStereographicC;
 import org.apache.sis.internal.referencing.Formulas;
-import org.apache.sis.internal.util.DoubleDouble;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.Workaround;
@@ -37,6 +36,7 @@ import org.apache.sis.measure.Latitude;
 import org.apache.sis.math.MathFunctions;
 
 import static java.lang.Math.*;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
 
 
 /**
@@ -185,8 +185,7 @@ public class PolarStereographic extends
             φ1 = -φ1;
         }
         φ1 = toRadians(φ1);  // May be anything in [-π/2 … 0] range.
-        final DoubleDouble ρ;
-        DoubleDouble ρF = null;    // Actually -ρF (compared to EPSG guide).
+        final Number ρ, ρF;  // This ρF is actually -ρF in EPSG guide.
         if (abs(φ1 + PI/2) < ANGULAR_TOLERANCE) {
             /*
              * Polar Stereographic (variant A)
@@ -199,19 +198,9 @@ public class PolarStereographic extends
              *    - the 't' factor, because it needs to be computed in the transform(…)
method.
              *
              * In the spherical case, should give ρ == 2.
-             *
-             * Opportunistically use double-double arithmetic below since this is what we
will store in the
-             * (de)normalization matrices. The extra precision that we get is not necessarily
significant,
-             * but we do that more in an attempt to reduce rounding errors in concatenations
of a sequence
-             * of MathTransforms (through matrix multiplications) than for map projection
precisions.
-             * Equivalent Java code for the following double-double arithmetic:
-             *
-             *     ρ = 2 / sqrt(pow(1+excentricity, 1+excentricity) * pow(1-excentricity,
1-excentricity));
              */
-            ρ = new DoubleDouble(pow(1+excentricity, 1+excentricity), 0);
-            ρ.multiply          (pow(1-excentricity, 1-excentricity), 0);
-            ρ.sqrt();
-            ρ.inverseDivide(2, 0);
+            ρ = verbatim(2 / sqrt(pow(1+excentricity, 1+excentricity) * pow(1-excentricity,
1-excentricity)));
+            ρF = null;
         } else {
             /*
              * Polar Stereographic (variant B or C)
@@ -229,21 +218,11 @@ public class PolarStereographic extends
              *   k₀ = ρ⋅√[…]/2  but we do not need that value.
              *
              * In the spherical case, should give ρ = 1 + sinφ1   (Synder 21-7 and 21-11).
-             *
-             * Equivalent Java code for the following double-double arithmetic:
-             *
-             *     final double mF = cos(φ1) / rν(sinφ1);
-             *     ρ = mF / expOfNorthing(φ1, excentricity*sinφ1);
-             *     if (variant == C) ρF = -mF;
              */
             final double sinφ1 = sin(φ1);
-            ρ = initializer.rν(sinφ1);
-            ρ.inverseDivide(cos(φ1), 0);
-            if (variant == C) {
-                ρF = new DoubleDouble(ρ);
-                ρF.negate();
-            }
-            ρ.divide(expOfNorthing(φ1, excentricity*sinφ1), 0);
+            final double mF = initializer.scaleAtφ(sinφ1, cos(φ1));
+            ρ = verbatim(mF / expOfNorthing(φ1, excentricity*sinφ1));
+            ρF = (variant == C) ? verbatim(-mF) : null;
         }
         /*
          * At this point, all parameters have been processed. Now process to their
@@ -253,8 +232,10 @@ public class PolarStereographic extends
         denormalize.convertBefore(0, ρ, null);
         denormalize.convertBefore(1, ρ, ρF);
         if (isNorth) {
-            context.getMatrix(true).convertAfter(1, -1, null);
-            denormalize.convertBefore(1, -1, null);
+            final Number reverseSign = verbatim(-1);
+            final MatrixSIS normalize = context.getMatrix(true);
+            normalize  .convertAfter (1, reverseSign, null);
+            denormalize.convertBefore(1, reverseSign, null);
         }
     }
 

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=1693564&r1=1693563&r2=1693564&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] Fri Jul 31 10:42:33 2015
@@ -31,6 +31,7 @@ import org.apache.sis.util.Workaround;
 import static java.lang.Math.*;
 import static org.apache.sis.math.MathFunctions.asinh;
 import static org.apache.sis.math.MathFunctions.atanh;
+import static org.apache.sis.internal.util.DoubleDouble.verbatim;
 
 
 /**
@@ -129,7 +130,7 @@ public class TransverseMercator extends
          * Compute B  =  (n4/64 + n2/4 + 1) / (n + 1)
          * Opportunistically uses double-double arithmetic since we use it anyway for denormalization
matrix.
          */
-        final DoubleDouble B = new DoubleDouble(n);
+        final DoubleDouble B = verbatim(n);
         B.add(1);
         B.inverseDivide(1, n4/64 + n2/4);
         /*
@@ -160,7 +161,7 @@ public class TransverseMercator extends
          */
         final double Q = asinh(tan(φ0)) - excentricity * atanh(excentricity * sin(φ0));
         final double β = atan(sinh(Q));
-        final DoubleDouble M0 = new DoubleDouble(β, 0);
+        final DoubleDouble M0 = verbatim(β);
         M0.add(h1 * sin(2*β), 0);
         M0.add(h2 * sin(4*β), 0);
         M0.add(h3 * sin(6*β), 0);

Modified: sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/LambertConicConformalTest.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -17,6 +17,7 @@
 package org.apache.sis.referencing.operation.projection;
 
 import java.util.Random;
+import java.math.BigDecimal;
 import org.opengis.util.FactoryException;
 import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.TransformException;
@@ -28,6 +29,7 @@ import org.apache.sis.internal.referenci
 import org.apache.sis.internal.referencing.provider.LambertConformalMichigan;
 import org.apache.sis.referencing.operation.transform.CoordinateDomain;
 import org.apache.sis.parameter.Parameters;
+import org.apache.sis.internal.util.DoubleDouble;
 import org.apache.sis.test.DependsOnMethod;
 import org.apache.sis.test.DependsOn;
 import org.apache.sis.test.TestUtilities;
@@ -53,6 +55,22 @@ import static org.junit.Assert.*;
 @DependsOn(ConformalProjectionTest.class)
 public final strictfp class LambertConicConformalTest extends MapProjectionTestCase {
     /**
+     * Verifies the value of the constant used in <cite>"Lambert Conic Conformal (2SP
Belgium)"</cite> projection.
+     *
+     * @see #testLambertConicConformalBelgium()
+     */
+    @Test
+    public void verifyBelgeConstant() {
+        final DoubleDouble BELGE_A = (DoubleDouble) LambertConicConformal.belgeA();
+        BigDecimal a = new BigDecimal(BELGE_A.value);
+        a = a.add     (new BigDecimal(BELGE_A.error));
+        a = a.multiply(new BigDecimal("57.29577951308232087679815481410517"));  // Conversion
from radians to degrees.
+        a = a.multiply(new BigDecimal(60 * 60));                                // Conversion
from degrees to seconds.
+        a = a.add     (new BigDecimal("29.2985"));                              // The standard
value.
+        assertTrue(Math.abs(a.doubleValue()) < 1E-31);
+    }
+
+    /**
      * Creates a new instance of {@link LambertConicConformal}. See the class javadoc for
an explanation
      * about why we ask only for the latitude of origin and not the standard parallels.
      *
@@ -199,7 +217,7 @@ public final strictfp class LambertConic
      * @see org.opengis.test.referencing.ParameterizedTransformTest#testLambertConicConformal1SP()
      */
     @Test
-    @DependsOnMethod("testLambertConicConformal2SP")
+    @DependsOnMethod({"testLambertConicConformal2SP", "verifyBelgeConstant"})
     public void testLambertConicConformalBelgium() throws FactoryException, TransformException
{
         createGeoApiTest(new LambertConformalBelgium()).testLambertConicConformalBelgium();
     }

Modified: sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java?rev=1693564&r1=1693563&r2=1693564&view=diff
==============================================================================
--- sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] (original)
+++ sis/branches/JDK8/core/sis-utility/src/main/java/org/apache/sis/internal/util/DoubleDouble.java
[UTF-8] Fri Jul 31 10:42:33 2015
@@ -237,6 +237,20 @@ public final class DoubleDouble extends
     }
 
     /**
+     * Uses the given value verbatim, without inferring an error term for double-double arithmetic.
+     * We use this method when the value has been computed using transcendental functions
(cosine,
+     * logarithm, <i>etc.</i>) in which case there is no way we can infer a meaningful
error term.
+     *
+     * <p>We use this method both for readability and for making easier to search where
such thing occur.</p>
+     *
+     * @param  value The value to wrap in a {@code DoubleDouble} instance.
+     * @return A {@code DoubleDouble} containing exactly the given value, without error term.
+     */
+    public static DoubleDouble verbatim(final double value) {
+        return new DoubleDouble(value, 0);
+    }
+
+    /**
      * Returns a new {@code DoubleDouble} instance initialized to the conversion factor
      * from radians to angular degrees.
      *



Mime
View raw message