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 468894b82c Use of FMA in series expansion of map projections.
Documentation fixes.
468894b82c is described below
commit 468894b82c294a557433fde891097c958d6b2c5a
Author: Martin Desruisseaux <martin.desruisse...@geomatys.com>
AuthorDate: Sat Dec 31 10:48:31 2022 +0100
Use of FMA in series expansion of map projections.
Documentation fixes.
---
.../operation/projection/AuthalicConversion.java | 18 ++++++++++-----
.../operation/projection/ConformalProjection.java | 27 ++++++++++++----------
.../projection/LambertAzimuthalEqualArea.java | 12 +---------
.../operation/projection/MeridianArcBased.java | 12 +++++-----
.../operation/transform/LinearTransform1D.java | 15 +++++++-----
.../operation/projection/SinusoidalTest.java | 2 +-
6 files changed, 44 insertions(+), 42 deletions(-)
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
index 8c681ea562..902268897b 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/AuthalicConversion.java
@@ -16,6 +16,7 @@
*/
package org.apache.sis.referencing.operation.projection;
+import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.referencing.Resources;
import static java.lang.Math.*;
@@ -37,11 +38,11 @@ import static org.apache.sis.math.MathFunctions.atanh;
* However, this {@code AuthalicConversion} is <strong>not</strong>
necessarily an authalic projection.
* Subclasses may want to use the latitude conversion formulas for other
purposes.
*
- * <h3>References</h3>
+ * <h2>References</h2>
* <p>Lee, L. P. "The Nomenclature and Classification of Map Projections."
Empire Survey Review 7, 190-200, 1944.</p>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.2
+ * @version 1.4
* @since 0.8
*/
abstract class AuthalicConversion extends NormalizedProjection {
@@ -141,9 +142,9 @@ abstract class AuthalicConversion extends
NormalizedProjection {
* 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.
*/
- c2β = 4./7 * e6 + 3./5 * e4 + 2./3 * e2;
- c4β = -3028./2835 * e6 + -23./45 * e4;
- c6β = 1522./2835 * e6;
+ c2β = fma(e2, 2./3, fma(e4, 3./5, e6*( 4./7)));
+ c4β = fma(e4, -23./45, e6*(-3028./2835));
+ c6β = e6*( 1522./2835);
qmPolar = qm(1);
useIterations = (eccentricity >= ECCENTRICITY_THRESHOLD);
}
@@ -252,7 +253,12 @@ abstract class AuthalicConversion extends
NormalizedProjection {
* Snyder 3-18, but rewriten using trigonometric identities in order
to avoid
* multiple calls to sin(double) method.
*/
- double φ = β + cos(β)*sinβ*(c2β + sinβ2*(c4β + sinβ2*c6β));
+ double φ;
+ if (Formulas.USE_FMA) {
+ φ = fma(fma(fma(sinβ2, c6β, c4β), sinβ2, c2β), cos(β)*sinβ, β);
+ } else {
+ φ = β + 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.
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 f23d690b2a..2e02096618 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
@@ -16,6 +16,7 @@
*/
package org.apache.sis.referencing.operation.projection;
+import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.referencing.Resources;
import static java.lang.Math.*;
@@ -54,7 +55,7 @@ import static java.lang.Math.*;
* parallel is a pole, the polar Stereographic results.” (Snyder, page
105)</div>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.4
* @since 0.6
*/
abstract class ConformalProjection extends NormalizedProjection {
@@ -94,12 +95,10 @@ abstract class ConformalProjection extends
NormalizedProjection {
*
* <blockquote>sin(2χ)⋅(c₂ + cos(2χ)⋅(c₄ + cos(2χ)⋅(c₆ +
cos(2χ)⋅c₈)))</blockquote>
*
- * <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>
+ * <h4>Serialization note</h4>
+ * we do not strictly need to serialize those fields because they could be
computed after deserialization.
+ * But we serialize them anyway in order to simplify a little bit this
class (it allows us to keep those
+ * fields final) and because different version of SIS computes those
coefficients in slightly different ways.
*/
private final double c2χ, c4χ, c6χ, c8χ;
@@ -140,10 +139,10 @@ abstract class ConformalProjection extends
NormalizedProjection {
* For each line, we add the smallest values first in order to reduce
rounding errors.
* The smallest values are the ones using the eccentricity raised to
the highest power.
*/
- c2χ = -73./ 2016 * e8 + 1./40 * e6 + 5./24 * e4 + 1./2 * e2;
- c4χ = 233./ 6720 * e8 + 29./120 * e6 + 7./24 * e4;
- c6χ = 81./ 280 * e8 + 7./30 * e6;
- c8χ = 4279./20160 * e8;
+ c2χ = fma(e2, 1./2, fma(e4, 5./24, fma(e6, 1./40, e8*( -73./2016))));
+ c4χ = fma(e4, 7./24, fma(e6, 29./120, e8*( 233./6720)));
+ c6χ = fma(e6, 7./30, e8*( 81./280));
+ c8χ = e8*(4279./20160);
}
/**
@@ -211,7 +210,11 @@ abstract class ConformalProjection extends
NormalizedProjection {
*/
final double sin_2φ = sin(2*φ);
final double cos_2φ = cos(2*φ);
- φ += sin_2φ * (c2χ + cos_2φ * (c4χ + cos_2φ * (c6χ + cos_2φ * c8χ)));
+ if (Formulas.USE_FMA) {
+ φ = fma(sin_2φ, fma(cos_2φ, fma(cos_2φ, fma(cos_2φ, c8χ, c6χ),
c4χ), c2χ), φ);
+ } else {
+ φ += sin_2φ * (c2χ + cos_2φ * (c4χ + cos_2φ * (c6χ + cos_2φ *
c8χ)));
+ }
/*
* Note: a previous version checked if the value of the smallest term
c8χ⋅sin(8φ) was smaller than
* the iteration tolerance. But this was not reliable enough. We use
now a hard coded threshold
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
index c6ea7da71a..5599e640e8 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/projection/LambertAzimuthalEqualArea.java
@@ -42,7 +42,7 @@ import static
org.apache.sis.internal.referencing.provider.LambertAzimuthalEqual
*
* @author Martin Desruisseaux (Geomatys)
* @author Rémi Maréchal (Geomatys)
- * @version 1.2
+ * @version 1.4
* @since 1.2
*/
public class LambertAzimuthalEqualArea extends AuthalicConversion {
@@ -133,16 +133,6 @@ public class LambertAzimuthalEqualArea extends
AuthalicConversion {
}
}
- /**
- * Creates a new projection initialized to the same parameters than the
given one.
- */
- LambertAzimuthalEqualArea(final LambertAzimuthalEqualArea other) {
- super(other);
- sinβ0 = other.sinβ0;
- cosβ0 = other.cosβ0;
- polar = other.polar;
- }
-
/**
* Returns the names of additional internal parameters which need to be
taken in account when
* comparing two {@code LambertAzimuthalEqualArea} projections or
formatting them in debug mode.
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 4c253979c1..688471ec9a 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
@@ -114,12 +114,12 @@ 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
*/
- cf0 = fma(e10, -441./0x10000, fma(e8, -175./0x4000, fma(e6,
-5./0x100, fma(e4, -3./0x40, fma(e2, -1./4, 1)))));
- cf1 = fma(e10, 1575./0x20000, fma(e8, 175./0x4000, fma(e6,
5./0x100, fma(e4, 3./0x40, e2*(-3./4)))));
- cf2 = fma(e10, -2625./0x8000, fma(e8, 175./0x6000, fma(e6,
5120./0x60000, e4*(-15./0x20))));
- cf3 = fma(e10, 735./0x800, fma(e8, 2240./0x60000, e6*
(-35./0x60)));
- cf4 = fma(e10, -2205./0x1000, e8*(-315./0x400));
- // cf5 = e10* (693./0x20000) omitted for now (not yet used).
+ cf0 = fma(e2, -1./4, fma(e4, -3./0x40, fma(e6, -5./0x100, fma(e8,
-175./0x4000, e10*( -441./0x10000))))) + 1;
+ cf1 = fma(e2, -3./4, fma(e4, 3./0x40, fma(e6, 5./0x100, fma(e8,
175./0x4000, e10*( 1575./0x20000)))));
+ cf2 = fma(e4, -15./0x20, fma(e6, 5120./0x60000, fma(e8,
175./0x6000, e10*(-2625./0x8000))));
+ cf3 = fma(e6, -35./0x60, fma(e8,
2240./0x60000, e10*( 735./0x800)));
+ cf4 = fma(e8,
-315./0x400, e10*(-2205./0x1000));
+ // cf5 = /* omitted for now (not yet used) */
e10*( 693./0x20000);
/*
* Coefficients for inverse transform derived from Snyder 3-26 and
EPSG guidance notes:
*
diff --git
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java
index 18414ac790..606d0d5f55 100644
---
a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java
+++
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform1D.java
@@ -80,6 +80,9 @@ class LinearTransform1D extends AbstractMathTransform1D
implements LinearTransfo
/**
* Error terms of {@link #scale} and {@link #offset} used in double-double
arithmetic.
+ * For performance reasons, those terms are not used in {@code
transform(…)} methods.
+ * They are used in {@link #getMatrix()} for making possible to use
double-double arithmetic
+ * during the creation of concatenated or inverse transforms.
*/
private final double scaleError, offsetError;
@@ -473,12 +476,12 @@ class LinearTransform1D extends AbstractMathTransform1D
implements LinearTransfo
doubleToRawLongBits(scaleError) ==
doubleToRawLongBits(that.scaleError) &&
doubleToRawLongBits(offsetError) ==
doubleToRawLongBits(that.offsetError);
/*
- * NOTE: 'LinearTransform1D' and 'ConstantTransform1D' are heavily
used by 'Category'
- * from 'org.apache.sis.coverage' package. It is essential for
Cateory to differenciate
- * various NaN values. Because 'equals' is used by
WeakHashSet.unique(Object) (which
- * is used by 'DefaultMathTransformFactory'), test for equality
cannot use the non-raw
- * doubleToLongBits method because it collapse all NaN into a
single canonical value.
- * The 'doubleToRawLongBits' method instead provides the needed
functionality.
+ * NOTE: `LinearTransform1D` and `ConstantTransform1D` are
extensively used by `SampleDimension`
+ * in `org.apache.sis.coverage` package. It is essential for
sample dimensions to differentiate
+ * various NaN values. Because `equals(…)` is used by
`WeakHashSet.unique(Object)` in turn used
+ * by `DefaultMathTransformFactory`, equality tests cannot use the
non-raw `doubleToLongBits(…)`
+ * method because it collapse all NaN into a single canonical
value.
+ * The `doubleToRawLongBits(…)` method instead provides the needed
functionality.
*/
}
return false;
diff --git
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
index 9d487a3617..3bf3efe63c 100644
---
a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
+++
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/projection/SinusoidalTest.java
@@ -123,7 +123,7 @@ public final class SinusoidalTest extends
MapProjectionTestCase {
createProjection(true);
final double delta = (1.0 / 60) / 1852; //
Approximately 1 metre.
derivativeDeltas = new double[] {delta, delta};
- tolerance = Formulas.LINEAR_TOLERANCE / 10;
+ tolerance = Formulas.LINEAR_TOLERANCE / 10000;
verifyDerivative(105, 30);
verifyDerivative(100, -60);
}
