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: Replace some calls to `Math.hypot` by `NormalizedProjection.fastHypot`, which is faster at the cost of accuracy. The accuracy penality seems to be within 1 ULP when x and y are in [-6 … +6] range (tested with random numbers). An analysis is available at https://issues.apache.org/jira/browse/NUMBERS-143
Date Mon, 13 Jul 2020 14:48:57 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 dd447f2  Replace some calls to `Math.hypot` by `NormalizedProjection.fastHypot`,
which is faster at the cost of accuracy. The accuracy penality seems to be within 1 ULP when
x and y are in [-6 … +6] range (tested with random numbers). An analysis is available at
https://issues.apache.org/jira/browse/NUMBERS-143
dd447f2 is described below

commit dd447f29b716bc5a69e0852c78622d48606474ef
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Jul 13 16:44:35 2020 +0200

    Replace some calls to `Math.hypot` by `NormalizedProjection.fastHypot`, which is faster
at the cost of accuracy.
    The accuracy penality seems to be within 1 ULP when x and y are in [-6 … +6] range (tested
with random numbers).
    An analysis is available at https://issues.apache.org/jira/browse/NUMBERS-143
---
 .../projection/LambertConicConformal.java          |  4 +--
 .../projection/ModifiedAzimuthalEquidistant.java   | 13 ++++----
 .../operation/projection/NormalizedProjection.java | 35 ++++++++++++++++++++++
 .../operation/projection/ObliqueMercator.java      | 13 ++++----
 .../operation/projection/ObliqueStereographic.java |  2 +-
 .../operation/projection/PolarStereographic.java   |  5 ++--
 6 files changed, 55 insertions(+), 17 deletions(-)

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 976a787..31d7117 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
@@ -473,8 +473,8 @@ public class LambertConicConformal extends ConformalProjection {
          * applied before the first non-linear one moved to the inverse of the "denormalize"
transform, and the
          * linear operations applied after the last non-linear one moved to the inverse of
the "normalize" transform.
          */
-        dstPts[dstOff  ] = atan2(x, y);                 // Really (x,y), not (y,x)
-        dstPts[dstOff+1] = -φ(pow(hypot(x, y), 1/n));   // Equivalent to φ(pow(hypot(x,y),
-1/n)) but more accurate for n>0.
+        dstPts[dstOff  ] = atan2(x, y);                     // Really (x,y), not (y,x)
+        dstPts[dstOff+1] = -φ(pow(fastHypot(x, y), 1/n));   // Equivalent to φ(pow(hypot(x,y),
-1/n)) but more accurate for n>0.
     }
 
 
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
index a3e3f18..6c9c386 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ModifiedAzimuthalEquidistant.java
@@ -219,8 +219,9 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant
{
     {
         final double x  = srcPts[srcOff  ];
         final double y  = srcPts[srcOff+1];
-        final double D  = hypot(x, y);                  // D = c′/ν₀, but division by
ν₀ is already done here.
-        final double D2 = D*D;
+        final double D2 = x*x + y*y;
+        final double D  = max(sqrt(D2), max(abs(x), abs(y)));       // See `NormalizedProjection.fastHypot(…)`.
+        // D = c′/ν₀, but division by ν₀ is already done here.
         /*
          * From ESPG guidance note:
          *
@@ -232,11 +233,11 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant
{
          * exactly (without epsilon) because even a very small value is sufficient for avoiding
NaN:
          * Since D ≥ max(|x|, |y|) we get x/D and y/D close to zero.
          */
-        double sinα = x;                                // x and y interchanged compared
to usual atan2(y, x).
-        double cosα = y;
+        double sinα = 0;
+        double cosα = 0;
         if (D != 0) {
-            sinα /= D;
-            cosα /= D;
+            sinα = x / D;                               // x and y interchanged compared
to usual atan2(y, x).
+            cosα = y / D;
         }
               double negA = Hp * cosα; negA *= negA;    // negA = −A  compared to EPSG
guidance note.
         final double B    = Bp * (1 + negA) * cosα;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
index 2f25267..707a648 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/NormalizedProjection.java
@@ -688,6 +688,41 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D
imple
     }
 
     /**
+     * Returns {@code sqrt(x² + y²)} for coordinate values on an ellipsoid of semi-major
axis length of 1.
+     * This method does not provides the accuracy guarantees offered by {@link Math#hypot(double,
double)}.
+     * However for values close to 1, this approximation seems to stay within 1 ULP of {@code
Math.hypot(…)}.
+     * We tested with random values in ranges up to [-6 … +6].
+     *
+     * <p>We define this method because {@link Math#hypot(double, double)} has been
measured with JMH as 6 times
+     * slower than {@link Math#sqrt(double)} on Java 14.  According posts on internet, the
same performance cost
+     * is observed in C/C++ too. Despite its cost, {@code hypot(…)} is generally recommended
because computing a
+     * hypotenuse from large magnitudes has accuracy problems. But in the context of {@code
NormalizedProjection}
+     * where semi-axis lengths are close to 1, input values should be (x,y) coordinates in
the [−1 … +1] range.
+     * The actual range may be greater (e.g. [−5 … +5]), but it still far from ranges
requiring protection against
+     * overflow.</p>
+     *
+     * <p>We may not need the full {@code Math.hypot(x,y)} accuracy in the context
of map projections on ellipsoids.
+     * However some projection formulas require that {@code fastHypot(x,y) ≥ max(|x|,|y|)},
otherwise normalizations
+     * such as {@code x/hypot(x,y)} could result in values larger than 1, which in turn result
in {@link Double#NaN}
+     * when given to {@link Math#asin(double)}. The assumption on x, y and {@code sqrt(x²+y²)}
relative magnitude is
+     * broken when x=0 and |y| ≤ 1.4914711209038602E-154 or conversely. This method contains
a check for making sure
+     * that this assumption is true all times.</p>
+     *
+     * <h4>When to use</h4>
+     * We reserve this method to ellipsoidal formulas where approximations are used anyway.
Implementations using
+     * exact formulas, such as spherical formulas, should use {@link Math#hypot(double, double)}
for its accuracy.
+     *
+     * @see <a href="https://issues.apache.org/jira/browse/NUMBERS-143">Investigate
Math.hypot for computing the absolute of a complex number</a>
+     * @see <a href="https://scicomp.stackexchange.com/questions/27758/is-there-any-point-to-using-hypot-for-sqrt1c2-0-le-c-le-1-for-real/27766">Is
+     *      there any point to using <code>hypot(1, c)</code> for <code>sqrt(1
+ c²)</code>, 0 ≤ c ≤ 1</a>
+     */
+    static double fastHypot(final double x, final double y) {
+        return max(sqrt(x*x + y*y), max(abs(x), abs(y)));
+        // TODO: consider using Math.fma on JDK9.
+        // Check also if we should do the same with plain x*x + y*y in subclasses.
+    }
+
+    /**
      * Converts a single coordinate in {@code srcPts} at the given offset and stores the
result
      * in {@code dstPts} at the given offset. In addition, opportunistically computes the
      * transform derivative if requested.
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 7a85dac..5f5fa34 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
@@ -353,12 +353,13 @@ public class ObliqueMercator extends ConformalProjection {
         final double dU_dλ = -B * (cosγ0 / T) * dV_dλ;
         final double dU_dφ = dQ_dφ * (sinγ0 + (sinγ0 + U)/(Q*Q) - U) / (2*T);
         final double dS_dφ = 0.5*dQ_dφ * (1 + 1/(Q*Q));
-        final double M = (S*cosγ0 + V*sinγ0);
-        final double L = hypot(dV_dλ, M);
-        final double P = L + dV_dλ;
-        final double D = (P*P + M*M);
-        final double dy_dλ = 2 * B * (dV_dλ * (sinγ0*P + (V - sinγ0*M) * M/L) + V*M)
/ D;      // Y = atan2(M, T)
-        final double dy_dφ = 2 * cosγ0 * dS_dφ * (P - M*M/L) / D;
+        final double M   = (S*cosγ0 + V*sinγ0);
+        final double L   = fastHypot(dV_dλ, M);
+        final double M_L = M / L;
+        final double P   = L + dV_dλ;
+        final double D   = (P*P + M*M);
+        final double dy_dλ = 2 * B * (dV_dλ * (sinγ0*P + (V - sinγ0*M) * M_L) + V*M)
/ D;    // Y = atan2(M, T)
+        final double dy_dφ = 2 * cosγ0 * dS_dφ * (P - M*M_L) / D;
         final double R = U*U - 1;
         return new Matrix2(dU_dλ/R, dU_dφ/R,
                            dy_dλ,   dy_dφ);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
index 77eae93..6b9dec9 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ObliqueStereographic.java
@@ -459,7 +459,7 @@ public class ObliqueStereographic extends NormalizedProjection {
         {
             final double x = srcPts[srcOff  ];
             final double y = srcPts[srcOff+1];
-            final double ρ = hypot(x, y);
+            final double ρ = fastHypot(x, y);
             final double λ, φ;
             if (abs(ρ) < ANGULAR_TOLERANCE) {
                 φ = χ0;
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 de89f86..0cd2180 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
@@ -206,7 +206,8 @@ public class PolarStereographic extends ConformalProjection {
              *
              * In the spherical case, should give ρ == 2.
              */
-            ρ = new DoubleDouble(2 / sqrt(pow(1+eccentricity, 1+eccentricity) * pow(1-eccentricity,
1-eccentricity)));
+            ρ = new DoubleDouble(2 / sqrt(pow(1+eccentricity, 1+eccentricity)
+                                        * pow(1-eccentricity, 1-eccentricity)));
             ρF = null;
         } else {
             /*
@@ -334,7 +335,7 @@ public class PolarStereographic extends ConformalProjection {
         final double x = srcPts[srcOff  ];
         final double y = srcPts[srcOff+1];
         dstPts[dstOff  ] = atan2(x, y);         // Really (x,y), not (y,x)
-        dstPts[dstOff+1] = -φ(hypot(x, y));
+        dstPts[dstOff+1] = -φ(fastHypot(x, y));
     }
 
 


Mime
View raw message