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));
}
|