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: Document better `GeodeticCalculator` accuracy considerations and replace the special cases in JUnit tests by a `KnownProblem` enumeration in an attempt to make clearer what the problem can be and how the tests handle them.
Date Mon, 18 Jan 2021 13:47:46 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 c65774a  Document better `GeodeticCalculator` accuracy considerations and replace
the special cases in JUnit tests by a `KnownProblem` enumeration in an attempt to make clearer
what the problem can be and how the tests handle them.
c65774a is described below

commit c65774ab7961b2bc0c9322083bb5a492de3118d7
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Jan 18 14:36:20 2021 +0100

    Document better `GeodeticCalculator` accuracy considerations and replace the special cases
in JUnit tests by a `KnownProblem` enumeration in an attempt to make clearer what the problem
can be and how the tests handle them.
---
 .../apache/sis/internal/referencing/Formulas.java  |   6 ++
 .../sis/referencing/GeodesicsOnEllipsoid.java      |  15 ++-
 .../apache/sis/referencing/GeodeticCalculator.java |   4 +-
 .../sis/referencing/GeodesicsOnEllipsoidTest.java  |  47 +++------
 .../sis/referencing/GeodeticCalculatorTest.java    | 105 ++++++++++++---------
 5 files changed, 94 insertions(+), 83 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 ee258e4..80c250a 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
@@ -43,6 +43,12 @@ public final class Formulas extends Static {
      * assuming that the unit of measurement is metre. This constant determines also
      * (indirectly) the minimum accuracy of iterative methods in map projections.
      *
+     * <h4>Maintenance</h4>
+     * If this value is modified, then all usages of this constant should be verified.
+     * Some usages may need to be compensated. For example {@code GeodesicsOnEllipsoid}
+     * uses a millimetric precision by dividing the tolerance by 10 or more. We way want
+     * to keep the same precision there even if {@code LINEAR_TOLERANCE} was made smaller.
+     *
      * @see #ANGULAR_TOLERANCE
      * @see org.apache.sis.internal.util.Numerics#COMPARISON_THRESHOLD
      */
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
index 52fd130..2fb94d0 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodesicsOnEllipsoid.java
@@ -88,16 +88,21 @@ class GeodesicsOnEllipsoid extends GeodeticCalculator {
 
     /**
      * Accuracy threshold for iterative computations, in radians.
-     * We take a finer accuracy than default SIS configuration in order to met the accuracy
of numbers
-     * published in Karney (2013). If this value is modified, the effect can be verified
by executing
-     * the {@code GeodesicsOnEllipsoidTest} methods that compare computed values against
Karney's tables.
-     * Remember to update {@link GeodeticCalculator} class javadoc if this value is changed.
+     * This accuracy must be at least {@value Formulas#ANGULAR_TOLERANCE} degrees (converted
to radians) for
+     * conformance with the accuracy reported in {@link GeodeticCalculator} class javadoc.
Actually we take
+     * a finer accuracy than above value in order to met the accuracy of numbers published
in Karney (2013),
+     * but this extra accuracy is not guaranteed because it is hard to achieve in all cases.
      *
      * <p><b>Note:</b> when the iteration loop detects that it reached
this requested accuracy, the loop
      * completes the iteration step which was in progress. Consequently the final accuracy
is one iteration
      * better than the accuracy computed from this value.</p>
+     *
+     * <h4>Maintenance</h4>
+     * If this value is modified, the effect can be verified by executing the {@code GeodesicsOnEllipsoidTest}
+     * methods that compare computed values against Karney's tables. The {@link GeodeticCalculator}
javadoc may
+     * need to be edited accordingly.
      */
-    static final double ITERATION_TOLERANCE = Formulas.ANGULAR_TOLERANCE * (PI/180) / 20;
+    static final double ITERATION_TOLERANCE = (Formulas.ANGULAR_TOLERANCE / 20) * (PI/180);
 
     /**
      * Difference between ending point and antipode of starting point for considering them
as nearly antipodal.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
index 6161778..06d29a7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/GeodeticCalculator.java
@@ -84,7 +84,7 @@ import static org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEqui
  *   </li>
  * </ul>
  *
- * <h2>Algorithms and accuracy</h2>
+ * <h2>Algorithms</h2>
  * {@code GeodeticCalculator} uses two set of formulas, depending if the figure of the Earth
  * {@linkplain Ellipsoid#isSphere() is a sphere} or an ellipsoid.
  * Publications relevant to this class are:
@@ -104,7 +104,9 @@ import static org.apache.sis.internal.referencing.provider.ModifiedAzimuthalEqui
  *       for the reference implementation.</li>
  * </ul>
  *
+ * <h2>Accuracy</h2>
  * {@code GeodeticCalculator} aims for a positional accuracy of one centimetre.
+ * The accuracy is often better (about one millimetre), but not everywhere.
  * Azimuthal accuracy corresponds to an error of one centimetre at a distance of one kilometer,
  * except for nearly antipodal points (less than 1° of longitude and latitude from antipode)
  * and points close to the poles where the azimuthal errors are larger.
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
index c39d1ec..75ccf4e 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodesicsOnEllipsoidTest.java
@@ -126,7 +126,7 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
          * {@code true} if iteration stopped before to reach the desired accuracy because
of limitation
          * in {@code double} precision. This field must be reset to {@code false} before
any new point.
          *
-         * @see #iterationReachedPrecisionLimit()
+         * @see #relaxIfConfirmed(KnownProblem)
          * @see #clear()
          */
         private boolean iterationReachedPrecisionLimit;
@@ -143,8 +143,7 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
         @Override void store(final String name, final double value) {
             if (name.equals("dα₁ ≪ α₁")) {
                 iterationReachedPrecisionLimit = true;
-            }
-            if (name.length() == 1) {
+            } else if (name.length() == 1) {
                 switch (name.charAt(0)) {
                     case 'x': x = value; break;
                     case 'y': y = value; break;
@@ -523,40 +522,24 @@ public final strictfp class GeodesicsOnEllipsoidTest extends GeodeticCalculatorT
     }
 
     /**
-     * Returns {@code true} if iteration stopped before to reach the desired accuracy because
of limitation in
-     * {@code double} precision. This problem may happen in the {@link GeodesicsOnEllipsoid#computeDistance()}
+     * Returns a value {@literal > 1} if iteration stopped before to reach the desired
accuracy because of limitation
+     * in {@code double} precision. This problem may happen in the {@link GeodesicsOnEllipsoid#computeDistance()}
      * method when {@literal dα₁ ≪ α₁}. If locale variable storage is enabled, this
situation is flagged by the
      * {@code "dα₁ ≪ α₁"} key. Otherwise we conservatively assume that this situation
occurred.
      */
     @Override
-    boolean iterationReachedPrecisionLimit() {
-        if (GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES) {
-            return testedEarth.iterationReachedPrecisionLimit;
-        } else {
-            // Conservative value in absence of information.
-            return true;
-        }
-    }
-
-    /**
-     * Tells whether failure to compute geodesic for the given data should cause the test
case to fail.
-     * This is invoked by {@link #compareAgainstDataset()}.
-     *
-     * @param  expected  a row from the {@code $SIS_DATA/Tests/GeodTest.dat} file.
-     * @return whether the JUnit test should fail.
-     */
-    @Override
-    boolean isFailure(final double[] expected) {
-        final double φ1 = expected[COLUMN_φ1];
-        final double φ2 = expected[COLUMN_φ2];
-        final double Δλ = expected[COLUMN_λ2] - expected[COLUMN_λ1];
-        if (Δλ > 90 && max(abs(φ1), abs(φ2)) < 2E-4) {
-            return true;                            // Ignore equatorial case in previous
version (not anymore).
-        }
-        if (Δλ > 179 && abs(φ1 + φ2) < 0.002) {
-            return false;                           // Ignore antipodal case.
+    double relaxIfConfirmed(final KnownProblem potentialProblem) {
+        if (potentialProblem == KnownProblem.ITERATION_REACHED_PRECISION_LIMIT) {
+            if (GeodesicsOnEllipsoid.STORE_LOCAL_VARIABLES) {
+                if (testedEarth.iterationReachedPrecisionLimit) {
+                    return 2;
+                }
+            } else {
+                // No information about whether the problem really occurred.
+                return 2;
+            }
         }
-        return super.isFailure(expected);
+        return super.relaxIfConfirmed(potentialProblem);
     }
 
     /**
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
index 1df1bb6..e7b62eb 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/GeodeticCalculatorTest.java
@@ -463,24 +463,25 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
                  * We execute only one test for each row instead than executing both tests,
                  * for making sure that `GeodeticCalculator` never see the expected values.
                  */
+                KnownProblem  potentialProblem = null;
                 final boolean isTestingInverse = random.nextBoolean();
-                final double cosφ1 = abs(cos(toRadians(expected[COLUMN_φ1])));        
 // For adjusting longitude tolerance.
-                final double cosφ2 = abs(cos(toRadians(expected[COLUMN_φ2])));
+                final double φ1    = expected[COLUMN_φ1];
+                final double φ2    = expected[COLUMN_φ2];
+                final double cosφ1 = abs(cos(toRadians(φ1)));       // For adjusting longitude
tolerance.
+                final double cosφ2 = abs(cos(toRadians(φ2)));
+                final double Δλ    = abs(expected[COLUMN_λ2] - expected[COLUMN_λ1]);
                 double linearTolerance, latitudeTolerance, longitudeTolerance, azimuthTolerance;
-                final boolean isLongArcOnEquator;
                 if (isSphere) {
                     /*
                      * When spherical formulas are used instead than ellipsoidal formulas,
an error up to 1% is expected
                      * in distance calculations (source: Wikipedia). A yet larger error is
observed for azimuth values,
                      * especially near poles and between antipodal points. Following are
empirical thresholds.
                      */
-                    isLongArcOnEquator = false;                                 // Not a
problem in spherical case.
                     linearTolerance    = expected[COLUMN_Δs] * 0.01;
                     latitudeTolerance  = toDegrees(linearTolerance / c.semiMajorAxis);
-                    longitudeTolerance = expected[COLUMN_φ2] > 89.5 ? 180 : latitudeTolerance
/ cosφ2;
+                    longitudeTolerance = φ2 > 89.5 ? 180 : latitudeTolerance / cosφ2;
                     azimuthTolerance   = 0.5;                                   // About
8.8 metres at distance of 1 km.
                     if (isTestingInverse) {
-                        final double Δλ = abs(expected[COLUMN_λ2] - expected[COLUMN_λ1]);
                              if (Δλ > 179) azimuthTolerance = 100;
                         else if (Δλ > 178) azimuthTolerance = 20;
                         else if (Δλ > 175) azimuthTolerance = 10;
@@ -498,18 +499,16 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
                     longitudeTolerance = Formulas.ANGULAR_TOLERANCE / cosφ2;
                     azimuthTolerance   = Formulas.LINEAR_TOLERANCE * (180/PI) / 10000;
                     if (isTestingInverse) {
-                        final double Δ = max(abs(180 - abs(expected[COLUMN_λ2] - expected[COLUMN_λ1])),
-                                                       abs(expected[COLUMN_φ1] + expected[COLUMN_φ2]));
-                        if (Δ < 1) {
+                        if (max(abs(180 - Δλ), abs(φ1 + φ2)) < 1) {
                             azimuthTolerance = 1 * (180/PI) / 10000;                // 1
meter for 10 km.
                         }
+                        if (Δλ > 90 && max(abs(φ1), abs(φ2)) < 2E-4) {
+                            potentialProblem = KnownProblem.ITERATION_REACHED_PRECISION_LIMIT;
+                        }
+                        if (Δλ > 179 && abs(φ1 + φ2) < 0.002) {
+                            potentialProblem = KnownProblem.NO_CONVERGENCE_ON_ANTIPODAL_POINTS;
+                        }
                     }
-                    /*
-                     * If the start and end points are both within about 3 meters from equator
and
-                     * their distance is more than 10,000 km, we may need to relax the tolerance.
-                     * This is an empirical adjustment based on test failures observation.
-                     */
-                    isLongArcOnEquator = Math.min(cosφ1, cosφ2) >= 0.9999999999999 &&
expected[COLUMN_Δs] > 1E+7;
                     assertEquals("Consistency with accuracy reported in Javadoc.", 0.001,
linearTolerance, 0.0005);
                 }
                 /*
@@ -536,21 +535,24 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
                     assertEquals("λ₂", expected[COLUMN_λ2], end.getOrdinate(1),     
longitudeTolerance);
                     assertEquals("α₂", expected[COLUMN_α2], c.getEndingAzimuth(),   
azimuthTolerance / cosφ2);
                     assertEquals("∆s", expected[COLUMN_Δs], c.getGeodesicDistance(), linearTolerance
*
-                                       (isLongArcOnEquator && iterationReachedPrecisionLimit()
? 2 : 1));
+                                                            relaxIfConfirmed(potentialProblem));
                     clear();
                 } catch (GeodeticException | AssertionError e) {
-                    if (!isTestingInverse || e instanceof AssertionError || isFailure(expected))
{
-                        out.printf("Test failure at line %d: %s%n"
-                                + "The values provided in the test file are:%n"
-                                + "(φ₁,λ₁) = %16.12f %16.12f%n"
-                                + "(φ₂,λ₂) = %16.12f %16.12f%n"
-                                + "The values computed by the geodesic calculator are:%n",
-                                reader.getLineNumber(), e.getLocalizedMessage(),
-                                expected[0], expected[1], expected[3], expected[4]);
-                        out.println(c);
-                        throw e;
+                    if (e instanceof GeodeticException) {
+                        if (potentialProblem == KnownProblem.NO_CONVERGENCE_ON_ANTIPODAL_POINTS)
{
+                            noConvergenceCount++;
+                            continue;
+                        }
                     }
-                    noConvergenceCount++;
+                    out.printf("Test failure at line %d: %s%n"
+                            + "The values provided in the test file are:%n"
+                            + "(φ₁,λ₁) = %16.12f %16.12f%n"
+                            + "(φ₂,λ₂) = %16.12f %16.12f%n"
+                            + "The values computed by the geodesic calculator are:%n",
+                            reader.getLineNumber(), e.getLocalizedMessage(),
+                            expected[0], expected[1], expected[3], expected[4]);
+                    out.println(c);
+                    throw e;
                 }
             }
         }
@@ -558,31 +560,44 @@ public strictfp class GeodeticCalculatorTest extends TestCase {
     }
 
     /**
-     * Returns {@code true} if iteration stopped before to reach the desired accuracy because
the correction
-     * to apply in iterative steps is smaller than what can be applied using IEEE 754 double
arithmetic.
-     * It never happen in the spherical case because there is no iteration in that case.
+     * Returns the factor by which to relax the linear tolerance, or 1 for no relaxation.
+     * This method is invoked after {@link GeodeticCalculator#computeDistance()} has been
invoked.
+     * It should check if the potential problem really occurred, and if yes return a value
greater
+     * than 1. If no problem (other than IEEE 754 rounding errors) is expected to occur,
then this
+     * method should return exactly 1.
      *
-     * <p>This method is invoked only in situations empirically identified as susceptible
to have this problem:</p>
-     * <uL>
-     *   <li>When start point and end point are both within about 3 meters from equator
-     *       and their distance is more than 10,000 km.</li>
-     * </ul>
+     * @param  potentialProblem  the problem that may happen, or {@code null} if none.
+     * @return factor by which to relax the linear tolerance threshold, or 1 if no relaxation.
      */
-    boolean iterationReachedPrecisionLimit() {
-        return false;
+    double relaxIfConfirmed(KnownProblem potentialProblem) {
+        return 1;
     }
 
     /**
-     * Tells whether failure to compute geodesic for the given data should cause the test
case to fail.
-     * The default implementation always return {@code true}. Subclass can override if some
points are
-     * known to fail.
+     * Known problems with our implementation of inverse geodetic with ellipsoidal formulas.
+     * Those problems may occur in {@link #compareAgainstDataset()} test when ellipsoidal
formulas
+     * are used (those problems do not occur with spherical formulas). This is an enumeration
+     * of cases where {@link GeodesicsOnEllipsoid#computeDistance()} does not converge.
      *
-     * @param  expected  a row from the {@code $SIS_DATA/Tests/GeodTest.dat} file.
-     *         Use {@code COLUMN_*} constant for accessing values by column indices.
-     * @return whether the JUnit test should fail.
+     * @see #compareAgainstDataset()
+     * @see #relaxIfConfirmed(KnownProblem)
+     * @see <a href="https://issues.apache.org/jira/browse/SIS-467">SIS-467</a>
      */
-    boolean isFailure(final double[] expected) {
-        return true;
+    enum KnownProblem {
+        /**
+         * Iteration stopped before to reach the desired accuracy because the correction
to apply in iterative
+         * steps is smaller than what can be applied using IEEE 754 double arithmetic. This
problem occurs in
+         * {@literal α₁ -= dα₁} expression when {@literal dα₁ ≪ α₁}, in which
case the subtraction has no effect.
+         * It has been observed on the equator between 2 close points. This is not considered
a failure because
+         * the precision is still in 1 cm precision target. We only need to relax the 1 mm
precision check.
+         */
+        ITERATION_REACHED_PRECISION_LIMIT,
+
+        /**
+         * Iteration failed with a "no convergence error". It sometime happens during distance
calculation between
+         * antipodal points.
+         */
+        NO_CONVERGENCE_ON_ANTIPODAL_POINTS;
     }
 
     /**


Mime
View raw message