You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2019/04/25 08:06:26 UTC
[sis] branch geoapi-4.0 updated: Complete the tests for
MeridionalDistance inverse method.
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 7e66753 Complete the tests for MeridionalDistance inverse method.
7e66753 is described below
commit 7e667532207af65ef6eebc4c286b755c72c279f4
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Thu Apr 25 09:58:36 2019 +0200
Complete the tests for MeridionalDistance inverse method.
---
.../operation/projection/Initializer.java | 2 +-
.../projection/MeridionalDistanceBased.java | 19 +++-
.../projection/MeridionalDistanceTest.java | 123 ++++++++++++++-------
3 files changed, 96 insertions(+), 48 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
index e5e8b8b..cb87dd4 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/Initializer.java
@@ -281,7 +281,7 @@ final class Initializer {
* Returns {@code b/a} where {@code a} is the semi-major axis length and {@code b} the semi-minor axis length.
* We retrieve this value from the eccentricity with {@code b/a = sqrt(1-ℯ²)}.
*
- * <p><b>Tip:</b> for ℯ₁ = [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)] (equation Snyder 3-24),
+ * <p><b>Tip:</b> for ℯ₁ = [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)] (Snyder 3-24),
* invoke {@link DoubleDouble#ratio_1m_1p()} on the returned value.</p>
*/
final DoubleDouble axisLengthRatio() {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
index 961e467..5024acc 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceBased.java
@@ -212,7 +212,9 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
* Computes the meridional distance (M) from equator to a given latitude φ.
* Special cases:
* <ul>
- * <li>If <var>φ</var> is 0°, then this method returns 0.</li>
+ * <li>If φ is 0°, then this method returns 0.</li>
+ * <li>If φ=+π/2, then this method returns a value slightly smaller than +π/2, depending on the eccentricity.</li>
+ * <li>If φ=-π/2, then this method returns a value slightly greater than -π/2, depending on the eccentricity.</li>
* </ul>
*
* @param φ latitude for which to compute the meridional distance, in radians.
@@ -221,9 +223,8 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
* @return meridional distance for the given latitude on an ellipsoid of semi-major axis of 1.
*/
final double meridianArc(final double φ, final double sinφ, final double cosφ) {
- final double sinφcosφ = cosφ * sinφ;
- final double sinφ2 = sinφ * sinφ;
- return c0*φ + sinφcosφ*(c1 + sinφ2*(c2 + sinφ2*(c3 + sinφ2*c4))); // TODO: use Math.fma with JDK9.
+ final double sinφ2 = sinφ * sinφ;
+ return c0*φ + cosφ*sinφ*(c1 + sinφ2*(c2 + sinφ2*(c3 + sinφ2*c4))); // TODO: use Math.fma with JDK9.
}
/**
@@ -244,14 +245,20 @@ abstract class MeridionalDistanceBased extends NormalizedProjection {
*
* @param distance meridian distance for which to compute the latitude.
* @return the latitude of given meridian distance, in radians.
- * @throws ProjectionException if the iteration does not converge.
*/
- final double inverse(final double distance) throws ProjectionException {
+ final double latitude(final double distance) {
double φ = distance / rld; // Rectifying latitude µ used as first approximation.
double sinφ = sin(φ);
double cosφ = cos(φ);
double sinφ2 = sinφ * sinφ;
φ += sinφ*cosφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4))); // Snyder 3-26.
+ /*
+ * We could improve accuracy by continuing from here with Newton's iterative method
+ * (see MeridionalDistanceTest.inverse(…) for implementation). However that method
+ * requires calls to meridianArc(double, …), which is itself an approximation based
+ * on series expansion. Consequently the accuracy of iterative method can not be
+ * better than `meridianArc(…)` accuracy.
+ */
return φ;
}
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
index a53ee70..993bfda 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MeridionalDistanceTest.java
@@ -21,6 +21,7 @@ import org.opengis.referencing.operation.MathTransform1D;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.referencing.operation.transform.AbstractMathTransform1D;
import org.apache.sis.referencing.operation.DefaultOperationMethod;
+import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.test.TestUtilities;
import org.apache.sis.test.DependsOnMethod;
import org.apache.sis.test.DependsOn;
@@ -54,27 +55,7 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
private MeridionalDistanceBased create(final boolean ellipsoidal) {
final DefaultOperationMethod provider = new org.apache.sis.internal.referencing.provider.Sinusoidal();
final Sinusoidal projection = new Sinusoidal(provider, parameters(provider, ellipsoidal));
- derivativeDeltas = new double[] {toRadians(0.01)};
tolerance = NormalizedProjection.ANGULAR_TOLERANCE; // = linear tolerance on a sphere of radius 1.
- transform = new AbstractMathTransform1D() {
- @Override public double transform (final double φ) {
- return projection.meridianArc(φ, sin(φ), cos(φ));
- }
- @Override public double derivative(final double φ) {
- final double sinφ = sin(φ);
- return projection.dM_dφ(sinφ*sinφ);
- }
- @Override public MathTransform1D inverse() {
- return new AbstractMathTransform1D() {
- @Override public double transform (final double M) throws TransformException {
- return projection.inverse(M);
- }
- @Override public double derivative(final double φ) throws TransformException {
- throw new TransformException("Unsupported");
- }
- };
- }
- };
return projection;
}
@@ -139,6 +120,43 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
}
/**
+ * Computes the latitude of given meridional distance using Newton's method. This method is described in
+ * <a href="https://en.wikipedia.org/wiki/Meridian_arc#The_inverse_meridian_problem_for_the_ellipsoid">Wikipedia</a>.
+ * This method depends on the forward method and can not be more accurate than that method. For testing purposes, we
+ * use the {@linkplain #referenceMoreAccurate must accurate method implemented in this class} regardless performance.
+ *
+ * @param m the meridional distance on ellipsoidal of semi-major axis length of 1.
+ * @return latitude for the given meridional distance, in radians.
+ */
+ private static double inverse(final MeridionalDistanceBased projection, final double m) throws TransformException {
+ final double ITERATION_TOLERANCE = 1E-13; // Should be less than `referenceMoreAccurate` accuracy.
+ final double e2 = projection.eccentricitySquared;
+ int i = NormalizedProjection.MAXIMUM_ITERATIONS;
+ double φ = m; // Starting point.
+ do {
+ /*
+ * Meridian radius of curvature: M(φ) = a(1 - ℯ²) / (1 - ℯ²sinφ)^(3/2).
+ * We compute the inverse of M.
+ */
+ final double sinφ = sin(φ);
+ double iM = 1 - e2*(sinφ*sinφ);
+ iM *= sqrt(iM);
+ iM /= (1 - e2);
+ /*
+ * From Wikipedia: Newton's method: φᵢ₊₁ = φᵢ - [m(φᵢ) - m] / M(φᵢ)
+ * In principle the numerator should be the derivative of m(φᵢ), but
+ * the meridian radius of curvature can be used instead.
+ */
+ final double t = (m - referenceMoreAccurate(projection, φ)) * iM;
+ φ += t;
+ if (abs(t) <= ITERATION_TOLERANCE) {
+ return φ;
+ }
+ } while (--i >= 0);
+ throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
+ }
+
+ /**
* Compares {@link MeridionalDistanceBased#meridianArc(double, double, double)} with formulas taken as references.
*/
@Test
@@ -172,48 +190,71 @@ public final strictfp class MeridionalDistanceTest extends MapProjectionTestCase
}
/**
- * Tests the {@link MeridionalDistanceBased#inverse(double)} method on an ellipsoid.
+ * Compares {@link MeridionalDistanceBased#latitude(double)} with iterative approach.
+ * On WGS 84 ellipsoid, the meridional distance <var>M</var> is between ±π/2 × 0.9983243 approximately.
*
- * @throws TransformException should never happen.
+ * @throws TransformException if the iterative method does not converge.
*/
@Test
@DependsOnMethod("compareWithReference")
- public void testInverse() throws TransformException {
- isDerivativeSupported = false; // For focusing on inverse transform.
- create(true);
- verifyInDomain(100);
+ public void compareInverse() throws TransformException {
+ final MeridionalDistanceBased projection = create(true);
+ final Random random = TestUtilities.createRandomNumberGenerator();
+ for (int i=0; i<50; i++) {
+ final double m = random.nextDouble() * (PI * 0.998) - (PI/2) * 0.998;
+ final double reference = inverse(projection, m);
+ final double actual = projection.latitude(m);
+ assertEquals("Implementation disagrees with reference.", reference, actual, 1E-10);
+ }
}
/**
- * Tests the {@link MeridionalDistanceBased#dM_dφ(double)} method on a sphere.
+ * Tests all methods, including {@link MeridionalDistanceBased#dM_dφ(double)}, on a sphere.
*
- * @throws TransformException should never happen.
+ * @throws TransformException if an iteration does not converge.
*/
@Test
- @DependsOnMethod("compareWithReference")
- public void testDerivativeOnSphere() throws TransformException {
- isInverseTransformSupported = false; // For focusing on derivative.
- create(false);
- verifyInDomain(20);
+ @DependsOnMethod("compareInverse")
+ public void testOnSphere() throws TransformException {
+ verifyInDomain(false, 20);
}
/**
- * Tests the {@link MeridionalDistanceBased#dM_dφ(double)} method on an ellipsoid.
+ * Tests all methods, including {@link MeridionalDistanceBased#dM_dφ(double)}, on an ellipsoid.
*
- * @throws TransformException should never happen.
+ * @throws TransformException if an iteration does not converge.
*/
@Test
- @DependsOnMethod("testDerivativeOnSphere")
- public void testDerivativeOnEllipsoid() throws TransformException {
- isInverseTransformSupported = false; // For focusing on derivative.
- create(true);
- verifyInDomain(100);
+ @DependsOnMethod("testOnSphere")
+ public void testOnEllipsoid() throws TransformException {
+ verifyInDomain(true, 100);
}
/**
* Verifies transform, inverse transform and derivative in the [-90 … 90]° latitude range.
*/
- private void verifyInDomain(final int numCoordinates) throws TransformException {
+ private void verifyInDomain(final boolean ellipsoidal, final int numCoordinates) throws TransformException {
+ final MeridionalDistanceBased projection = create(ellipsoidal);
+ derivativeDeltas = new double[] {toRadians(0.01)};
+ transform = new AbstractMathTransform1D() {
+ @Override public double transform (final double φ) {
+ return projection.meridianArc(φ, sin(φ), cos(φ));
+ }
+ @Override public double derivative(final double φ) {
+ final double sinφ = sin(φ);
+ return projection.dM_dφ(sinφ*sinφ);
+ }
+ @Override public MathTransform1D inverse() {
+ return new AbstractMathTransform1D() {
+ @Override public double transform (final double M) throws TransformException {
+ return projection.latitude(M);
+ }
+ @Override public double derivative(final double φ) throws TransformException {
+ throw new TransformException("Unsupported");
+ }
+ };
+ }
+ };
verifyInDomain(new double[] {toRadians(-89)},
new double[] {toRadians(+89)},
new int[] {numCoordinates},