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/29 16:34:21 UTC
[sis] branch geoapi-4.0 updated (851bd10 -> 5b271a0)
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a change to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git.
from 851bd10 Apply to ConformalProjection the same Snyder 3-35 formulas (trigonometric identities) than we did for MeridianArcBased. Include the coefficients in the serialization, since recomputing values after deserialization may not produce the exact same results than the values computed at construction time. Since we apply a Snyder formula now, those fields shouls be more stable than before.
new 3ec1609 Cleanup in test and benchmark methods.
new 6760894 Apply to EqualAreaProjection the same Snyder 3-35 formulas (trigonometric identities) than we did for MeridianArcBased. Include coefficients and other fields in serialization for the same reason than ConformalProjection (simplicity, stability). During this work, we identified a bug in the series expansion of φ(β) where the last term was multiplied by e⁸ instead of e⁶. Fixing that bug bring the series expansion to a precision where we no longer need iterations at least wi [...]
new 5b271a0 Include MeridianArcBased fields in serialization for the same simplicity and numerical stability reasons than ConformalProjection and EqualAreaProjection. Inconvenient is that the serialization is more at risk to break if the formulas implemented in those case are improved. However MathTransform serialization should not be used for long-term storage. It is okay for RMI purpose between JVMs running the same version of SIS.
The 3 revisions listed above as "new" are entirely new to this
repository and will be described in separate emails. The revisions
listed as "add" were already present in the repository and have only
been added to this reference.
Summary of changes:
.../operation/projection/AlbersEqualArea.java | 1 -
.../operation/projection/ConformalProjection.java | 14 +-
.../operation/projection/CylindricalEqualArea.java | 1 -
.../operation/projection/EqualAreaProjection.java | 202 ++++++++++-----------
.../operation/projection/MeridianArcBased.java | 109 +++++------
.../projection/ConformalProjectionTest.java | 134 ++++++--------
.../projection/EqualAreaProjectionTest.java | 127 +++++++++++++
.../projection/MercatorMethodComparison.java | 107 +++--------
.../sis/test/suite/ReferencingTestSuite.java | 1 +
9 files changed, 360 insertions(+), 336 deletions(-)
create mode 100644 core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EqualAreaProjectionTest.java
[sis] 01/03: Cleanup in test and benchmark methods.
Posted by de...@apache.org.
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 3ec160968286ea19a72eff75412e23f202586504
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Mon Apr 29 10:48:11 2019 +0200
Cleanup in test and benchmark methods.
---
.../projection/ConformalProjectionTest.java | 134 +++++++++------------
.../projection/MercatorMethodComparison.java | 107 +++++-----------
2 files changed, 89 insertions(+), 152 deletions(-)
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
index d024bef..651a031 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/ConformalProjectionTest.java
@@ -60,7 +60,7 @@ public final strictfp class ConformalProjectionTest extends TransformTestCase {
* }
*/
@Test
- public void testMath() {
+ public void verifyMath() {
assertEquals("Forward 0°N", 0, log(tan(PI/4)), TOLERANCE);
assertEquals("Inverse 0 m", 0, PI/2 - 2*atan(exp(0)), TOLERANCE);
assertEquals("Forward 90°S", NEGATIVE_INFINITY, log(tan(0)), TOLERANCE);
@@ -88,10 +88,29 @@ public final strictfp class ConformalProjectionTest extends TransformTestCase {
}
/**
- * Implementation of {@link #testExpOfNorthing()}.
- * The {@link #transform} field must have been set before this method is called.
+ * Tests the {@link ConformalProjection#expOfNorthing(double, double)} function.
+ *
+ * {@preformat text
+ * Forward: y = -log(t(φ))
+ * Inverse: φ = φ(exp(-y))
+ * }
+ */
+ @Test
+ @DependsOnMethod("verifyMath")
+ public void verifyNorthingKnownValues() {
+ verifyNorthingKnownValues(false); // Spherical case
+ verifyNorthingKnownValues(true); // Ellipsoidal case
+ }
+
+ /**
+ * Tests {@link ConformalProjection#expOfNorthing(double, double)} on some known values.
+ *
+ * @param ellipsoidal {@code true} for an ellipsoidal case, or {@code false} for a spherical case.
*/
- private void doTestExpOfNorthing() {
+ private void verifyNorthingKnownValues(final boolean ellipsoidal) {
+ transform = new NoOp(ellipsoidal);
+ tolerance = TOLERANCE;
+
assertEquals("f(NaN) = NaN", NaN, expOfNorthing(NaN), tolerance);
assertEquals("f( ±∞) = NaN", NaN, expOfNorthing(NEGATIVE_INFINITY), tolerance);
assertEquals("f( ±∞) = NaN", NaN, expOfNorthing(POSITIVE_INFINITY), tolerance);
@@ -135,91 +154,64 @@ public final strictfp class ConformalProjectionTest extends TransformTestCase {
/**
* Computes {@link ConformalProjection#expOfNorthing(double, double)} for the given latitude.
+ * The {@link #transform} field must have been set before this method is invoked.
*
* @param φ the latitude in radians.
* @return {@code Math.exp} of the Mercator projection of the given latitude.
- */
- private double expOfNorthing(final double φ) {
- return expOfNorthing((ConformalProjection) transform, φ);
- }
-
- /**
- * Computes {@link ConformalProjection#expOfNorthing(double, double)} for the given latitude.
*
- * @param projection the projection on which to invoke {@code expOfNorthing(…)}.
- * @param φ the latitude in radians.
- * @return {@code Math.exp} of the Mercator projection of the given latitude.
+ * @see #φ(double)
*/
- private static double expOfNorthing(final ConformalProjection projection, final double φ) {
+ private double expOfNorthing(final double φ) {
+ final ConformalProjection projection = (ConformalProjection) transform;
return projection.expOfNorthing(φ, projection.eccentricity * sin(φ));
}
/**
- * Tests the {@link ConformalProjection#expOfNorthing(double, double)} function.
- *
- * {@preformat text
- * Forward: y = -log(t(φ))
- * Inverse: φ = φ(exp(-y))
- * }
- */
- @Test
- @DependsOnMethod("testMath")
- public void testExpOfNorthing() {
- transform = new NoOp(false); // Spherical case
- tolerance = TOLERANCE;
- doTestExpOfNorthing();
- transform = new NoOp(true); // Ellipsoidal case
- doTestExpOfNorthing();
- }
-
- /**
* Tests the {@link ConformalProjection#dy_dφ(double, double)} method.
*
* @throws TransformException if an error occurred while projecting a point.
*/
@Test
- @DependsOnMethod("testExpOfNorthing")
+ @DependsOnMethod("verifyNorthingKnownValues")
public void test_dy_dφ() throws TransformException {
- tolerance = 1E-7;
- doTest_dy_dφ(new NoOp(false)); // Spherical case
- doTest_dy_dφ(new NoOp(true)); // Ellipsoidal case
+ test_dy_dφ(false); // Spherical case
+ test_dy_dφ(true); // Ellipsoidal case
}
/**
- * Implementation of {@link #test_dy_dφ()}.
+ * Tests the {@link ConformalProjection#dy_dφ(double, double)} method.
+ *
+ * @param ellipsoidal {@code true} for an ellipsoidal case, or {@code false} for a spherical case.
*/
- private void doTest_dy_dφ(final NoOp projection) throws TransformException {
+ private void test_dy_dφ(final boolean ellipsoidal) throws TransformException {
+ final NoOp projection = new NoOp(ellipsoidal);
transform = new AbstractMathTransform1D() {
@Override public double transform(final double φ) {
- return expOfNorthing(projection, φ);
+ return projection.expOfNorthing(φ, projection.eccentricity * sin(φ));
}
@Override public double derivative(final double φ) {
final double sinφ = sin(φ);
- return projection.dy_dφ(sinφ, cos(φ)) * expOfNorthing(projection, φ);
+ return projection.dy_dφ(sinφ, cos(φ)) * transform(φ);
}
};
- verifyInDomain(-89 * (PI/180), 89 * (PI/180)); // Verify from 85°S to 85°N.
- }
-
- /**
- * Convenience method invoking {@link TransformTestCase#verifyInDomain} for an 1D transform.
- */
- private void verifyInDomain(final double min, final double max) throws TransformException {
- derivativeDeltas = new double[] {2E-8};
isInverseTransformSupported = false;
- verifyInDomain(
- new double[] {min}, // Minimal value to test.
- new double[] {max}, // Maximal value to test.
- new int[] {100}, // Number of points to test.
- TestUtilities.createRandomNumberGenerator());
+ derivativeDeltas = new double[] {2E-8};
+ tolerance = 1E-7;
+ verifyInDomain(new double[] {-89 * (PI/180)}, // Minimal value to test.
+ new double[] {+89 * (PI/180)}, // Maximal value to test.
+ new int[] {100}, // Number of points to test.
+ TestUtilities.createRandomNumberGenerator());
}
/**
* Computes {@link ConformalProjection#φ(double)}.
+ * The {@link #transform} field must have been set before this method is invoked.
*
* @param expOfSouthing the reciprocal of the value returned by {@link #expOfNorthing(double)}.
* @return the latitude in radians.
* @throws ProjectionException if the iteration does not converge.
+ *
+ * @see #expOfNorthing(double)
*/
private double φ(final double expOfSouthing) throws ProjectionException {
return ((ConformalProjection) transform).φ(expOfSouthing);
@@ -234,21 +226,19 @@ public final strictfp class ConformalProjectionTest extends TransformTestCase {
* @throws ProjectionException if an error occurred while projecting a point.
*/
@Test
- @DependsOnMethod("testExpOfNorthing")
+ @DependsOnMethod("verifyNorthingKnownValues")
public void test_φ() throws ProjectionException {
- transform = new NoOp(false); // Spherical case
- tolerance = TOLERANCE;
- doTest_φ();
- transform = new NoOp(true); // Ellipsoidal case
- tolerance = NormalizedProjection.ITERATION_TOLERANCE;
- doTest_φ();
+ test_φ(false); // Spherical case
+ test_φ(true); // Ellipsoidal case
}
/**
- * Implementation of {@link #test_φ()}.
- * The {@link #transform} field must have been set before this method is called.
+ * Tests {@link ConformalProjection#φ(double)} on known values and as the reverse of {@code expOfNorthing(φ)}.
*/
- private void doTest_φ() throws ProjectionException {
+ private void test_φ(final boolean ellipsoidal) throws ProjectionException {
+ transform = new NoOp(ellipsoidal);
+ tolerance = ellipsoidal ? NormalizedProjection.ITERATION_TOLERANCE : TOLERANCE;
+
assertEquals("φ(NaN) = NaN", NaN, φ(NaN), tolerance);
assertEquals("φ( ∞) = -90°", -PI/2, φ(POSITIVE_INFINITY), tolerance);
assertEquals("φ( ∞) = -90°", -PI/2, φ(MAX_VALUE), tolerance);
@@ -282,21 +272,18 @@ public final strictfp class ConformalProjectionTest extends TransformTestCase {
* <ol>
* <li>φ values computed by an iterative method.</li>
* <li>φ values computed by the series expansion given by EPSG guide.</li>
- * <li>φ values computed by modified form of series expansion, using trigonometric identities.</li>
- * <li>φ values computed by the actual {@link ConformalProjection} implementation.</li>
+ * <li>φ values computed by the actual {@link ConformalProjection} implementation,
+ * which uses a modified form of series expansion using trigonometric identities.</li>
* </ol>
*
- * {@link ConformalProjection#φ(double)} which uses a mix of 1 and 3 in the above list.
- * See {@link MercatorMethodComparison} for a discussion.
- *
* @throws ProjectionException if an error occurred during computation of φ.
*
* @see MercatorMethodComparison
*/
@Test
- public void compareWithSeriesExpansion() throws ProjectionException {
+ public void compareReference() throws ProjectionException {
final ConformalProjection projection = new NoOp(true);
- final MercatorMethodComparison comparator = new MercatorMethodComparison(projection.eccentricitySquared);
+ final MercatorMethodComparison comparator = new MercatorMethodComparison(projection);
final Random random = TestUtilities.createRandomNumberGenerator();
final int numSamples = 2000;
for (int i=0; i<numSamples; i++) {
@@ -304,19 +291,16 @@ public final strictfp class ConformalProjectionTest extends TransformTestCase {
final double t = 1 / comparator.expOfNorthing(φ);
final double byIterativeMethod = comparator.byIterativeMethod(t);
final double bySeriesExpansion = comparator.bySeriesExpansion(t);
- final double byTrigoIdentities = comparator.usingTrigonometricIdentities(t);
final double byImplementation = projection.φ(t);
assertEquals("Iterative method", φ, byIterativeMethod, 1E-11);
assertEquals("Series expansion", φ, bySeriesExpansion, 1E-11);
- assertEquals("Trigo. identities", φ, byTrigoIdentities, 1E-11);
assertEquals("Implementation", φ, byImplementation, 1E-11);
/*
* Verify that the formulas modified with trigonometric identities give the same results
* than the original formulas. The main purpose of this test is to detect mistake during
* the application of identities.
*/
- assertEquals(byTrigoIdentities, bySeriesExpansion, 1E-15); // Tolerance threshold close to 1 ULP of 2π.
- assertEquals(projection.φ(t), byTrigoIdentities, 1E-15);
+ assertEquals(byImplementation, bySeriesExpansion, 1E-15); // Tolerance threshold close to 1 ULP of 2π.
}
}
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
index dd31382..342fab2 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/MercatorMethodComparison.java
@@ -23,7 +23,6 @@ import java.io.UncheckedIOException;
import org.apache.sis.io.TableAppender;
import org.apache.sis.math.Statistics;
import org.apache.sis.math.StatisticsFormat;
-import org.apache.sis.util.ArraysExt;
import org.apache.sis.internal.referencing.Resources;
import org.apache.sis.internal.metadata.ReferencingServices;
@@ -41,9 +40,9 @@ import static org.apache.sis.internal.util.StandardDateFormat.NANOS_PER_SECOND;
* <li>an iterative process for solving equation (7-9) from Snyder, initially implemented by USGS.</li>
* </ul>
*
- * In our measurements, both the iterative process (USGS) and the series expansion (EPSG) have the
- * same accuracy when applied on the WGS84 ellipsoid. However the EPSG formula is 2 times faster.
- * On the other hand, accuracy of the EPSG formula decreases when we increase the eccentricity,
+ * In our measurements, both the iterative process (USGS) and the series expansion (EPSG) have the same accuracy
+ * when applied on the WGS84 ellipsoid. However the EPSG formula is 2 times faster on Java 8 (less on more recent
+ * Java versions). On the other hand, accuracy of the EPSG formula decreases when we increase the eccentricity,
* while the iterative process keeps its accuracy (at the cost of more iterations).
* For the Earth (eccentricity of about 0.082) the errors are less than 0.01 millimetres.
* But the errors become centimetric (for a hypothetical planet of the size of the Earth)
@@ -55,7 +54,7 @@ import static org.apache.sis.internal.util.StandardDateFormat.NANOS_PER_SECOND;
* win in about 50% of cases. But as we increase the eccentricity, the iterative method wins more often.</p>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 0.6
+ * @version 1.0
* @since 0.6
* @module
*/
@@ -77,29 +76,19 @@ public final class MercatorMethodComparison { // No 'strictfp' keyword here si
private final double c2χ, c4χ, c6χ, c8χ;
/**
- * Creates a new instance for the eccentricity of the WGS84 ellipsoid, which is approximately 0.08181919084262157.
- * Reminder: the eccentricity of a sphere is 0.
+ * The actual implementation to compare with.
*/
- public MercatorMethodComparison() {
- this(0.00669437999014133); // Squared eccentricity.
- }
+ private final ConformalProjection implementation;
/**
* Creates a new instance for the same eccentricity than the given projection.
*
* @param projection the projection from which to take the eccentricity.
*/
- public MercatorMethodComparison(final NormalizedProjection projection) {
- this(projection.eccentricitySquared);
- }
-
- /**
- * Creates a new instance for the given squared eccentricity.
- *
- * @param e2 the square of the eccentricity.
- */
- public MercatorMethodComparison(final double e2) {
- eccentricity = sqrt(e2);
+ MercatorMethodComparison(final ConformalProjection projection) {
+ implementation = projection;
+ eccentricity = projection.eccentricity;
+ final double e2 = projection.eccentricitySquared;
final double e4 = e2 * e2;
final double e6 = e2 * e4;
final double e8 = e4 * e4;
@@ -129,30 +118,6 @@ public final class MercatorMethodComparison { // No 'strictfp' keyword here si
}
/**
- * Same formula than {@link #bySeriesExpansion(double)}, but replacing some sine by trigonometric identities.
- * The identities used are:
- *
- * <ul>
- * <li>sin(2⋅x) = 2⋅sin(x)⋅cos(x)</li>
- * <li>sin(3⋅x) = (3 - 4⋅sin²(x))⋅sin(x)</li>
- * <li>sin(4⋅x) = (4 - 8⋅sin²(x))⋅sin(x)⋅cos(x)</li>
- * </ul>
- *
- * @param t the {@code expOfSouthing} parameter value.
- * @return the latitude (in radians) for the given parameter.
- */
- public double usingTrigonometricIdentities(final double t) {
- final double χ = PI/2 - 2*atan(t);
- final double sin2χ = sin(2*χ);
- final double sin_cos2χ = cos(2*χ) * sin2χ;
- final double sin_sin2χ = sin2χ * sin2χ;
- return c8χ * (4 - 8*sin_sin2χ)*sin_cos2χ
- + c6χ * (3 - 4*sin_sin2χ)*sin2χ
- + c4χ * 2*sin_cos2χ
- + c2χ * sin2χ + χ;
- }
-
- /**
* Computes φ using the iterative method used by USGS.
* This is the second part of the {@link ConformalProjection#φ(double)} method.
*
@@ -193,34 +158,28 @@ public final class MercatorMethodComparison { // No 'strictfp' keyword here si
* @throws ProjectionException if an error occurred during the calculation of φ.
*/
public void printAccuracyComparison(final int numSamples) throws ProjectionException {
- compare(null, numSamples, null);
+ compare(numSamples, null);
}
/**
* Implementation of {@link #printAccuracyComparison(int)} and {@link #printErrorForExcentricities(double,double)},
* optionally with a comparison with {@link ConformalProjection}.
*/
- private void compare(final ConformalProjection projection, final int numSamples, final TableAppender summarize)
- throws ProjectionException
- {
+ private void compare(final int numSamples, final TableAppender summarize) throws ProjectionException {
final Statistics iterativeMethodErrors = new Statistics("Iterative method error");
final Statistics seriesExpansionErrors = new Statistics("Series expansion error");
- final Statistics usingTrigoIdentErrors = new Statistics("Using trigonometric identities");
- final Statistics abstractLambertErrors = new Statistics("'ConformalProjection' error");
+ final Statistics implementationErrors = new Statistics("Implementation error");
final Random random = new Random();
for (int i=0; i<numSamples; i++) {
final double φ = random.nextDouble() * PI - PI/2;
final double t = 1 / expOfNorthing(φ);
final double byIterativeMethod = byIterativeMethod(t);
final double bySeriesExpansion = bySeriesExpansion(t);
- final double usingTrigoIdent = usingTrigonometricIdentities(t);
+ final double byImplementation = implementation.φ(t);
iterativeMethodErrors.accept(abs(φ - byIterativeMethod) / NormalizedProjection.ITERATION_TOLERANCE);
seriesExpansionErrors.accept(abs(φ - bySeriesExpansion) / NormalizedProjection.ITERATION_TOLERANCE);
- usingTrigoIdentErrors.accept(abs(φ - usingTrigoIdent) / NormalizedProjection.ITERATION_TOLERANCE);
- if (projection != null) {
- abstractLambertErrors.accept(abs(φ - projection.φ(t)) / NormalizedProjection.ITERATION_TOLERANCE);
- }
+ implementationErrors .accept(abs(φ - byImplementation) / NormalizedProjection.ITERATION_TOLERANCE);
}
/*
* At this point we finished to collect the statistics for the eccentricity of this particular
@@ -238,12 +197,8 @@ public final class MercatorMethodComparison { // No 'strictfp' keyword here si
Statistics[] stats = new Statistics[] {
iterativeMethodErrors,
seriesExpansionErrors,
- usingTrigoIdentErrors,
- abstractLambertErrors
+ implementationErrors,
};
- if (projection == null) {
- stats = ArraysExt.remove(stats, 2, 1);
- }
out.println("Comparison of different ways to compute φ for eccentricity " + eccentricity + '.');
out.println("Values are in units of " + NormalizedProjection.ITERATION_TOLERANCE + " radians (about "
+ round(toDegrees(NormalizedProjection.ITERATION_TOLERANCE) * 60 * ReferencingServices.NAUTICAL_MILE * 1000)
@@ -288,8 +243,9 @@ public final class MercatorMethodComparison { // No 'strictfp' keyword here si
crossThreshold = true;
table.appendHorizontalSeparator();
}
- final MercatorMethodComparison alt = new MercatorMethodComparison(eccentricity * eccentricity);
- alt.compare(null, 10000, table);
+ final double semiMinor = sqrt(1 - eccentricity * eccentricity);
+ final MercatorMethodComparison alt = new MercatorMethodComparison(new NoOp(1, semiMinor));
+ alt.compare(10000, table);
}
table.appendHorizontalSeparator();
try {
@@ -321,13 +277,13 @@ public final class MercatorMethodComparison { // No 'strictfp' keyword here si
}
final long t2 = System.nanoTime();
for (int i=0; i<t.length; i++) {
- s2 += usingTrigonometricIdentities(t[i]);
+ s2 += implementation.φ(t[i]);
}
final long t3 = System.nanoTime();
final float c = (t1 - t0) / 100f;
- out.println("Iterative method: " + ((t1 - t0) / (float) NANOS_PER_SECOND) + " seconds (" + round((t1 - t0) / c) + "%).");
- out.println("Series expansion: " + ((t2 - t1) / (float) NANOS_PER_SECOND) + " seconds (" + round((t2 - t1) / c) + "%).");
- out.println("Trigonometric identities: " + ((t3 - t2) / (float) NANOS_PER_SECOND) + " seconds (" + round((t3 - t2) / c) + "%).");
+ out.println("Iterative method: " + ((t1 - t0) / (float) NANOS_PER_SECOND) + " seconds (" + round((t1 - t0) / c) + "%).");
+ out.println("Series expansion: " + ((t2 - t1) / (float) NANOS_PER_SECOND) + " seconds (" + round((t2 - t1) / c) + "%).");
+ out.println("Implementation: " + ((t3 - t2) / (float) NANOS_PER_SECOND) + " seconds (" + round((t3 - t2) / c) + "%).");
out.println("Mean φ values: " + (s0 / t.length) + ", "
+ (s1 / t.length) + " and "
+ (s2 / t.length) + ".");
@@ -348,29 +304,26 @@ public final class MercatorMethodComparison { // No 'strictfp' keyword here si
out.println("Comparison of the errors for a sphere.");
out.println("The errors should be almost zero:");
out.println();
- ConformalProjection projection = new NoOp(false);
- MercatorMethodComparison c = new MercatorMethodComparison(projection);
- c.compare(projection, 10000, null);
+ MercatorMethodComparison c = new MercatorMethodComparison(new NoOp(false));
+ c.compare(10000, null);
out.println();
out.println("Comparison of the errors for the WGS84 eccentricity.");
out.println("The 'ConformalProjection' errors should be the same than the series expansion errors:");
out.println();
- projection = new NoOp(true);
- c = new MercatorMethodComparison(projection);
- c.compare(projection, 1000000, null);
+ c = new MercatorMethodComparison(new NoOp(true));
+ c.compare(1000000, null);
out.println();
out.println("Comparison of the errors for the eccentricity of an imaginary ellipsoid.");
out.println("The 'ConformalProjection' errors should be the close to the iterative method errors:");
out.println();
- projection = new NoOp(100, 95);
- c = new MercatorMethodComparison(projection);
- c.compare(projection, 1000000, null);
+ c = new MercatorMethodComparison(new NoOp(100, 95));
+ c.compare(1000000, null);
out.println();
out.println("Benchmarks");
- c = new MercatorMethodComparison();
+ c = new MercatorMethodComparison(new NoOp(true));
for (int i=0; i<4; i++) {
System.gc();
Thread.sleep(1000);
[sis] 03/03: Include MeridianArcBased fields in serialization for
the same simplicity and numerical stability reasons than
ConformalProjection and EqualAreaProjection. Inconvenient is that the
serialization is more at risk to break if the formulas implemented in those
case are improved. However MathTransform serialization should not be used
for long-term storage. It is okay for RMI purpose between JVMs running the
same version of SIS.
Posted by de...@apache.org.
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 5b271a0a011fd1559f5548ef972f3dc7e09feb04
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Mon Apr 29 14:39:10 2019 +0200
Include MeridianArcBased fields in serialization for the same simplicity and numerical stability reasons than ConformalProjection and EqualAreaProjection.
Inconvenient is that the serialization is more at risk to break if the formulas implemented in those case are improved. However MathTransform serialization
should not be used for long-term storage. It is okay for RMI purpose between JVMs running the same version of SIS.
---
.../operation/projection/MeridianArcBased.java | 109 +++++++++------------
1 file changed, 44 insertions(+), 65 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
index 81a8010..4f005ea 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/MeridianArcBased.java
@@ -16,8 +16,6 @@
*/
package org.apache.sis.referencing.operation.projection;
-import java.io.IOException;
-import java.io.ObjectInputStream;
import org.apache.sis.internal.util.DoubleDouble;
import static java.lang.Math.*;
@@ -40,72 +38,33 @@ abstract class MeridianArcBased extends NormalizedProjection {
/**
* For cross-version compatibility.
*/
- private static final long serialVersionUID = -1390453290591193455L;
+ private static final long serialVersionUID = 6105123637473385555L;
/**
- * Coefficients for the formula implemented by the {@link #distance(double, double, double)} method.
- * Values are computed by {@link #computeCoefficients(double)} in constructor or after deserialization.
- * The values depends on the form of equation implemented by {@code distance(…)}. We do not use the
- * form commonly found in publications since a few algebraic operations allow to replace the
- * sin(2φ), sin(4φ), sin(6φ) and sin(8φ) terms by only sin²(φ), which is faster to calculate.
- * See {@link #computeCoefficients(double)} comments in method body for more information.
+ * Coefficients for the forward transform implemented by the {@link #distance(double, double, double)} method.
+ * Values are computed by the constructor and depend on the form of equation implemented by {@code distance(…)}.
+ * We do not use the formulas in the form published by EPSG or Snyder since a few algebraic operations allow to
+ * replace the sin(2φ), sin(4φ), sin(6φ) and sin(8φ) terms by sin²(φ), which is faster to calculate.
+ * See comments in constructor body for more information.
*/
- private transient double c0, c1, c2, c3, c4;
+ private final double cf0, cf1, cf2, cf3, cf4;
/**
- * Coefficients used in inverse transform. Same comment than for {@link #c0}… applies.
+ * Coefficients used in inverse transform. Same comment than for {@link #cf0}… applies.
*/
- private transient double ci1, ci2, ci3, ci4;
+ private final double ci1, ci2, ci3, ci4;
/**
* Denominator of <cite>rectifying latitude</cite> equation. The rectifying latitude is computed by
* µ = M/(1 – ℯ²/4 – 3ℯ⁴/64 – 5ℯ⁶/256 – …) (Snyder 7-19 with a=1).
*/
- private transient double rld;
+ private final double rµ;
/**
* Creates a new normalized projection from the parameters computed by the given initializer.
*/
MeridianArcBased(final Initializer initializer) {
super(initializer);
- DoubleDouble e1 = initializer.axisLengthRatio();
- e1.ratio_1m_1p();
- computeCoefficients(e1.doubleValue());
- }
-
- /**
- * Creates a new projection initialized to the same parameters than the given one.
- */
- MeridianArcBased(final MeridianArcBased other) {
- super(other);
- c0 = other.c0;
- c1 = other.c1;
- c2 = other.c2;
- c3 = other.c3;
- c4 = other.c4;
- ci1 = other.ci1;
- ci2 = other.ci2;
- ci3 = other.ci3;
- ci4 = other.ci4;
- rld = other.rld;
- }
-
- /**
- * Restores transient fields after deserialization.
- */
- private void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException {
- in.defaultReadObject();
- double ba = sqrt(1 - eccentricitySquared);
- computeCoefficients((1 - ba) / (1 + ba));
- }
-
- /**
- * Computes the coefficients in the series expansions from the {@link #eccentricitySquared} value.
- * This method shall be invoked after {@code MeridianArcBased} construction or deserialization.
- *
- * @param ei value of [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)] (Snyder 3-24).
- */
- private void computeCoefficients(final double ei) {
final double e2 = eccentricitySquared;
final double e4 = e2*e2;
final double e6 = e2*e4;
@@ -155,11 +114,11 @@ abstract class MeridianArcBased extends NormalizedProjection {
* of 2, which result in exact representations in IEEE 754 double type. Derivation from Kawase formula can be
* viewed at: https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods
*/
- c0 = -441./0x10000 * e10 + -175./0x4000 * e8 + -5./0x100 * e6 + -3./0x40 * e4 + -1./4 * e2 + 1;
- c1 = 1575./0x20000 * e10 + 175./0x4000 * e8 + 5./0x100 * e6 + 3./0x40 * e4 + -3./4 * e2;
- c2 = -2625./0x8000 * e10 + 175./0x6000 * e8 + 5120./0x60000 * e6 + -15./0x20 * e4;
- c3 = 735./0x800 * e10 + 2240./0x60000 * e8 + -35./0x60 * e6;
- c4 = -2205./0x1000 * e10 + -315./0x400 * e8;
+ cf0 = -441./0x10000 * e10 + -175./0x4000 * e8 + -5./0x100 * e6 + -3./0x40 * e4 + -1./4 * e2 + 1;
+ cf1 = 1575./0x20000 * e10 + 175./0x4000 * e8 + 5./0x100 * e6 + 3./0x40 * e4 + -3./4 * e2;
+ cf2 = -2625./0x8000 * e10 + 175./0x6000 * e8 + 5120./0x60000 * e6 + -15./0x20 * e4;
+ cf3 = 735./0x800 * e10 + 2240./0x60000 * e8 + -35./0x60 * e6;
+ cf4 = -2205./0x1000 * e10 + -315./0x400 * e8;
// c6 = 693./0x20000 * e10 omitted for now (not yet used).
/*
* Coefficients for inverse transform derived from Snyder 3-26 and EPSG guidance notes:
@@ -176,7 +135,10 @@ abstract class MeridianArcBased extends NormalizedProjection {
* µ is the rectifying latitude.
* Derivation of coefficients used by this class are provided in the above-cited spreadsheet.
*/
- rld = (1 - 1./4*e2 - 3./64*e4 - 5./256*e6); // Part of Snyder 7-19 for computing rectifying latitude.
+ rµ = (1 - 1./4*e2 - 3./64*e4 - 5./256*e6); // Part of Snyder 7-19 for computing rectifying latitude.
+ DoubleDouble e1 = initializer.axisLengthRatio();
+ e1.ratio_1m_1p();
+ final double ei = e1.doubleValue(); // Equivalent to [1 - √(1 - ℯ²)] / [1 + √(1 - ℯ²)] (Snyder 3-24).
final double ei2 = ei*ei;
final double ei3 = ei*ei2;
final double ei4 = ei2*ei2;
@@ -187,6 +149,23 @@ abstract class MeridianArcBased extends NormalizedProjection {
}
/**
+ * Creates a new projection initialized to the same parameters than the given one.
+ */
+ MeridianArcBased(final MeridianArcBased other) {
+ super(other);
+ cf0 = other.cf0;
+ cf1 = other.cf1;
+ cf2 = other.cf2;
+ cf3 = other.cf3;
+ cf4 = other.cf4;
+ ci1 = other.ci1;
+ ci2 = other.ci2;
+ ci3 = other.ci3;
+ ci4 = other.ci4;
+ rµ = other.rµ;
+ }
+
+ /**
* Computes the distance (M) along meridian arc from equator to a given latitude φ.
* Special cases:
*
@@ -196,7 +175,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
* <li>If φ=-π/2, then this method returns a value slightly greater than -π/2, depending on the eccentricity.</li>
* </ul>
*
- * This value is related to <cite>rectifying latitude</cite> by µ = (π/2)⋅(M/{@linkplain #c0}) (derived from Snyder 3-20).
+ * This value is related to <cite>rectifying latitude</cite> by µ = (π/2)⋅(M/{@linkplain #cf0}) (derived from Snyder 3-20).
*
* @param φ latitude for which to compute the distance, in radians.
* @param sinφ value of sin(φ).
@@ -205,7 +184,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
*/
final double distance(final double φ, final double sinφ, final double cosφ) {
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.
+ return cf0*φ + cosφ*sinφ*(cf1 + sinφ2*(cf2 + sinφ2*(cf3 + sinφ2*cf4))); // TODO: use Math.fma with JDK9.
}
/**
@@ -214,11 +193,11 @@ abstract class MeridianArcBased extends NormalizedProjection {
* @return the derivative at the specified latitude.
*/
final double dM_dφ(final double sinφ2) {
- return ((((7 - 8*sinφ2)*c4 - 6*c3) * sinφ2
- + 5*c3 - 4*c2) * sinφ2
- + 3*c2 - 2*c1) * sinφ2
- + c1
- + c0;
+ return ((((7 - 8*sinφ2)*cf4 - 6*cf3) * sinφ2
+ + 5*cf3 - 4*cf2) * sinφ2
+ + 3*cf2 - 2*cf1) * sinφ2
+ + cf1
+ + cf0;
}
/**
@@ -228,7 +207,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
* @return the latitude of given meridian distance, in radians.
*/
final double latitude(final double distance) {
- double φ = distance / rld; // Rectifying latitude µ used as first approximation.
+ double φ = distance / rµ; // Rectifying latitude µ used as first approximation.
double sinφ = sin(φ);
double sinφ2 = sinφ * sinφ;
φ += cos(φ)*sinφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4))); // Snyder 3-26.
[sis] 02/03: Apply to EqualAreaProjection the same Snyder 3-35 formulas (trigonometric identities) than we did for MeridianArcBased. Include coefficients and other fields in serialization for the same reason than ConformalProjection (simplicity, stability). During this work, we identified a bug in the series expansion of φ(β) where the last term was multiplied by e⁸ instead of e⁶. Fixing that bug bring the series expansion to a precision where we no longer need iterations at least with WGS84 ellipsoid.
Posted by de...@apache.org.
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 6760894f16bca4bfcfbd9e12c8d0daf914b9a089
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Mon Apr 29 14:10:24 2019 +0200
Apply to EqualAreaProjection the same Snyder 3-35 formulas (trigonometric identities) than we did for MeridianArcBased.
Include coefficients and other fields in serialization for the same reason than ConformalProjection (simplicity, stability).
During this work, we identified a bug in the series expansion of φ(β) where the last term was multiplied by e⁸ instead of e⁶.
Fixing that bug bring the series expansion to a precision where we no longer need iterations at least with WGS84 ellipsoid.
---
.../operation/projection/AlbersEqualArea.java | 1 -
.../operation/projection/ConformalProjection.java | 14 +-
.../operation/projection/CylindricalEqualArea.java | 1 -
.../operation/projection/EqualAreaProjection.java | 202 ++++++++++-----------
.../projection/EqualAreaProjectionTest.java | 127 +++++++++++++
.../sis/test/suite/ReferencingTestSuite.java | 1 +
6 files changed, 227 insertions(+), 119 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
index c2a137a..c7de81a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AlbersEqualArea.java
@@ -161,7 +161,6 @@ public class AlbersEqualArea extends EqualAreaProjection {
denormalize.convertBefore(0, rn, null); rn.negate();
denormalize.convertBefore(1, rn, ρ0); rn.inverseDivide(-1);
normalize.convertAfter(0, rn, null);
- super.computeCoefficients();
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
index 35f0151..c2a9045 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/ConformalProjection.java
@@ -23,10 +23,9 @@ import static java.lang.Math.*;
/**
* Base class of {@link LambertConicConformal}, {@link Mercator} and {@link PolarStereographic} projections.
- * All those projections have in common the property of being <cite>conformal</cite>, i.e. they preserve
- * angles locally. However we do not put this base class in public API because not all conformal projections
- * will extend this base class. For example {@link TransverseMercator}, despite being conformal, uses very
- * different formulas.
+ * Those projections have in common the property of being <cite>conformal</cite>, i.e. they preserve angles locally.
+ * However we do not put this base class in public API because not all conformal projections extend this base class.
+ * For example the {@link TransverseMercator} projection, despite being conformal, uses very different formulas.
*
* <p>Note that no projection can be both conformal and equal-area. This restriction is implemented in class
* hierarchy with {@link ConformalProjection} and {@link EqualAreaProjection} being two distinct classes.</p>
@@ -85,8 +84,9 @@ abstract class ConformalProjection extends NormalizedProjection {
private final boolean useIterations;
/**
- * Coefficients in the series expansion of the inverse projection, depending only on {@linkplain #eccentricity
- * eccentricity} value. The series expansion is published under the following form, where χ is the conformal latitude:
+ * Coefficients of the first terms in the series expansion of the inverse projection.
+ * Values of those coefficients depend only on {@linkplain #eccentricity eccentricity} value.
+ * The series expansion is published under the following form, where χ is the <cite>conformal latitude</cite>:
*
* <blockquote>c₂⋅sin(2χ) + c₄⋅sin(4χ) + c₆⋅sin(6χ) + c₈⋅sin(8χ)</blockquote>
*
@@ -200,7 +200,7 @@ abstract class ConformalProjection extends NormalizedProjection {
* Note that the φ value computed by the line below is called χ in EPSG guide.
* We name it φ in our code because we will modify that value in-place in order to get φ.
*/
- double φ = (PI/2) - 2*atan(expOfSouthing); // at this point == χ
+ double φ = (PI/2) - 2*atan(expOfSouthing); // at this point == χ (conformal latitude).
/*
* Add a correction for the flattened shape of the Earth. The correction can be represented by an
* infinite series. Here, we apply only the first 4 terms. Those terms are given by §1.3.3 in the
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
index 6cc9fb6..7fda325 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/CylindricalEqualArea.java
@@ -168,7 +168,6 @@ public class CylindricalEqualArea extends EqualAreaProjection {
ik.divide(k0);
denormalize.convertAfter(0, k0, null);
denormalize.convertAfter(1, ik, null);
- super.computeCoefficients();
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
index 26c6579..3999e54 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/EqualAreaProjection.java
@@ -16,8 +16,6 @@
*/
package org.apache.sis.referencing.operation.projection;
-import java.io.IOException;
-import java.io.ObjectInputStream;
import org.apache.sis.internal.referencing.Resources;
import static java.lang.Math.*;
@@ -25,10 +23,16 @@ import static org.apache.sis.math.MathFunctions.atanh;
/**
- * Provides formulas common to Equal Area projections.
+ * Base class of {@link AlbersEqualArea} and {@link CylindricalEqualArea} projections.
+ * Those projections have in common the property of being <cite>equal-area</cite>.
+ * However we do not put this base class in public API because not all equal-area projections extend this base class.
+ * For example the {@link Sinusoidal} projection, despite being equal-area, uses different formulas.
+ *
+ * <p>Note that no projection can be both conformal and equal-area. This restriction is implemented in class
+ * hierarchy with {@link ConformalProjection} and {@link EqualAreaProjection} being two distinct classes.</p>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.0
* @since 0.8
* @module
*/
@@ -36,56 +40,48 @@ abstract class EqualAreaProjection extends NormalizedProjection {
/**
* For cross-version compatibility.
*/
- private static final long serialVersionUID = -6175270149094989517L;
+ private static final long serialVersionUID = 2221537810038553082L;
/**
- * {@code false} for using the original formulas as published by Snyder, or {@code true} for using formulas
- * modified using trigonometric identities. The use of trigonometric identities is for reducing the amount
- * of calls to the {@link Math#sin(double)} and similar methods. Some identities used are:
- *
- * <ul>
- * <li>sin(2β) = 2⋅sinβ⋅cosβ</li>
- * <li>sin(4β) = (2 - 4⋅sin²β)⋅sin(2β)</li>
- * <li>sin(8β) = 4⋅sin(2β)⋅(cos²β - sin²β)⋅(8⋅cos⁴β - 8⋅cos²β + 1)</li>
- * </ul>
- *
- * Note that since this boolean is static final, the compiler should exclude the code in the branch that is never
- * executed (no need to comment-out that code).
+ * The threshold value of {@link #eccentricity} at which we consider the accuracy of the
+ * series expansion insufficient. This threshold is determined empirically with the help
+ * of the {@code EqualAreaProjectionTest.searchThreshold()} method in the test directory.
*
- *
- * <p><b>BENCHMARK AND ANALYSIS:</b>
- * as of July 2016, benchmarking shows small benefit in using trigonometric identities for {@code EqualAreaProjection}
- * (contrarily to {@link ConformalProjection} where we measured a greater benefit). This may be because in this class,
- * the series expansion is unconditionally followed by iterative method in order to reach the centimetric precision.
- * We observe that the original series expansion allows convergence in only one iteration, while the formulas using
- * trigonometric identifies often requires two iterations.</p>
- *
- * @todo Redo the benchmark with JDK9, since {@code sin} seems much faster.
+ * <p>Note: WGS84 eccentricity is about 0.082.</p>
*/
- private static final boolean ALLOW_TRIGONOMETRIC_IDENTITIES = true;
+ private static final double ECCENTRICITY_THRESHOLD = 0.1;
/**
- * Coefficients in the series expansion of the inverse projection,
- * depending only on {@linkplain #eccentricity eccentricity} value.
- * The series expansion is of the following form:
+ * Coefficients of the first terms in the series expansion of the inverse projection.
+ * Values of those coefficients depend only on {@linkplain #eccentricity eccentricity} value.
+ * The series expansion is published under the following form, where β is the <cite>authalic latitude</cite>:
*
- * <blockquote>φ = ci₂⋅sin(2β) + ci₄⋅sin(4β) + ci₈⋅sin(8β)</blockquote>
+ * <blockquote>φ = c₂⋅sin(2β) + c₄⋅sin(4β) + c₈⋅sin(6β)</blockquote>
*
- * This {@code EqualAreaProjection} class uses those coefficients in {@link #φ(double)}.
+ * However we rewrite above series expansion for taking advantage of trigonometric identities.
+ * The equation become (with different <var>c</var> values):
*
- * <p><strong>Consider those fields as final!</strong> They are not final only for sub-class
- * constructors convenience and for the purpose of {@link #readObject(ObjectInputStream)}.</p>
+ * <blockquote>sin(2β)⋅(c₂ + cos(2β)⋅(c₄ + cos(2β)⋅c₆))</blockquote>
*
- * @see #computeCoefficients()
+ * <div class="note"><b>Serialization note:</b>
+ * we do not strictly need to serialize those fields since they could be computed after deserialization.
+ * Bu we serialize them anyway in order to simplify a little bit this class (it allows us to keep those
+ * fields final) and because values computed after deserialization could be slightly different than the
+ * ones computed after construction if a future version of the constructor uses the double-double values
+ * provided by {@link Initializer}.</div>
*/
- private transient double ci2, ci4, ci8;
+ private final double c2β, c4β, c6β;
/**
* Value of {@link #qm(double)} function (part of Snyder equation (3-12)) at pole (sinφ = 1).
- *
- * @see #computeCoefficients()
*/
- private transient double qmPolar;
+ private final double qmPolar;
+
+ /**
+ * {@code true} if the {@link #eccentricity} value is greater than or equals to the eccentricity threshold,
+ * in which case the {@link #φ(double)} method will need to use an iterative method.
+ */
+ private final boolean useIterations;
/**
* Creates a new normalized projection from the parameters computed by the given initializer.
@@ -94,31 +90,37 @@ abstract class EqualAreaProjection extends NormalizedProjection {
*/
EqualAreaProjection(final Initializer initializer) {
super(initializer);
- }
-
- /**
- * Computes the coefficients in the series expansions from the {@link #eccentricitySquared} value.
- * This method shall be invoked after {@code EqualAreaProjection} construction or deserialization.
- */
- void computeCoefficients() {
final double e2 = eccentricitySquared;
- final double e4 = e2 * e2;
- final double e6 = e2 * e4;
- ci2 = 517/5040. * e6 + 31/180. * e4 + 1/3. * e2;
- ci4 = 251/3780. * e6 + 23/360. * e4;
- ci8 = 761/45360. * e6;
+ final double e4 = e2 * e2;
+ final double e6 = e2 * e4;
/*
- * When rewriting equations using trigonometric identities, some constants appear.
- * For example sin(2β) = 2⋅sinβ⋅cosβ, so we can factor out the 2 constant into the
- * into the corresponding 'c' field.
+ * Published equation is of the following form:
+ *
+ * φ = β + A⋅sin(2β) + B⋅sin(4β) + C⋅sin(6β) (part of Snyder 3-34)
+ *
+ * We rewrite as:
+ *
+ * φ = β + sin(2β)⋅(A′ + cos(2β)⋅(B′ + cos(2β)⋅C′)) (part of Snyder 3-35)
+ * A′ = A - C
+ * B′ = 2B - 4D
+ * C′ = 4C
+ *
+ * Published coefficients are:
+ *
+ * A = 517/5040 ⋅ℯ⁶ + 31/180 ⋅ℯ⁴ + ℯ²/3
+ * B = 251/3780 ⋅ℯ⁶ + 23/360 ⋅ℯ⁴
+ * C = 761/45360 ⋅ℯ⁶
+ *
+ * We replace A B C D by A′ B′ C′ D′ using Snyder 3-35, then replace sin2β by 2sinβ⋅cosβ
+ * and cos2β by 1 – 2sin²β. The result are coefficients below. For each line, we add the
+ * smallest values first in order to reduce rounding errors. The algebraic changes are
+ * in https://svn.apache.org/repos/asf/sis/analysis/Map%20projection%20formulas.ods.
*/
- if (ALLOW_TRIGONOMETRIC_IDENTITIES) {
- // Multiplication by powers of 2 does not bring any additional rounding error.
- ci2 *= 2;
- ci4 *= 8;
- ci8 *= 64;
- }
+ c2β = 4./7 * e6 + 3./5 * e4 + 2./3 * e2;
+ c4β = -3028./2835 * e6 + -23./45 * e4;
+ c6β = 1522./2835 * e6;
qmPolar = qm(1);
+ useIterations = (eccentricity >= ECCENTRICITY_THRESHOLD);
}
/**
@@ -129,10 +131,11 @@ abstract class EqualAreaProjection extends NormalizedProjection {
*/
EqualAreaProjection(final EqualAreaProjection other) {
super(other);
- ci2 = other.ci2;
- ci4 = other.ci4;
- ci8 = other.ci8;
+ c2β = other.c2β;
+ c4β = other.c4β;
+ c6β = other.c6β;
qmPolar = other.qmPolar;
+ useIterations = other.useIterations;
}
/**
@@ -204,49 +207,36 @@ abstract class EqualAreaProjection extends NormalizedProjection {
*/
final double φ(final double y) throws ProjectionException {
final double sinβ = y / qmPolar;
+ final double sinβ2 = sinβ * sinβ;
final double β = asin(sinβ);
- double φ;
- if (!ALLOW_TRIGONOMETRIC_IDENTITIES) {
- φ = ci8 * sin(8*β)
- + ci4 * sin(4*β)
- + ci2 * sin(2*β)
- + β; // Snyder 3-18
- } else {
- /*
- * Same formula than above, but rewriten using trigonometric identities in order to avoid
- * multiple calls to sin(double) method. The cost is only one sqrt(double) method call.
- */
- final double sin2_β = sinβ*sinβ; // = sin²β
- final double cos2_β = 1 - sin2_β; // = cos²β
- final double t2β = sinβ * sqrt(cos2_β); // = sin(2β) / 2
- final double t4β = 0.5 - sin2_β; // = sin(4β) / ( 4⋅sin(2β))
- final double t8β = (cos2_β - sin2_β)*(cos2_β*cos2_β - cos2_β + 1./8); // = sin(8β) / (32⋅sin(2β))
-
- φ = (ci8*t8β + ci4*t4β + ci2) * t2β + β;
- }
/*
- * At this point φ is close to the desired value, but may have an error of a few centimetres.
- * Use the iterative method for reaching the last part of missing accuracy. Usually this loop
- * will perform exactly one iteration, no more, because φ is already quite close to the result.
- *
- * Mathematical note: Snyder 3-16 gives q/(1-ℯ²) instead of y in the calculation of Δφ below.
- * For Cylindrical Equal Area projection, Snyder 10-17 gives q = (qPolar⋅sinβ), which simplifies
- * as y.
- *
- * For Albers Equal Area projection, Snyder 14-19 gives q = (C - ρ²n²/a²)/n, which we rewrite
- * as q = (C - ρ²)/n (see comment in AlbersEqualArea.inverseTransform(…) for the mathematic).
- * The y value given to this method is y = (C - ρ²) / (n⋅(1-ℯ²)) = q/(1-ℯ²), the desired value.
+ * Snyder 3-18, but rewriten using trigonometric identities in order to avoid
+ * multiple calls to sin(double) method.
*/
- for (int i=0; i<MAXIMUM_ITERATIONS; i++) {
- final double sinφ = sin(φ);
- final double cosφ = cos(φ);
- final double ℯsinφ = eccentricity * sinφ;
- final double ome = 1 - ℯsinφ*ℯsinφ;
- final double Δφ = ome*ome/(2*cosφ) * (y - sinφ/ome - atanh(ℯsinφ)/eccentricity);
- φ += Δφ;
- if (abs(Δφ) <= ITERATION_TOLERANCE) {
- return φ;
+ double φ = β + cos(β)*sinβ*(c2β + sinβ2*(c4β + sinβ2*c6β));
+ if (useIterations) {
+ /*
+ * Mathematical note: Snyder 3-16 gives q/(1-ℯ²) instead of y in the calculation of Δφ below.
+ * For Cylindrical Equal Area projection, Snyder 10-17 gives q = (qPolar⋅sinβ), which simplifies
+ * as y.
+ *
+ * For Albers Equal Area projection, Snyder 14-19 gives q = (C - ρ²n²/a²)/n, which we rewrite
+ * as q = (C - ρ²)/n (see comment in AlbersEqualArea.inverseTransform(…) for the mathematic).
+ * The y value given to this method is y = (C - ρ²) / (n⋅(1-ℯ²)) = q/(1-ℯ²), the desired value.
+ */
+ for (int i=0; i<MAXIMUM_ITERATIONS; i++) {
+ final double sinφ = sin(φ);
+ final double cosφ = cos(φ);
+ final double ℯsinφ = eccentricity * sinφ;
+ final double ome = 1 - ℯsinφ*ℯsinφ;
+ final double Δφ = ome*ome/(2*cosφ) * (y - sinφ/ome - atanh(ℯsinφ)/eccentricity);
+ φ += Δφ;
+ if (abs(Δφ) <= ITERATION_TOLERANCE) {
+ return φ;
+ }
}
+ } else if (!Double.isNaN(φ)) {
+ return φ;
}
/*
* In the Albers Equal Area discussion, Snyder said that above algorithm does not converge if
@@ -275,12 +265,4 @@ abstract class EqualAreaProjection extends NormalizedProjection {
// Value should have converged but did not.
throw new ProjectionException(Resources.format(Resources.Keys.NoConvergence));
}
-
- /**
- * Restores transient fields after deserialization.
- */
- private void readObject(final ObjectInputStream in) throws IOException, ClassNotFoundException {
- in.defaultReadObject();
- computeCoefficients();
- }
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EqualAreaProjectionTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EqualAreaProjectionTest.java
new file mode 100644
index 0000000..40953a4
--- /dev/null
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/EqualAreaProjectionTest.java
@@ -0,0 +1,127 @@
+/*
+ * Licensed to the Apache Software Foundation (ASF) under one or more
+ * contributor license agreements. See the NOTICE file distributed with
+ * this work for additional information regarding copyright ownership.
+ * The ASF licenses this file to You under the Apache License, Version 2.0
+ * (the "License"); you may not use this file except in compliance with
+ * the License. You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+package org.apache.sis.referencing.operation.projection;
+
+import java.util.Random;
+import org.apache.sis.parameter.Parameters;
+import org.apache.sis.internal.referencing.provider.LambertCylindricalEqualArea;
+import org.apache.sis.referencing.operation.DefaultOperationMethod;
+import org.apache.sis.internal.util.Constants;
+import org.apache.sis.test.TestUtilities;
+import org.apache.sis.test.DependsOn;
+import org.junit.Test;
+
+import static java.lang.StrictMath.*;
+import static org.apache.sis.math.MathFunctions.atanh;
+import static org.junit.Assert.assertEquals;
+
+
+/**
+ * Tests {@link EqualAreaProjection}.
+ *
+ * @author Martin Desruisseaux (Geomatys)
+ * @version 1.0
+ * @since 1.0
+ * @module
+ */
+@DependsOn(NormalizedProjectionTest.class)
+public final strictfp class EqualAreaProjectionTest extends MapProjectionTestCase {
+ /**
+ * Creates the projection to be tested.
+ *
+ * @param ellipsoidal {@code false} for a sphere, or {@code true} for WGS84 ellipsoid.
+ * @return a test instance of the projection.
+ */
+ private EqualAreaProjection create(final boolean ellipsoidal) {
+ final DefaultOperationMethod provider = new LambertCylindricalEqualArea();
+ final CylindricalEqualArea projection = new CylindricalEqualArea(provider, parameters(provider, ellipsoidal));
+ tolerance = NormalizedProjection.ANGULAR_TOLERANCE; // = linear tolerance on a sphere of radius 1.
+ return projection;
+ }
+
+ /**
+ * Computes φ using equation given in EPSG guidance notes, which is also from Snyder book.
+ * We use this equation as a reference for testing validity of other forms.
+ *
+ * @param y in the cylindrical case, this is northing on the normalized ellipsoid.
+ * @return the latitude in radians.
+ */
+ private static double reference(final EqualAreaProjection projection, final double y) {
+ final double e = projection.eccentricity;
+ final double e2 = projection.eccentricitySquared;
+ final double e4 = e2 * e2;
+ final double e6 = e2 * e4;
+ final double c2β = 517./5040 * e6 + 31./180 * e4 + 1./3 * e2;
+ final double c4β = 251./3780 * e6 + 23./360 * e4;
+ final double c6β = 761./45360 * e6;
+ final double qmp = (1/(1 - e*e) + atanh(e)/e);
+ final double sinβ = y / qmp;
+ final double β = asin(sinβ);
+ return c6β * sin(6*β)
+ + c4β * sin(4*β)
+ + c2β * sin(2*β)
+ + β; // Snyder 3-18
+ }
+
+ /**
+ * Compares {@link EqualAreaProjection#φ(double)} with formula taken as references.
+ *
+ * @throws ProjectionException if the function does not converge.
+ */
+ @Test
+ public void compareWithReference() throws ProjectionException {
+ final EqualAreaProjection projection = create(true);
+ final Random random = TestUtilities.createRandomNumberGenerator();
+ for (int i=0; i<100; i++) {
+ final double y = random.nextDouble() * 3 - 1.5;
+ final double reference = reference(projection, y);
+ final double actual = projection.φ(y);
+ assertEquals(reference, actual, NormalizedProjection.ITERATION_TOLERANCE);
+ }
+ }
+
+ /**
+ * Searches a value for {@link EqualAreaProjection#ECCENTRICITY_THRESHOLD}.
+ * This method is not part of test suite. Steps to enable:
+ *
+ * <ol>
+ * <li>In {@link EqualAreaProjection#φ(double)} method, for {@code useIterations} to {@code false}.</li>
+ * <li>Add a {@link Test} annotation on this method.
+ * </ol>
+ *
+ * @throws ProjectionException if {@link EqualAreaProjection#φ(double)} did not converge.
+ */
+ @SuppressWarnings("UseOfSystemOutOrSystemErr")
+ public void searchThreshold() throws ProjectionException {
+ tolerance = NormalizedProjection.ANGULAR_TOLERANCE;
+ final DefaultOperationMethod provider = new LambertCylindricalEqualArea();
+ final Parameters parameters = parameters(provider, true);
+ for (double e = 0.05; e <= 0.2; e += 0.001) {
+ final double a = parameters.parameter(Constants.SEMI_MAJOR).doubleValue();
+ parameters.parameter(Constants.SEMI_MINOR).setValue(a * sqrt(1 - e*e));
+ final CylindricalEqualArea projection = new CylindricalEqualArea(provider, parameters);
+ for (double y = -1.25; y <= 1.25; y += 0.01) {
+ final double reference = reference(projection, y);
+ final double actual = projection.φ(y);
+ if (abs(actual - reference) > NormalizedProjection.ANGULAR_TOLERANCE) {
+ System.out.println("Error exceeds tolerance threshold at eccentricity " + e);
+ return;
+ }
+ }
+ }
+ }
+}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
index 43e533b..85d9c85 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java
@@ -173,6 +173,7 @@ import org.junit.BeforeClass;
org.apache.sis.referencing.operation.projection.PolarStereographicTest.class,
org.apache.sis.referencing.operation.projection.ObliqueStereographicTest.class,
org.apache.sis.referencing.operation.projection.ObliqueMercatorTest.class,
+ org.apache.sis.referencing.operation.projection.EqualAreaProjectionTest.class,
org.apache.sis.referencing.operation.projection.CylindricalEqualAreaTest.class,
org.apache.sis.referencing.operation.projection.AlbersEqualAreaTest.class,
org.apache.sis.referencing.operation.projection.MeridianArcTest.class,