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: Relax the check for max values in `NormalizedProjection.fastHypot` because it is not always needed. Documentation has been added for warning callers about possible unexpected results. Test cases have been added for the 2 projections where not testing may be at risk.
Date Tue, 14 Jul 2020 16:37:54 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 712bbb1  Relax the check for max values in `NormalizedProjection.fastHypot` because
it is not always needed. Documentation has been added for warning callers about possible unexpected
results. Test cases have been added for the 2 projections where not testing may be at risk.
712bbb1 is described below

commit 712bbb1de39041fb210fa9fad75a5ce609ee808a
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Tue Jul 14 18:34:34 2020 +0200

    Relax the check for max values in `NormalizedProjection.fastHypot` because it is not always
needed.
    Documentation has been added for warning callers about possible unexpected results.
    Test cases have been added for the 2 projections where not testing may be at risk.
---
 .../projection/ModifiedAzimuthalEquidistant.java   | 10 +++++---
 .../operation/projection/NormalizedProjection.java | 17 ++++++++++---
 .../operation/projection/ObliqueStereographic.java | 11 ++++----
 .../projection/AzimuthalEquidistantTest.java       | 29 ++++++++++++++++++++++
 .../projection/MapProjectionTestCase.java          | 27 ++++++++++++++++++++
 .../projection/ObliqueStereographicTest.java       | 18 +++++++++++++-
 6 files changed, 98 insertions(+), 14 deletions(-)

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 6c9c386..652e4f2 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
@@ -220,8 +220,7 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant
{
         final double x  = srcPts[srcOff  ];
         final double y  = srcPts[srcOff+1];
         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.
+        final double D  = sqrt(D2);                     // D = c′/ν₀, but division by
ν₀ is already done here.
         /*
          * From ESPG guidance note:
          *
@@ -231,7 +230,12 @@ public class ModifiedAzimuthalEquidistant extends AzimuthalEquidistant
{
          *
          * But we rewrite in a way that avoid the use of trigonometric functions. We test
(D != 0)
          * 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.
+         * 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(…)).
+         * 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().
          */
         double sinα = 0;
         double cosα = 0;
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 707a648..3e1560c 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
@@ -701,12 +701,21 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D
imple
      * 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.
+     * <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 contains
a check for making sure
-     * that this assumption is true all times.</p>
+     * 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
@@ -717,7 +726,7 @@ public abstract class NormalizedProjection extends AbstractMathTransform2D
imple
      *      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)));
+        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.
     }
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 6b9dec9..4839bb5 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
@@ -461,17 +461,16 @@ public class ObliqueStereographic extends NormalizedProjection {
             final double y = srcPts[srcOff+1];
             final double ρ = fastHypot(x, y);
             final double λ, φ;
-            if (abs(ρ) < ANGULAR_TOLERANCE) {
-                φ = χ0;
+            if (ρ == 0) {         // Exact comparison is okay here. Values > 0 (even
tiny) work with complete formula.
                 λ = 0.0;
+                φ = χ0;
             } else {
                 final double c    = 2*atan(ρ);
                 final double cosc = cos(c);
                 final double sinc = sin(c);
-                final double ct   = ρ * cosχ0*cosc - y*sinχ0*sinc;
-                final double t    = x * sinc;
-                φ = asin(cosc*sinχ0 + y*sinc*cosχ0 / ρ);
-                λ = atan2(t, ct);
+                λ = atan2(x * sinc,
+                          cosc*cosχ0*ρ - y*sinc*sinχ0  );
+                φ = asin( cosc*sinχ0   + y*sinc*cosχ0/ρ);
             }
             dstPts[dstOff]   = λ;
             dstPts[dstOff+1] = φ;
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
index c716b81..a87d5bf 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/AzimuthalEquidistantTest.java
@@ -20,6 +20,7 @@ import org.opengis.util.FactoryException;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.internal.referencing.provider.MapProjection;
 import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.parameter.Parameters;
 import org.apache.sis.test.DependsOn;
 import org.junit.Test;
 
@@ -42,6 +43,19 @@ public strictfp class AzimuthalEquidistantTest extends MapProjectionTestCase
{
     }
 
     /**
+     * Creates a new instance of {@link AzimuthalEquidistant} for a sphere or an ellipsoid.
+     * The new instance is stored in the inherited {@link #transform} field.
+     */
+    private void createNormalizedProjection() {
+        final MapProjection op = method();
+        final Parameters p = Parameters.castOrWrap(op.getParameters().createValue());
+        p.parameter("semi_major").setValue(WGS84_A);
+        p.parameter("semi_minor").setValue(WGS84_B);
+        p.parameter("Latitude of natural origin").setValue(0);
+        transform = new ObliqueStereographic(op, p);
+    }
+
+    /**
      * Tests the forward and inverse projection using test point given in Snyder page 337.
      * The Snyder's test uses a sphere of radius R=3 and a center at 40°N and 100°W.
      * The test in this class modify the longitude to 10°W for avoiding to mix wraparound
@@ -159,4 +173,19 @@ public strictfp class AzimuthalEquidistantTest extends MapProjectionTestCase
{
         verifyDerivative(27, 20);
         verifyDerivative(40, 25);
     }
+
+    /**
+     * Tests {@link AzimuthalEquidistant#inverseTransform(double[], int, double[], int)}
with input
+     * coordinates close to zero. The tested method implementation has an indetermination
at D = 0,
+     * so we test its behavior close to that indetermination point.
+     *
+     * @throws TransformException if an error occurred while projecting the coordinate.
+     */
+    @Test
+    public void testValuesNearZero() throws TransformException {
+        createNormalizedProjection();
+        transform = transform.inverse();
+        tolerance = 1E-15;
+        verifyValuesNearZero(0, 0);
+    }
 }
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 fb61dc8..90f73f9 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
@@ -16,6 +16,7 @@
  */
 package org.apache.sis.referencing.operation.projection;
 
+import java.util.Random;
 import org.opengis.util.FactoryException;
 import org.opengis.referencing.datum.Ellipsoid;
 import org.opengis.referencing.operation.MathTransform;
@@ -30,6 +31,7 @@ import org.apache.sis.referencing.operation.transform.MathTransformTestCase;
 import org.apache.sis.referencing.operation.transform.MathTransformFactoryMock;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
 import org.apache.sis.referencing.datum.GeodeticDatumMock;
+import org.apache.sis.test.TestUtilities;
 
 import static java.lang.Double.isNaN;
 import static java.lang.StrictMath.*;
@@ -232,4 +234,29 @@ abstract strictfp class MapProjectionTestCase extends MathTransformTestCase
{
         }
         verifyInDomain(domain, randomSeed);
     }
+
+    /**
+     * 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}.
+     *
+     * @throws TransformException if an error occurred while projecting the coordinate.
+     */
+    final void verifyValuesNearZero(final double... expected) throws TransformException {
+        final double[] source = new double[2];
+        verifyTransform(source, expected);        // Test (0,0).
+        for (int i=0; i<source.length; i++) {
+            source[i] = 1E-20;                    verifyTransform(source, expected);
+            source[i] = 1E-100;                   verifyTransform(source, expected);
+            source[i] = 1.4914711209038602E-154;  verifyTransform(source, expected);
+            source[i] = 2.2227587494850775E-162;  verifyTransform(source, expected);
+            source[i] = Double.MIN_NORMAL;        verifyTransform(source, expected);
+            source[i] = Double.MIN_VALUE;         verifyTransform(source, expected);
+            source[i] = 0;
+        }
+        final Random r = TestUtilities.createRandomNumberGenerator();
+        for (int i=0; i<25; i++) {
+            source[i & 1] = Math.scalb(r.nextDouble() - 0.5, r.nextInt(800) + Double.MIN_EXPONENT);
+        }
+    }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java
index a4050df..c4ea9bc 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ObliqueStereographicTest.java
@@ -41,7 +41,7 @@ import static org.junit.Assert.*;
  *
  * @author  Rémi Maréchal (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.1
  * @since   0.7
  * @module
  */
@@ -381,4 +381,20 @@ public final strictfp class ObliqueStereographicTest extends MapProjectionTestCa
         tolerance = 0.01;
         verifyTransform(new double[] {44, 73}, new double[] {3320416.75, 632668.43});
     }
+
+    /**
+     * Tests {@link ObliqueStereographic.Spherical#inverseTransform(double[], int, double[],
int)} with input
+     * coordinates close to zero. The tested method implementation has an indetermination
at ρ = 0, so we test
+     * its behavior close to that indetermination point. We test both coordinates, but the
main coordinate of
+     * interest is <var>y</var> because it is used in a formula containing y/ρ
where y and ρ both tend to zero.
+     *
+     * @throws TransformException if an error occurred while projecting the coordinate.
+     */
+    @Test
+    public void testValuesNearZero() throws TransformException {
+        createNormalizedProjection(false);
+        transform = transform.inverse();
+        tolerance = 1E-15;
+        verifyValuesNearZero(0, φ0);
+    }
 }


Mime
View raw message