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
commit 61aa86927fc6b5e207312c71a490a974f113edcd
Author: Martin Desruisseaux <martin.desruisse...@geomatys.com>
AuthorDate: Sat Jan 29 16:36:14 2022 +0100
Move `fastHypot` in a location where `RotatedPole` can use.
This is mostly for documenting why we use that function.
---
.../apache/sis/internal/referencing/Formulas.java | 48 +++++++++++++++++++++-
.../projection/LambertConicConformal.java | 1 +
.../projection/ModifiedAzimuthalEquidistant.java | 2 +-
.../operation/projection/NormalizedProjection.java | 46 +--------------------
.../operation/projection/ObliqueMercator.java | 1 +
.../operation/projection/ObliqueStereographic.java | 1 +
.../operation/projection/PolarStereographic.java | 1 +
.../operation/projection/package-info.java | 2 +-
.../operation/transform/RotatedPole.java | 3 +-
.../projection/MapProjectionTestCase.java | 2 +-
10 files changed, 57 insertions(+), 50 deletions(-)
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
index 80c250a..ce7408c 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/Formulas.java
@@ -33,7 +33,7 @@ import static
org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE
* do not want to expose publicly those arbitrary values (or at least not in a
too direct way).
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.2
* @since 0.4
* @module
*/
@@ -216,4 +216,50 @@ public final class Formulas extends Static {
*/
return semiMajorAxis / (semiMajorAxis - semiMinorAxis);
}
+
+ /**
+ * 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>
+ *
+ * <h4>Caution</h4>
+ * 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 does not check for such cases;
+ * it is caller responsibility to add this check is necessary, for example
as below:
+ *
+ * {@preformat java
+ * double D = max(fastHypot(x, y), max(abs(x), abs(y)));
+ * }
+ *
+ * According JMH, above check is 1.65 time slower than {@code fastHypot}
without checks.
+ * We define this {@code fastHypot(…)} method for tracing where {@code
sqrt(x² + y²)} is used,
+ * so we can verify if it is used in context where the inaccuracy is
acceptable.
+ *
+ * <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.
+ *
+ * @param x one side of the triangle. Should be approximately in the
[-1 … +1] range.
+ * @param y other side of the triangle. Should be approximately in the
[-1 … +1] range.
+ * @return hypotenuse, not smaller than {@code max(|x|,|y|)} unless the
values are less than 1.5E-154.
+ *
+ * @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>
+ */
+ public static double fastHypot(final double x, final double y) {
+ return sqrt(x*x + y*y);
+ }
}
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 31d7117..037e770 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
@@ -41,6 +41,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.referencing.Formulas.fastHypot;
/**
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 fd62083..121e82c 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
@@ -244,7 +244,7 @@ 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.
*
- * Note: the D ≥ max(|x|,|y|) assumption may not be always true (see
NormalizedProjection.fastHypot(…)).
+ * Note: the D ≥ max(|x|,|y|) assumption may not be always true (see
`Formulas.fastHypot(…)`).
* Consequently sin(α) or cos(α) may be slightly greater than 1.
However they are multiplied by terms
* involving eccentricity, which are smaller than 1. An empirical
verification is done with cos(φ₀) = 1
* in AzimuthalEquidistantTest.testValuesNearZero().
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 3e88c0a..43c7066 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
@@ -127,7 +127,7 @@ import static java.lang.Math.*;
* @author André Gosselin (MPO)
* @author Rueben Schulz (UBC)
* @author Rémi Maréchal (Geomatys)
- * @version 1.1
+ * @version 1.2
*
* @see ContextualParameters
* @see <a href="https://mathworld.wolfram.com/MapProjection.html";>Map
projections on MathWorld</a>
@@ -687,50 +687,6 @@ public abstract class NormalizedProjection extends
AbstractMathTransform2D imple
return θ;
}
- /**
- * 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>
- *
- * <h4>Caution</h4>
- * 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 does not a check for such cases;
- * it is caller responsibility to add this check is necessary, for example
as below:
- *
- * {@preformat java
- * double D = max(fastHypot(x, y), max(abs(x), abs(y)));
- * }
- *
- * According JMH, above check is 1.65 time slower than {@code fastHypot}
without checks.
- * We define this {@code fastHypot(…)} method for tracing where {@code
sqrt(x² + y²)} is used,
- * so we can verify if it is used in context where the inaccuracy is
acceptable.
- *
- * <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 sqrt(x*x + y*y);
- // TODO: consider using Math.fma on JDK9.
- // Check also if we should do the same with plain x*x + y*y in
subclasses.
- }
-
/*
* TODO: consider adding a sqrt1ms(x) method for sqrt(1 - x*x), which
could be implemented as sqrt(fma(x, -x, 1)).
* The use of Math.fma(…) in this context would be valuable especially
when x is close to 1 (to be verified).
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 8080965..c38110c 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
@@ -33,6 +33,7 @@ import org.apache.sis.util.Workaround;
import static java.lang.Math.*;
import static org.apache.sis.math.MathFunctions.atanh;
+import static org.apache.sis.internal.referencing.Formulas.fastHypot;
import static
org.apache.sis.internal.referencing.provider.ObliqueMercatorCenter.*;
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 4839bb5..8a96322 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
@@ -33,6 +33,7 @@ import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.util.Workaround;
import static java.lang.Math.*;
+import static org.apache.sis.internal.referencing.Formulas.fastHypot;
import static
org.apache.sis.internal.referencing.provider.ObliqueStereographic.*;
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 0cd2180..ff765a1 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
@@ -38,6 +38,7 @@ import org.apache.sis.measure.Latitude;
import org.apache.sis.math.MathFunctions;
import static java.lang.Math.*;
+import static org.apache.sis.internal.referencing.Formulas.fastHypot;
/**
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
index 255d2c5..8a0001d 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/package-info.java
@@ -159,7 +159,7 @@
* @author Rémi Maréchal (Geomatys)
* @author Adrian Custer (Geomatys)
* @author Matthieu Bastianelli (Geomatys)
- * @version 1.1
+ * @version 1.2
*
* @see <a href="https://mathworld.wolfram.com/MapProjection.html";>Map
projections on MathWorld</a>
*
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java
index 513d4cd..cd5a385 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/RotatedPole.java
@@ -43,6 +43,7 @@ import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.Debug;
import static java.lang.Math.*;
+import static org.apache.sis.internal.referencing.Formulas.fastHypot;
/**
@@ -311,7 +312,7 @@ public class RotatedPole extends AbstractMathTransform2D
implements Serializable
*/
double xr = cosφp * z - sinφp * x;
double zr = -cosφp * x - sinφp * z;
- double R = sqrt(xr*xr + y*y); // The slower hypot(…) is not
needed because values are close to 1.
+ double R = fastHypot(xr, y); // The slower hypot(…) is not
needed because values are close to 1.
dstPts[dstOff] = atan2(y, xr);
dstPts[dstOff+1] = atan2(zr, R);
if (!derivate) {
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
index 90f73f9..f45342d 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MapProjectionTestCase.java
@@ -238,7 +238,7 @@ abstract strictfp class MapProjectionTestCase extends
MathTransformTestCase {
/**
* Tests coordinates close to zero. Callers must set the transform and
tolerance threshold before to invoke
* this method. This method tests (among others) the
1.4914711209038602E-154 value, which is the threshold
- * documented in {@link NormalizedProjection#fastHypot}.
+ * documented in {@link
org.apache.sis.internal.referencing.Formulas#fastHypot(double, double)}.
*
* @throws TransformException if an error occurred while projecting the
coordinate.
*/
