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 2022/12/30 17:35:33 UTC
[sis] 02/02: More attemps to fix accuracy problems, in particular in `LinearTransform1D`. Include deprecation of `InterpolatedMolodenskyTransform`.
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 17ac6d41c2a3aa68b7c63fb0b8aa5fb42c1f6484
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Fri Dec 30 18:34:44 2022 +0100
More attemps to fix accuracy problems, in particular in `LinearTransform1D`.
Include deprecation of `InterpolatedMolodenskyTransform`.
---
.../apache/sis/internal/referencing/Formulas.java | 10 +-
.../sis/internal/referencing/package-info.java | 2 +-
.../provider/MolodenskyInterpolation.java | 3 +
.../operation/CoordinateOperationFinder.java | 19 ++-
.../operation/matrix/GeneralMatrix.java | 5 +
.../sis/referencing/operation/matrix/Matrices.java | 41 ++---
.../sis/referencing/operation/package-info.java | 2 +-
.../operation/projection/MeridianArcBased.java | 8 +-
.../transform/AbstractLinearTransform.java | 14 +-
.../operation/transform/ConstantTransform1D.java | 13 +-
.../operation/transform/ContextualParameters.java | 18 +-
.../transform/EllipsoidToCentricTransform.java | 1 +
.../transform/ExponentialTransform1D.java | 1 +
.../operation/transform/IdentityTransform1D.java | 4 +-
.../transform/InterpolatedMolodenskyTransform.java | 4 +
.../InterpolatedMolodenskyTransform2D.java | 3 +
.../operation/transform/LinearInterpolator1D.java | 8 +-
.../operation/transform/LinearTransform1D.java | 185 ++++++++++++++++-----
.../operation/transform/MathTransforms.java | 8 +-
.../operation/transform/PoleRotation.java | 1 +
.../operation/transform/ProjectiveTransform.java | 27 ++-
.../operation/CoordinateOperationFinderTest.java | 4 +-
.../referencing/operation/matrix/MatricesTest.java | 19 ++-
.../org/apache/sis/internal/util/Numerics.java | 8 -
.../org/apache/sis/measure/AbstractConverter.java | 2 +-
.../org/apache/sis/measure/LinearConverter.java | 27 ++-
.../main/java/org/apache/sis/measure/Units.java | 1 +
27 files changed, 305 insertions(+), 133 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 d6c41e07f1..2face4aff3 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
@@ -33,7 +33,7 @@ import static org.apache.sis.internal.metadata.ReferencingServices.NAUTICAL_MILE
* do not want to expose publicly those arbitrary values (or at least not in a too direct way).
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
* @since 0.4
*/
public final class Formulas extends Static {
@@ -98,6 +98,14 @@ public final class Formulas extends Static {
*/
public static final int MAXIMUM_ITERATIONS = 18;
+ /**
+ * Whether to use {@link Math#fma(double, double, double)} for performance reasons.
+ * We do not use this flag when the goal is to get better accuracy rather than performance.
+ * Use of FMA brings performance benefits on machines having hardware support,
+ * but come at a high cost on older machines without hardware support.
+ */
+ public static final boolean USE_FMA = true;
+
/**
* Do not allow instantiation of this class.
*/
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java
index 30f5f420b2..7d28e5d2fd 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/package-info.java
@@ -24,7 +24,7 @@
* may change in incompatible ways in any future version without notice.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
* @since 0.3
*/
package org.apache.sis.internal.referencing;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java
index 3f3a264f57..33f3eb118f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/MolodenskyInterpolation.java
@@ -42,8 +42,11 @@ import org.apache.sis.referencing.operation.transform.InterpolatedMolodenskyTran
* @author Martin Desruisseaux (Geomatys)
* @version 0.7
* @since 0.7
+ *
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-500">Deprecate (for removal) InterpolatedMolodenskyTransform</a>
*/
@XmlTransient
+@Deprecated(since="1.4", forRemoval=true)
public final class MolodenskyInterpolation extends FranceGeocentricInterpolation {
/**
* Serial number for inter-operability with different versions.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
index c0e620c9ab..436735e4d3 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/CoordinateOperationFinder.java
@@ -49,6 +49,7 @@ import org.apache.sis.internal.referencing.provider.GeographicToGeocentric;
import org.apache.sis.internal.referencing.provider.GeocentricToGeographic;
import org.apache.sis.internal.referencing.provider.GeocentricAffine;
import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.internal.util.Constants;
import org.apache.sis.measure.Units;
import org.apache.sis.metadata.iso.citation.Citations;
@@ -62,6 +63,7 @@ import org.apache.sis.referencing.cs.CoordinateSystems;
import org.apache.sis.referencing.datum.BursaWolfParameters;
import org.apache.sis.referencing.datum.DefaultGeodeticDatum;
import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.resources.Vocabulary;
@@ -113,7 +115,7 @@ import static org.apache.sis.util.Utilities.equalsIgnoreMetadata;
* </ul>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see DefaultCoordinateOperationFactory#createOperation(CoordinateReferenceSystem, CoordinateReferenceSystem, CoordinateOperationContext)
*
@@ -871,9 +873,9 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
* This "epoch shift" is in units of `targetCRS`.
*/
final Unit<Time> targetUnit = targetCS.getAxis(0).getUnit().asType(Time.class);
- double epochShift = sourceDatum.getOrigin().getTime() -
- targetDatum.getOrigin().getTime();
- epochShift = Units.MILLISECOND.getConverterTo(targetUnit).convert(epochShift);
+ DoubleDouble epochShift = new DoubleDouble(sourceDatum.getOrigin().getTime());
+ epochShift.subtract(new DoubleDouble(targetDatum.getOrigin().getTime()));
+ epochShift = DoubleDouble.castOrCopy(Units.MILLISECOND.getConverterTo(targetUnit).convert(epochShift));
/*
* Check axis directions. The method `swapAndScaleAxes` should returns a matrix of size 2×2.
* The element at index (0,0) may be +1 if source and target axes are in the same direction,
@@ -883,15 +885,16 @@ public class CoordinateOperationFinder extends CoordinateOperationRegistry {
*
* The "epoch shift" previously computed is a translation. Consequently, it is added to element (0,1).
*/
- final Matrix matrix;
+ final MatrixSIS matrix;
try {
- matrix = CoordinateSystems.swapAndScaleAxes(sourceCS, targetCS);
+ matrix = MatrixSIS.castOrCopy(CoordinateSystems.swapAndScaleAxes(sourceCS, targetCS));
} catch (IllegalArgumentException | IncommensurableException exception) {
throw new OperationNotFoundException(notFoundMessage(sourceCRS, targetCRS), exception);
}
final int translationColumn = matrix.getNumCol() - 1; // Paranoiac check: should always be 1.
- final double translation = matrix.getElement(0, translationColumn);
- matrix.setElement(0, translationColumn, translation + epochShift);
+ final DoubleDouble translation = DoubleDouble.castOrCopy(matrix.getNumber(0, translationColumn));
+ translation.add(epochShift);
+ matrix.setNumber(0, translationColumn, translation);
return asList(createFromAffineTransform(AXIS_CHANGES, sourceCRS, targetCRS, matrix));
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
index 0fba129a0c..b603a1e18d 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/GeneralMatrix.java
@@ -90,12 +90,17 @@ class GeneralMatrix extends MatrixSIS implements ExtendedPrecisionMatrix {
* If {@code setToIdentity} is {@code true}, then the elements
* on the diagonal (<var>j</var> == <var>i</var>) are set to 1.
*
+ * <p>Do not invoke this constructor directly (except by {@link NonSquareMatrix} constructor) unless
+ * the matrix is known to be square. If this is not the case, invoke a factory method instead.</p>
+ *
* @param numRow number of rows.
* @param numCol number of columns.
* @param setToIdentity {@code true} for initializing the matrix to the identity matrix,
* or {@code false} for leaving it initialized to zero.
* @param precision 1 for normal precision, or 2 for extended precision.
* No other value is allowed (this is not verified).
+ *
+ * @see #createExtendedPrecision(int, int, boolean)
*/
GeneralMatrix(final int numRow, final int numCol, final boolean setToIdentity, final int precision) {
ensureValidSize(numRow, numCol);
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
index c1f50e882c..b20281f0aa 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrices.java
@@ -253,7 +253,7 @@ public final class Matrices extends Static {
* enough because callers in other package may perform additional arithmetic operations
* on it (for example org.apache.sis.referencing.cs.CoordinateSystems.swapAndScaleAxes).
*/
- final MatrixSIS matrix = new GeneralMatrix(dstAxes.length+1, srcAxes.length+1, false, 2);
+ final MatrixSIS matrix = GeneralMatrix.createExtendedPrecision(dstAxes.length+1, srcAxes.length+1, false);
/*
* Maps source axes to destination axes. If no axis is moved (for example if the user
* want to transform (NORTH,EAST) to (SOUTH,EAST)), then source and destination index
@@ -404,10 +404,9 @@ public final class Matrices extends Static {
* then an exception will be thrown.</li>
* </ul>
*
- * <div class="note"><b>Example:</b>
- * it is legal to transform from (<i>easting</i>, <i>northing</i>, <i>up</i>) to
- * (<i>easting</i>, <i>northing</i>) — this is the first above case — but illegal
- * to transform (<i>easting</i>, <i>northing</i>) to (<i>easting</i>, <i>up</i>).</div>
+ * For example, it is legal to transform from (<i>easting</i>, <i>northing</i>, <i>up</i>)
+ * to (<i>easting</i>, <i>northing</i>) — this is the first above case — but illegal
+ * to transform (<i>easting</i>, <i>northing</i>) to (<i>easting</i>, <i>up</i>).
*
* <h4>Example</h4>
* The following method call:
@@ -444,8 +443,10 @@ public final class Matrices extends Static {
/*
* createTransform(…) may fail if the arrays contain two axes with the same direction, for example
* AxisDirection.OTHER. This check prevents that failure for the common case of an identity transform.
+ * The returned matrix must use extended precision for reason documented in `createTransform(…)`.
*/
- return Matrices.createIdentity(srcAxes.length + 1);
+ final int n = srcAxes.length + 1;
+ return new GeneralMatrix(n, n, true, 2);
}
return createTransform(null, srcAxes, null, dstAxes, false);
}
@@ -647,13 +648,13 @@ public final class Matrices extends Static {
* The 'stride' and 'length' values will be used for computing indices in that array.
* The DoubleDouble temporary object is used only if the array contains error terms.
*/
- final int stride = sourceDimensions;
- final int length = sourceDimensions * targetDimensions;
- final double[] sources = getExtendedElements(subMatrix);
- final DoubleDouble transfer = (sources.length > length) ? new DoubleDouble() : null;
- final MatrixSIS matrix = createZero(targetDimensions-- + expansion,
- sourceDimensions-- + expansion,
- transfer != null);
+ final int stride = sourceDimensions;
+ final int length = sourceDimensions * targetDimensions;
+ final double[] sources = getExtendedElements(subMatrix);
+ final var transfer = (sources.length > length) ? new DoubleDouble() : null;
+ final MatrixSIS matrix = createZero(targetDimensions-- + expansion,
+ sourceDimensions-- + expansion,
+ transfer != null);
/*
* Following code processes from upper row to lower row.
* First, set the diagonal elements on leading new dimensions.
@@ -788,13 +789,13 @@ public final class Matrices extends Static {
if (numRow == srcRow && numCol == srcCol) {
return matrix;
}
- final int stride = srcCol;
- final int length = srcCol * srcRow;
- final double[] sources = getExtendedElements(matrix);
- final DoubleDouble transfer = (sources.length > length) ? new DoubleDouble() : null;
- final MatrixSIS resized = createZero(numRow, numCol, transfer != null);
- final int copyRow = Math.min(--numRow, --srcRow);
- final int copyCol = Math.min(--numCol, --srcCol);
+ final int stride = srcCol;
+ final int length = srcCol * srcRow;
+ final double[] sources = getExtendedElements(matrix);
+ final var transfer = (sources.length > length) ? new DoubleDouble() : null;
+ final MatrixSIS resized = createZero(numRow, numCol, transfer != null);
+ final int copyRow = Math.min(--numRow, --srcRow);
+ final int copyCol = Math.min(--numCol, --srcCol);
for (int j=copyRow; j<numRow; j++) {
resized.setElement(j, j, 1);
}
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java
index 23a4dbf177..6732e798ed 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/package-info.java
@@ -91,7 +91,7 @@
* for example by specifying the area of interest.
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 1.3
+ * @version 1.4
* @since 0.6
*/
@XmlSchema(location = "http://schemas.opengis.net/gml/3.2.1/coordinateOperations.xsd",
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 b9c784e698..4c253979c1 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,8 @@
*/
package org.apache.sis.referencing.operation.projection;
-import org.apache.sis.internal.util.Numerics;
import org.apache.sis.internal.util.DoubleDouble;
+import org.apache.sis.internal.referencing.Formulas;
import static java.lang.Math.*;
@@ -186,7 +186,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
*/
final double distance(final double φ, final double sinφ, final double cosφ) {
final double sinφ2 = sinφ * sinφ;
- if (Numerics.USE_FMA) {
+ if (Formulas.USE_FMA) {
return fma(cosφ*sinφ, fma(sinφ2, fma(sinφ2, fma(sinφ2, cf4, cf3), cf2), cf1), cf0*φ);
} else {
return cf0*φ + cosφ*sinφ*(cf1 + sinφ2*(cf2 + sinφ2*(cf3 + sinφ2*cf4)));
@@ -199,7 +199,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
* @return the derivative at the specified latitude.
*/
final double dM_dφ(final double sinφ2) {
- if (Numerics.USE_FMA) {
+ if (Formulas.USE_FMA) {
return fma(fma(fma(fma(fma(-8, sinφ2, 7),
cf4, -6*cf3), sinφ2,
5*cf3 - 4*cf2), sinφ2,
@@ -225,7 +225,7 @@ abstract class MeridianArcBased extends NormalizedProjection {
double cosφ = cos(φ);
double sinφ = sin(φ);
double sinφ2 = sinφ * sinφ;
- if (Numerics.USE_FMA) {
+ if (Formulas.USE_FMA) {
φ = fma(cosφ*sinφ, fma(sinφ2, fma(sinφ2, fma(sinφ2, ci4, ci3), ci2), ci1), φ);
} else {
φ += cosφ*sinφ*(ci1 + sinφ2*(ci2 + sinφ2*(ci3 + sinφ2*ci4))); // Snyder 3-26.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
index 3533809e00..e41b38a184 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
@@ -27,6 +27,7 @@ import org.opengis.referencing.operation.NoninvertibleTransformException;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.internal.referencing.provider.Affine;
import org.apache.sis.internal.referencing.Resources;
+import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.internal.util.Numerics;
import org.apache.sis.util.ComparisonMode;
import org.apache.sis.util.resources.Errors;
@@ -34,8 +35,8 @@ import org.opengis.util.FactoryException;
/**
- * Base class of linear transforms. For efficiency reasons, this transform implements itself the matrix
- * to be returned by {@link #getMatrix()}.
+ * Base class of linear transforms.
+ * For efficiency reasons, this transform implements itself the matrix to be returned by {@link #getMatrix()}.
*
* <p>Subclasses need to implement the following methods:</p>
* <ul>
@@ -44,7 +45,7 @@ import org.opengis.util.FactoryException;
* </ul>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.1
+ * @version 1.4
* @since 0.6
*/
abstract class AbstractLinearTransform extends AbstractMathTransform implements LinearTransform, Matrix, Serializable {
@@ -60,6 +61,7 @@ abstract class AbstractLinearTransform extends AbstractMathTransform implements
*
* @see #inverse()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
volatile LinearTransform inverse;
/**
@@ -239,7 +241,11 @@ abstract class AbstractLinearTransform extends AbstractMathTransform implements
for (int i=0; i<srcDim; i++) {
final double e = getElement(j, i);
if (e != 0) { // See the comment in ProjectiveTransform for the purpose of this test.
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java
index 47388b935a..9b72b07479 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ConstantTransform1D.java
@@ -30,7 +30,7 @@ import java.util.Arrays;
* those values to NaN. Overwriting NaN by the constant value is required by {@link org.apache.sis.coverage}.</p>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.5
+ * @version 1.4
* @since 0.5
*/
final class ConstantTransform1D extends LinearTransform1D {
@@ -42,20 +42,21 @@ final class ConstantTransform1D extends LinearTransform1D {
/**
* A transform for the positive zero constant.
*/
- static final ConstantTransform1D ZERO = new ConstantTransform1D(0);
+ static final ConstantTransform1D ZERO = new ConstantTransform1D(0, 0);
/**
* A transform for the one constant.
*/
- static final ConstantTransform1D ONE = new ConstantTransform1D(1);
+ static final ConstantTransform1D ONE = new ConstantTransform1D(1, 0);
/**
* Constructs a new constant transform.
*
- * @param offset the {@code offset} term in the linear equation.
+ * @param offset the {@code offset} term in the linear equation.
+ * @param offsetError error term of {@code offset} for double-double arithmetic.
*/
- ConstantTransform1D(final double offset) {
- super(0, offset);
+ ConstantTransform1D(final double offset, final double offsetError) {
+ super(0, offset, 0, offsetError);
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java
index 22891e1f70..9c6c8ecc71 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ContextualParameters.java
@@ -124,7 +124,7 @@ import static java.util.logging.Logger.getLogger;
* Serialization should be used only for short term storage or RMI between applications running the same SIS version.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see org.apache.sis.referencing.operation.projection.NormalizedProjection
* @see AbstractMathTransform#getContextualParameters()
@@ -187,6 +187,7 @@ public class ContextualParameters extends Parameters implements Serializable {
*
* @see #getDescriptor()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private final ParameterDescriptorGroup descriptor;
/**
@@ -198,6 +199,7 @@ public class ContextualParameters extends Parameters implements Serializable {
*
* @see #getMatrix(MatrixRole)
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private Matrix normalize, denormalize;
/**
@@ -210,6 +212,7 @@ public class ContextualParameters extends Parameters implements Serializable {
* @see #parameter(String)
* @see #freeze()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private ParameterValue<?>[] values;
/**
@@ -300,19 +303,6 @@ public class ContextualParameters extends Parameters implements Serializable {
isFrozen = true;
}
- /**
- * Creates a matrix for a linear step of the transforms chain.
- * It is important that the matrices created here are instances of {@link MatrixSIS}, in order
- * to allow {@link #getMatrix(MatrixRole)} to return the reference to the (de)normalize matrices.
- */
- private static MatrixSIS linear(final String name, final Integer size) {
- if (size == null) {
- throw new IllegalArgumentException(Errors.format(Errors.Keys.MissingValueForProperty_1, name));
- }
- final int n = size + 1;
- return Matrices.create(n, n, ExtendedPrecisionMatrix.IDENTITY);
- }
-
/**
* Creates a new and frozen {@code ContextualParameters} for the inverse of this operation.
* An optional {@code mapper} can be specified for setting the parameter values.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java
index af458d6e1f..8ce552ab11 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/EllipsoidToCentricTransform.java
@@ -214,6 +214,7 @@ public class EllipsoidToCentricTransform extends AbstractMathTransform implement
* construction time). In addition this field is part of serialization form in order to preserve the
* references graph.</div>
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private final AbstractMathTransform inverse;
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java
index 247acb0f25..03512753b2 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ExponentialTransform1D.java
@@ -83,6 +83,7 @@ final class ExponentialTransform1D extends AbstractMathTransform1D implements Se
* The inverse of this transform. Created only when first needed. Serialized in order to avoid
* rounding error if this transform is actually the one which was created from the inverse.
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private MathTransform1D inverse;
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java
index 44872e4d40..6dbc23b56e 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/IdentityTransform1D.java
@@ -22,7 +22,7 @@ package org.apache.sis.referencing.operation.transform;
* This class is a special case of {@link LinearTransform1D} optimized for speed.
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.5
+ * @version 1.4
* @since 0.5
*/
final class IdentityTransform1D extends LinearTransform1D {
@@ -40,7 +40,7 @@ final class IdentityTransform1D extends LinearTransform1D {
* Constructs a new identity transform.
*/
private IdentityTransform1D() {
- super(1, 0);
+ super(1, 0, 0, 0);
}
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java
index 0a308e9876..1008bce1df 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform.java
@@ -59,13 +59,17 @@ import org.apache.sis.util.Debug;
* But in the {@code InterpolatedMolodenskyTransform} case, the interpolated translations are rather the
* ({@linkplain #tX}, {@linkplain #tY}, {@linkplain #tZ}) parameters of a Molodensky transformation.</p>
*
+ * @deprecated This operation method is non-standard, of little use and has greater errors than intended.
+ *
* @author Martin Desruisseaux (Geomatys)
* @version 1.0
*
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-500">Deprecate (for removal) InterpolatedMolodenskyTransform</a>
* @see InterpolatedGeocentricTransform
*
* @since 0.7
*/
+@Deprecated(since="1.4", forRemoval=true)
public class InterpolatedMolodenskyTransform extends MolodenskyFormula {
/**
* Serial number for inter-operability with different versions.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java
index 8b341bc1ce..5cc19d8ce7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedMolodenskyTransform2D.java
@@ -33,7 +33,10 @@ import org.apache.sis.referencing.datum.DatumShiftGrid;
* @author Martin Desruisseaux (Geomatys)
* @version 0.7
* @since 0.7
+ *
+ * @see <a href="https://issues.apache.org/jira/browse/SIS-500">Deprecate (for removal) InterpolatedMolodenskyTransform</a>
*/
+@Deprecated(since="1.4", forRemoval=true)
final class InterpolatedMolodenskyTransform2D extends InterpolatedMolodenskyTransform implements MathTransform2D {
/**
* Serial number for compatibility with different versions.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java
index b18873bf7d..c3faf5bdd7 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearInterpolator1D.java
@@ -54,7 +54,7 @@ import org.apache.sis.util.resources.Errors;
* @author Johann Sorel (Geomatys)
* @author Rémi Maréchal (Geomatys)
* @author Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.4
*
* @see MathTransforms#interpolate(double[], double[])
*
@@ -116,15 +116,15 @@ final class LinearInterpolator1D extends AbstractMathTransform1D implements Seri
* return a one-dimensional affine transform instead of an interpolator.
* We need to perform this check before the sign reversal applied after this loop.
*/
- double value;
+ double value, tolerance;
int i = 0;
do {
if (++i >= n) {
return LinearTransform1D.create(slope, offset);
}
value = values[i];
- } while (Numerics.epsilonEqual(value, Math.fma(i, slope, offset),
- Math.max(Math.abs(value), as) * Numerics.COMPARISON_THRESHOLD));
+ tolerance = Math.max(Math.abs(value), as) * Numerics.COMPARISON_THRESHOLD;
+ } while (Numerics.epsilonEqual(value, Math.fma(i, slope, offset), tolerance));
/*
* If the values are in decreasing order, reverse their sign so we get increasing order.
* We will multiply the results by -1 after the transformation.
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 72e8e9db50..18414ac790 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
@@ -26,8 +26,11 @@ import org.opengis.referencing.operation.TransformException;
import org.opengis.referencing.operation.NoninvertibleTransformException;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.Matrix1;
-import org.apache.sis.referencing.operation.matrix.Matrix2;
+import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
import org.apache.sis.internal.referencing.provider.Affine;
+import org.apache.sis.internal.referencing.Formulas;
+import org.apache.sis.internal.util.DoubleDouble;
+import org.apache.sis.util.resources.Errors;
import org.apache.sis.util.ComparisonMode;
/*
* We really want to use doubleToRawLongBits, not doubleToLongBits, because the
@@ -46,14 +49,15 @@ import static java.lang.Double.doubleToRawLongBits;
* is faster. This kind of transform is extensively used by {@code org.apache.sis.coverage.grid.GridCoverage2D}.
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 0.7
+ * @version 1.4
*
* @see LogarithmicTransform1D
* @see ExponentialTransform1D
*
* @since 0.5
*/
-class LinearTransform1D extends AbstractMathTransform1D implements LinearTransform, Serializable {
+@SuppressWarnings("CloneInNonCloneableClass")
+class LinearTransform1D extends AbstractMathTransform1D implements LinearTransform, ExtendedPrecisionMatrix, Serializable {
/**
* Serial number for inter-operability with different versions.
*/
@@ -62,7 +66,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
/**
* A transform that just reverse the sign of input values.
*/
- static final LinearTransform1D NEGATE = new LinearTransform1D(-1, 0);
+ static final LinearTransform1D NEGATE = new LinearTransform1D(-1, 0, 0, 0);
/**
* The value which is multiplied to input values.
@@ -74,6 +78,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
*/
final double offset;
+ /**
+ * Error terms of {@link #scale} and {@link #offset} used in double-double arithmetic.
+ */
+ private final double scaleError, offsetError;
+
/**
* The inverse of this transform. Created only when first needed.
*/
@@ -84,14 +93,36 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* Instances should be created using the {@linkplain #create(double, double) factory method},
* which may returns optimized implementations for some particular argument values.
*
- * @param scale the {@code scale} term in the linear equation.
- * @param offset the {@code offset} term in the linear equation.
+ * @param scale the {@code scale} term in the linear equation.
+ * @param offset the {@code offset} term in the linear equation.
+ * @param scaleError error term of {@code scale} for double-double arithmetic.
+ * @param offsetError error term of {@code offset} for double-double arithmetic.
*
* @see #create(double, double)
*/
- protected LinearTransform1D(final double scale, final double offset) {
- this.scale = scale;
- this.offset = offset;
+ LinearTransform1D(final double scale, final double offset, final double scaleError, final double offsetError) {
+ this.scale = scale;
+ this.offset = offset;
+ this.scaleError = scaleError;
+ this.offsetError = offsetError;
+ }
+
+ /**
+ * Constructs a new linear transform and remember error terms.
+ *
+ * @param scale the {@code scale} term in the linear equation.
+ * @param offset the {@code offset} term in the linear equation.
+ * @return the linear transform for the given scale and offset.
+ */
+ static LinearTransform1D create(final DoubleDouble scale, final DoubleDouble offset) {
+ if (scale.error == 0) {
+ if (offset.error == 0) {
+ return create(scale.value, offset.value);
+ } else if (scale.value == 0) {
+ return new ConstantTransform1D(offset.value, offset.error);
+ }
+ }
+ return new LinearTransform1D(scale.value, offset.value, scale.error, offset.error);
}
/**
@@ -111,9 +142,9 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
if (scale == 0) {
if (offset == 0) return ConstantTransform1D.ZERO;
if (offset == 1) return ConstantTransform1D.ONE;
- return new ConstantTransform1D(offset);
+ return new ConstantTransform1D(offset, 0);
}
- return new LinearTransform1D(scale, offset);
+ return new LinearTransform1D(scale, offset, 0, 0);
}
/**
@@ -129,6 +160,43 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
return tr;
}
+ /**
+ * Implementation of Matrix API. No need for clone because this matrix is immutable.
+ */
+ @SuppressWarnings("CloneDoesntCallSuperClone")
+ @Override public final Matrix clone() {return this;}
+ @Override public final Matrix getMatrix() {return this;}
+ @Override public final int getNumRow() {return 2;}
+ @Override public final int getNumCol() {return 2;}
+
+ /** Unsupported operation because this matrix shall be immutable. */
+ @Override public void setElement(int row, int column, double value) {
+ throw new UnsupportedOperationException(Errors.format(Errors.Keys.UnmodifiableObject_1, Matrix.class));
+ }
+
+ /**
+ * Retrieves the value at the specified row and column of the matrix.
+ * Row and column indices can be only 0 or 1.
+ */
+ @Override
+ public final double getElement(final int row, final int column) {
+ if (((row | column) & ~1) != 0) {
+ throw new IndexOutOfBoundsException();
+ }
+ return (row == 0) ? ((column == 0) ? scale : offset) : column;
+ }
+
+ /**
+ * Returns a copy of all matrix elements followed by the error terms for extended-precision arithmetic.
+ * Matrix elements are returned in a flat, row-major (column indices vary fastest) array.
+ *
+ * @return a copy of matrix elements followed by error terms.
+ */
+ @Override
+ public final double[] getExtendedElements() {
+ return new double[] {scale, offset, 0, 1, scaleError, offsetError, 0, 0};
+ }
+
/**
* Returns the parameter descriptors for this math transform.
*/
@@ -149,19 +217,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
return Affine.parameters(getMatrix());
}
- /**
- * Returns this transform as an affine transform matrix.
- */
- @Override
- public Matrix getMatrix() {
- return new Matrix2(scale, offset, 0, 1);
- }
-
/**
* Creates the inverse transform of this object.
*/
@Override
- public LinearTransform1D inverse() throws NoninvertibleTransformException {
+ public synchronized LinearTransform1D inverse() throws NoninvertibleTransformException {
if (inverse == null) {
/*
* Note: we do not perform the following optimization, because MathTransforms.linear(…)
@@ -173,7 +233,15 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
*/
if (scale != 0) {
final LinearTransform1D inverse;
- inverse = create(1/scale, -offset/scale);
+ if (DoubleDouble.DISABLED) {
+ inverse = create(1/scale, -offset/scale);
+ } else {
+ final var sd = new DoubleDouble(scale, scaleError);
+ final var od = new DoubleDouble(sd);
+ od.inverseDivide(-offset, -offsetError);
+ sd.inverseDivide(1, 0);
+ inverse = create(sd, od);
+ }
inverse.inverse = this;
this.inverse = inverse;
} else {
@@ -187,20 +255,19 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* Returns {@code true} since this transform is affine.
*/
@Override
- public boolean isAffine() {
+ public final boolean isAffine() {
return true;
}
/**
* Tests whether this transform does not move any points.
- *
- * <div class="note"><b>Note:</b> this method should always returns {@code false}, since
+ * This method should always returns {@code false}, because
* {@code MathTransforms.linear(…)} should have created specialized implementations for identity cases.
* Nevertheless we perform the full check as a safety, in case someone instantiated this class directly
- * instead of using a factory method.</div>
+ * instead of using a factory method.
*/
@Override
- public boolean isIdentity() {
+ public final boolean isIdentity() {
return offset == 0 && scale == 1;
}
@@ -211,7 +278,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* @return the derivative at the given point.
*/
@Override
- public Matrix derivative(final DirectPosition point) throws TransformException {
+ public final Matrix derivative(final DirectPosition point) throws TransformException {
return new Matrix1(scale);
}
@@ -222,7 +289,7 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* @return the derivative at the given point.
*/
@Override
- public double derivative(final double value) {
+ public final double derivative(final double value) {
return scale;
}
@@ -230,8 +297,12 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
* Transforms the specified value.
*/
@Override
- public double transform(double value) {
- return offset + scale * value;
+ public double transform(final double value) {
+ if (Formulas.USE_FMA) {
+ return Math.fma(value, scale, offset);
+ } else {
+ return offset + scale * value;
+ }
}
/**
@@ -245,7 +316,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
final boolean derivate)
{
if (dstPts != null) {
- dstPts[dstOff] = offset + scale*srcPts[srcOff];
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff] = Math.fma(srcPts[srcOff], scale, offset);
+ } else {
+ dstPts[dstOff] = offset + scale*srcPts[srcOff];
+ }
}
return derivate ? new Matrix1(scale) : null;
}
@@ -260,13 +335,21 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
{
if (srcPts != dstPts || srcOff >= dstOff) {
while (--numPts >= 0) {
- dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ }
}
} else {
srcOff += numPts;
dstOff += numPts;
while (--numPts >= 0) {
- dstPts[--dstOff] = offset + scale * srcPts[--srcOff];
+ if (Formulas.USE_FMA) {
+ dstPts[--dstOff] = Math.fma(srcPts[--srcOff], scale, offset);
+ } else {
+ dstPts[--dstOff] = offset + scale * srcPts[--srcOff];
+ }
}
}
}
@@ -282,13 +365,21 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
{
if (srcPts != dstPts || srcOff >= dstOff) {
while (--numPts >= 0) {
- dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ }
}
} else {
srcOff += numPts;
dstOff += numPts;
while (--numPts >= 0) {
- dstPts[--dstOff] = (float) (offset + scale * srcPts[--srcOff]);
+ if (Formulas.USE_FMA) {
+ dstPts[--dstOff] = (float) Math.fma(srcPts[--srcOff], scale, offset);
+ } else {
+ dstPts[--dstOff] = (float) (offset + scale * srcPts[--srcOff]);
+ }
}
}
}
@@ -303,7 +394,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
final float [] dstPts, int dstOff, int numPts)
{
while (--numPts >= 0) {
- dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = (float) Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = (float) (offset + scale * srcPts[srcOff++]);
+ }
}
}
@@ -316,7 +411,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
final double[] dstPts, int dstOff, int numPts)
{
while (--numPts >= 0) {
- dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ if (Formulas.USE_FMA) {
+ dstPts[dstOff++] = Math.fma(srcPts[srcOff++], scale, offset);
+ } else {
+ dstPts[dstOff++] = offset + scale * srcPts[srcOff++];
+ }
}
}
@@ -348,7 +447,11 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
*/
@Override
protected int computeHashCode() {
- return Long.hashCode(doubleToRawLongBits(offset) + 31*doubleToRawLongBits(scale)) ^ super.computeHashCode();
+ return Long.hashCode(super.computeHashCode() ^
+ (doubleToRawLongBits(offset) +
+ doubleToRawLongBits(scale) * 31 +
+ doubleToRawLongBits(offsetError) * 37 +
+ doubleToRawLongBits(scaleError) * 7));
}
/**
@@ -365,8 +468,10 @@ class LinearTransform1D extends AbstractMathTransform1D implements LinearTransfo
}
} else if (super.equals(object, mode)) {
final LinearTransform1D that = (LinearTransform1D) object;
- return doubleToRawLongBits(this.scale) == doubleToRawLongBits(that.scale) &&
- doubleToRawLongBits(this.offset) == doubleToRawLongBits(that.offset);
+ return doubleToRawLongBits(scale) == doubleToRawLongBits(that.scale) &&
+ doubleToRawLongBits(offset) == doubleToRawLongBits(that.offset) &&
+ 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
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
index e9a8bae8a2..ac484de9ef 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java
@@ -36,6 +36,7 @@ import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
import org.apache.sis.referencing.operation.matrix.AffineTransforms2D;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.referencing.operation.matrix.Matrices;
+import org.apache.sis.internal.util.DoubleDouble;
import org.apache.sis.util.OptionalCandidate;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.ArraysExt;
@@ -57,7 +58,7 @@ import org.apache.sis.util.Static;
* GeoAPI factory interfaces instead.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.3
+ * @version 1.4
*
* @see MathTransformFactory
*
@@ -198,7 +199,10 @@ public final class MathTransforms extends Static {
if (Matrices.isAffine(matrix)) {
switch (sourceDimension) {
case 1: {
- return linear(matrix.getElement(0,0), matrix.getElement(0,1));
+ final MatrixSIS m = MatrixSIS.castOrCopy(matrix);
+ return LinearTransform1D.create(
+ DoubleDouble.castOrCopy(m.getNumber(0,0)),
+ DoubleDouble.castOrCopy(m.getNumber(0,1)));
}
case 2: {
if (matrix instanceof ExtendedPrecisionMatrix) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java
index 3fd77bdec3..a445b93e44 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/PoleRotation.java
@@ -126,6 +126,7 @@ public class PoleRotation extends AbstractMathTransform2D implements Serializabl
*
* @see #inverse()
*/
+ @SuppressWarnings("serial") // Not statically typed as Serializable.
private MathTransform2D inverse;
/**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
index 1ae2d841df..ed8506e906 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
@@ -23,6 +23,7 @@ import org.apache.sis.internal.referencing.DirectPositionView;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
+import org.apache.sis.internal.referencing.Formulas;
import org.apache.sis.util.ArgumentChecks;
@@ -35,7 +36,7 @@ import org.apache.sis.util.ArgumentChecks;
* lines in the source is preserved in the output.</p>
*
* @author Martin Desruisseaux (IRD, Geomatys)
- * @version 1.2
+ * @version 1.4
*
* @see java.awt.geom.AffineTransform
*
@@ -281,7 +282,11 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
* getting 2D points from 3D points. In such case, the fact that the excluded dimensions had
* NaN values should not force the retained dimensions to get NaN values.
*/
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
@@ -345,7 +350,11 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
for (int i=0; i<srcDim; i++) {
final double e = elt[mix++];
if (e != 0) { // See comment in transform(double[], ...)
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
@@ -381,7 +390,11 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
for (int i=0; i<srcDim; i++) {
final double e = elt[mix++];
if (e != 0) { // See comment in transform(double[], ...)
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
@@ -416,7 +429,11 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
for (int i=0; i<srcDim; i++) {
final double e = elt[mix++];
if (e != 0) { // See comment in transform(double[], ...)
- sum += srcPts[srcOff + i] * e;
+ if (Formulas.USE_FMA) {
+ sum = Math.fma(srcPts[srcOff + i], e, sum);
+ } else {
+ sum += srcPts[srcOff + i] * e;
+ }
}
}
buffer[j] = sum;
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java
index 0fe369f82f..e897d9feab 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/CoordinateOperationFinderTest.java
@@ -663,7 +663,7 @@ public final class CoordinateOperationFinderTest extends MathTransformTestCase {
assertInstanceOf("operation", Conversion.class, operation);
transform = operation.getMathTransform();
- tolerance = 1E-12;
+ tolerance = 2E-12;
verifyTransform(new double[] {
// December 31, 1899 at 12:00 UTC in seconds.
CommonCRS.Temporal.DUBLIN_JULIAN.datum().getOrigin().getTime() / 1000
@@ -990,7 +990,7 @@ public final class CoordinateOperationFinderTest extends MathTransformTestCase {
0, 0, 0, 1
}), ((LinearTransform) transform).getMatrix(), 1E-12);
- tolerance = 1E-12;
+ tolerance = 2E-12;
verifyTransform(new double[] {
-5, -8, CommonCRS.Temporal.DUBLIN_JULIAN.datum().getOrigin().getTime() / 1000
}, new double[] {
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java
index 64bc647065..58d455e271 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/matrix/MatricesTest.java
@@ -41,7 +41,7 @@ import static org.opengis.referencing.cs.AxisDirection.*;
*
* @author Martin Desruisseaux (Geomatys)
* @author Johann Sorel (Geomatys)
- * @version 1.1
+ * @version 1.4
* @since 0.4
*/
@DependsOn({
@@ -106,6 +106,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {NORTH, EAST, UP});
+ assertExtendedPrecision(matrix);
assertTrue ("isAffine", matrix.isAffine());
assertTrue ("isIdentity", matrix.isIdentity());
assertEquals("numRow", 4, matrix.getNumRow());
@@ -128,6 +129,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {WEST, UP, SOUTH});
+ assertExtendedPrecision(matrix);
assertTrue ("isAffine", matrix.isAffine());
assertFalse("isIdentity", matrix.isIdentity());
assertEquals("numRow", 4, matrix.getNumRow());
@@ -155,6 +157,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {DOWN, NORTH});
+ assertExtendedPrecision(matrix);
assertFalse("isIdentity", matrix.isIdentity());
assertEquals("numRow", 3, matrix.getNumRow());
assertEquals("numCol", 4, matrix.getNumCol());
@@ -180,6 +183,7 @@ public final class MatricesTest extends TestCase {
new AxisDirection[] {NORTH, EAST, UP},
new AxisDirection[] {DOWN, DOWN});
+ assertExtendedPrecision(matrix);
assertFalse("isIdentity", matrix.isIdentity());
assertEquals("numRow", 3, matrix.getNumRow());
assertEquals("numCol", 4, matrix.getNumCol());
@@ -242,6 +246,19 @@ public final class MatricesTest extends TestCase {
}
}
+ /**
+ * Asserts that the given matrix uses extended precision. This is mandatory for all matrices
+ * returned by {@link Matrices#createTransform(AxisDirection[], AxisDirection[])} because
+ * {@code CoordinateSystems.swapAndScaleAxes(…)} will modify those matrices in-place with
+ * the assumption that they accept extended precision.
+ *
+ * @param matrix the matrix to test.
+ */
+ private static void assertExtendedPrecision(final Matrix matrix) {
+ assertTrue(matrix instanceof GeneralMatrix);
+ assertTrue(((GeneralMatrix) matrix).isExtendedPrecision());
+ }
+
/**
* Tests {@link Matrices#createTransform(Envelope, Envelope)}.
* This method tests the example given in {@code Matrices.createTransform(…)} javadoc.
diff --git a/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java b/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java
index f4f2c9d68a..7c2e038a4c 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/internal/util/Numerics.java
@@ -74,14 +74,6 @@ public final class Numerics extends Static {
boxed = -value; CACHE.put(boxed, boxed);
}
- /**
- * Whether to use {@link Math#fma(double, double, double)} for performance reasons.
- * We do not use this flag when the goal is to get better accuracy rather than performance.
- * Use of FMA brings performance benefits on machines having hardware support,
- * but come at a high cost on older machines without hardware support.
- */
- public static final boolean USE_FMA = true;
-
/**
* Maximum number of rows or columns in Apache SIS matrices. We define a maximum because SIS is expected to work
* mostly with small matrices, because their sizes are related to the number of dimensions in coordinate systems.
diff --git a/core/sis-utility/src/main/java/org/apache/sis/measure/AbstractConverter.java b/core/sis-utility/src/main/java/org/apache/sis/measure/AbstractConverter.java
index 6005d0aef9..e719e227ca 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/measure/AbstractConverter.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/measure/AbstractConverter.java
@@ -64,7 +64,7 @@ abstract class AbstractConverter implements UnitConverter, Serializable {
/**
* If the conversion can be represented by a polynomial equation, returns the coefficients of that equation.
- * Otherwise returns {@code null}.
+ * Otherwise returns {@code null}. This is the implementation of {@link Units#coefficients(UnitConverter)}.
*/
Number[] coefficients() {
return isIdentity() ? new Number[0] : null;
diff --git a/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java b/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java
index d75cc6884a..355b2b495a 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/measure/LinearConverter.java
@@ -26,6 +26,7 @@ import org.apache.sis.util.LenientComparable;
import org.apache.sis.math.MathFunctions;
import org.apache.sis.math.Fraction;
import org.apache.sis.internal.util.Numerics;
+import org.apache.sis.internal.util.DoubleDouble;
/**
@@ -239,6 +240,7 @@ final class LinearConverter extends AbstractConverter implements LenientComparab
/**
* Returns the coefficient of this linear conversion.
+ * Coefficients are offset and scale, in that order.
*/
@Override
@SuppressWarnings("fallthrough")
@@ -258,11 +260,13 @@ final class LinearConverter extends AbstractConverter implements LenientComparab
* package to perform more accurate calculations.
*/
private static Number ratio(final double value, final double divisor) {
- final int numerator = (int) value;
- if (numerator == value) {
- final int denominator = (int) divisor;
- if (denominator == divisor) {
- return (denominator == 1) ? numerator : new Fraction(numerator, denominator);
+ if (value != 0) {
+ final int numerator = (int) value;
+ if (numerator == value) {
+ final int denominator = (int) divisor;
+ if (denominator == divisor) {
+ return (denominator == 1) ? numerator : new Fraction(numerator, denominator);
+ }
}
}
return value / divisor;
@@ -278,15 +282,20 @@ final class LinearConverter extends AbstractConverter implements LenientComparab
/**
* Applies the linear conversion on the given value. This method uses {@link BigDecimal} arithmetic if
- * the given value is an instance of {@code BigDecimal}, or IEEE 754 floating-point arithmetic otherwise.
- *
- * <p>Apache SIS rarely uses {@link BigDecimal} arithmetic, so providing an efficient implementation of
- * this method is not a goal.</p>
+ * the given value is an instance of {@code BigDecimal}, or double-double arithmetic if the value is an
+ * instance of {@link DoubleDouble}, or IEEE 754 floating-point arithmetic otherwise.
*/
@Override
public Number convert(Number value) {
ArgumentChecks.ensureNonNull("value", value);
if (!isIdentity()) {
+ if (value instanceof DoubleDouble) {
+ var dd = new DoubleDouble((DoubleDouble) value);
+ dd.multiply(scale);
+ dd.add(offset);
+ dd.divide(divisor);
+ return dd;
+ }
if (value instanceof BigInteger) {
value = new BigDecimal((BigInteger) value);
}
diff --git a/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java b/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java
index 7327047740..1bec091089 100644
--- a/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java
+++ b/core/sis-utility/src/main/java/org/apache/sis/measure/Units.java
@@ -1654,6 +1654,7 @@ public final class Units extends Static {
*
* @since 0.8
*/
+ @OptionalCandidate
@SuppressWarnings("fallthrough")
public static Number[] coefficients(final UnitConverter converter) {
if (converter != null) {