#### Design note

+ * A previous version of {@code WraparoundTransform} had no period. It was expecting coordinates + * normalized in the [0 … 1] range. Coordinates in [0 … 360]° range were divided by 360 using an + * affine transforms before {@code WraparoundTransform} and multiplied by 360 using another affine + * transform after {@code WraparoundTransform}. That approach allowed to delegate more work to the + * affine transforms which can efficiently be combined with other affine transforms, but it caused + * more rounding errors. + */ + private final double period; + + /** * Creates a new transform with a wraparound behavior in the given dimension. * Input and output values in the wraparound dimension shall be normalized in - * the [0 … 1] range. + * the [−p/2 … +p/2] range where p is the period (e.g. 360°). */ - private WraparoundTransform(final int dimension, final int wraparoundDimension) { + private WraparoundTransform(final int dimension, final int wraparoundDimension, final double period) { this.dimension = dimension; this.wraparoundDimension = wraparoundDimension; + this.period = period; } /** - * Returns an instance with the given number of dimensions while keeping {@link #wraparoundDimension} unchanged. + * Returns an instance with the number of dimensions compatible with the given matrix, while keeping + * {@link #wraparoundDimension} and {@link #period} unchanged. * If no instance can be created for the given number of dimensions, then this method returns {@code null}. + * + * @param other matrix defining a transform which will be applied before or after {@code this}. + * @param applyOtherFirst whether the transform defined by the matrix will be applied before {@code this}. */ - private WraparoundTransform redim(final int n) { + private WraparoundTransform redim(final boolean applyOtherFirst, final Matrix other) { + final int n = (applyOtherFirst ? other.getNumRow() : other.getNumCol()) - 1; if (n == dimension) return this; if (n >= wraparoundDimension) return null; - return create(n, wraparoundDimension); + return create(n, wraparoundDimension, period); } /** @@ -136,7 +166,7 @@ public final class WraparoundTransform extends AbstractMathTransform { * on an inspection of source and target CRS axes. For a method making decision based on a domain of use, * see {@link #forDomainOfUse forDomainOfUse(…)} instead.

* - * @param op the coordinate operation for which to get the math transform. + * @param op the coordinate operation for which to get the math transform. * @return the math transform for the given coordinate operation. * @throws TransformException if a coordinate can not be computed. */ @@ -193,42 +223,51 @@ public final class WraparoundTransform extends AbstractMathTransform { private static MathTransform concatenate(final MathTransform tr, final int wraparoundDimension, final CoordinateSystem targetCS, final DirectPosition median) throws TransformException { - final double span = WraparoundAdjustment.range(targetCS, wraparoundDimension); - if (!(span > 0 && span != Double.POSITIVE_INFINITY)) { + final double period = WraparoundAdjustment.range(targetCS, wraparoundDimension); + if (!(period > 0 && period != Double.POSITIVE_INFINITY)) { return tr; } - double start; - if (median != null) { - try { - start = median.getOrdinate(wraparoundDimension) + -0.5 * span; - } catch (BackingStoreException e) { - // Some implementations compute coordinates only when first needed. - throw e.unwrapOrRethrow(TransformException.class); - } - if (!Double.isFinite(start)) return tr; + double m; + if (median == null) { + final CoordinateSystemAxis axis = targetCS.getAxis(wraparoundDimension); + m = (axis.getMinimumValue() + axis.getMaximumValue()) / 2; + } else try { + m = median.getOrdinate(wraparoundDimension); + } catch (BackingStoreException e) { + // Some implementations compute coordinates only when first needed. + throw e.unwrapOrRethrow(TransformException.class); + } + if (Double.isFinite(m)) { + /* + * Round the median to a value having an exact representation in base 2 using about 10 bits. + * The intent is to reduce the risk of rounding errors with add/subtract operations. + */ + final int power = 10 - Math.getExponent(m); + m = Math.scalb(Math.rint(Math.scalb(m, power)), -power); + } else if (median == null) { + /* + * May happen if `WraparoundAdjustment.range(…)` recognized a longitude axis + * despite the `CoordinateSystemAxis` not declarining minimum/maximum values. + * Use 0 as the range center (e.g. center of [-180 … 180]° longitude range). + */ + m = 0; } else { - start = targetCS.getAxis(wraparoundDimension).getMinimumValue(); - if (!Double.isFinite(start)) { - /* - * May happen if `WraparoundAdjustment.range(…)` recognized a longitude axis - * despite the `CoordinateSystemAxis` not declarining minimum/maximum values. - * Use 0 as the range center (e.g. center of [-180 … 180]° longitude range). - */ - start = -0.5 * span; - } + // Invalid median value. Assume caller means "no wrap". + return tr; } final int dimension = tr.getTargetDimensions(); - final MatrixSIS m = Matrices.createIdentity(dimension + 1); - m.setElement(wraparoundDimension, wraparoundDimension, span); - m.setElement(wraparoundDimension, dimension, start); - final MathTransform denormalize = MathTransforms.linear(m); - MathTransform wraparound = create(dimension, wraparoundDimension); - /* - * Do not use the 3-arguments method because we need to - * control the order in which concatenation is applied. - */ - wraparound = MathTransforms.concatenate(denormalize.inverse(), - MathTransforms.concatenate(wraparound, denormalize)); + MathTransform wraparound = create(dimension, wraparoundDimension, period); + if (m != 0) { + final double[] t = new double[dimension]; + t[wraparoundDimension] = m; + final MathTransform denormalize = MathTransforms.translation(t); + /* + * Do not use the 3-arguments method because we need to + * control the order in which concatenation is applied. + */ + wraparound = MathTransforms.concatenate(denormalize.inverse(), + MathTransforms.concatenate(wraparound, denormalize)); + } return MathTransforms.concatenate(tr, wraparound); } @@ -253,29 +292,19 @@ public final class WraparoundTransform extends AbstractMathTransform { } /** - * Whether to allow discontinuities in transform derivatives. Strictly speaking the derivative - * at 1, 2, 3… should be infinite because the coordinate value jumps "instantaneously" from an - * integer value to 0. However in practice we use derivatives as linear approximations around - * small regions, not for calculations requiring strict mathematical values. An infinite value - * goes against the approximation goal. Furthermore whether a source coordinate is an integer - * value or not is subject to rounding errors, which may cause unpredictable behavior if - * discontinuities are allowed. - */ - private static final boolean ALLOW_DISCONTINUITIES = false; - - /** * Gets the derivative of this transform at a point. + * + *

#### Mathematical note

+ * Strictly speaking the derivative at (n + 0.5) × {@link #period} where n is an integer + * should be infinite because the coordinate value jumps "instantaneously" from any value to ±{@link #period}/2. + * However in practice we use derivatives as linear approximations around small regions, not for calculations + * requiring strict mathematical values. An infinite value goes against the approximation goal. + * Furthermore whether a source coordinate is an integer value or not is subject to rounding errors, + * which may cause unpredictable behavior if those infinite values were provided. */ @Override public Matrix derivative(final DirectPosition point) { - final MatrixSIS derivative = Matrices.createIdentity(dimension); - if (ALLOW_DISCONTINUITIES) { - final double v = point.getOrdinate(wraparoundDimension); - if (v == Math.floor(v)) { - derivative.setElement(wraparoundDimension, wraparoundDimension, Double.NEGATIVE_INFINITY); - } - } - return derivative; + return Matrices.createIdentity(dimension); } /** @@ -283,23 +312,15 @@ public final class WraparoundTransform extends AbstractMathTransform { * and optionally computes the transform derivative at that location. */ @Override - public Matrix transform(final double[] srcPts, final int srcOff, - final double[] dstPts, final int dstOff, final boolean derivate) + public Matrix transform(final double[] srcPts, int srcOff, + final double[] dstPts, int dstOff, final boolean derivate) { - double v = srcPts[srcOff]; - v -= Math.floor(v); if (dstPts != null) { System.arraycopy(srcPts, srcOff, dstPts, dstOff, dimension); - dstPts[dstOff + wraparoundDimension] = v; - } - if (!derivate) { - return null; + dstOff += wraparoundDimension; + dstPts[dstOff] = Math.IEEEremainder(dstPts[dstOff], period); } - final MatrixSIS derivative = Matrices.createIdentity(dimension); - if (ALLOW_DISCONTINUITIES && v == 0) { - derivative.setElement(wraparoundDimension, wraparoundDimension, Double.NEGATIVE_INFINITY); - } - return derivative; + return derivate ? Matrices.createIdentity(dimension) : null; } /** @@ -312,7 +333,7 @@ public final class WraparoundTransform extends AbstractMathTransform { System.arraycopy(srcPts, srcOff, dstPts, dstOff, numPts * dimension); dstOff += wraparoundDimension; while (--numPts >= 0) { - dstPts[dstOff] -= Math.floor(dstPts[dstOff]); + dstPts[dstOff] = Math.IEEEremainder(dstPts[dstOff], period); dstOff += dimension; } } @@ -327,7 +348,7 @@ public final class WraparoundTransform extends AbstractMathTransform { System.arraycopy(srcPts, srcOff, dstPts, dstOff, numPts * dimension); dstOff += wraparoundDimension; while (--numPts >= 0) { - dstPts[dstOff] -= Math.floor(dstPts[dstOff]); + dstPts[dstOff] = (float) Math.IEEEremainder(dstPts[dstOff], period); dstOff += dimension; } } @@ -373,66 +394,85 @@ public final class WraparoundTransform extends AbstractMathTransform { if (equals(other, null)) { return this; } - /* - * If two `WraparoundTransform` instances are separated only by a `LinearTransform`, - * then that linear transform is moved before or after this `WraparoundTransform` - * for increasing the chances to concatenate it using a matrix multiplication. - */ - if (applyOtherFirst) { - final List steps = MathTransforms.getSteps(other); - if (steps.size() == 2) { - final MathTransform middle = steps.get(1); - Matrix matrix = MathTransforms.getMatrix(middle); - if (matrix != null) try { - final MathTransformsOrFactory mf = MathTransformsOrFactory.wrap(factory); - boolean modified = false; - MathTransform step2 = this; - final MathTransform after = movable(matrix, mf); - if (after != null) { - /* - * Update the middle matrix with everything that we could not put in `after`. - * Usually the matrix is square before the multiplication. But if it was not the case, - * the new matrix will have the same number of columns (source coordinates) but a new - * number of rows (target coordinates). The result should be a square matrix. - */ - final Matrix remaining = Matrices.multiply(MathTransforms.getMatrix(after.inverse()), matrix); - final WraparoundTransform redim = redim(remaining.getNumRow() - 1); + final List steps = MathTransforms.getSteps(other); + final int count = steps.size(); + if (count >= 2) { + Matrix middle = MathTransforms.getMatrix(steps.get(applyOtherFirst ? count - 1 : 0)); + if (middle != null) try { + /* + * We have a matrix between this `WraparoundTransform` and something else, + * potentially another `WraparoundTransform`. Try to move as many rows as + * possible outside that `middle` matrix. Ideally we will be able to move + * the matrix completely, which increase the chances to multiply (outside + * this method) with another matrix. + */ + final MathTransformsOrFactory mf = MathTransformsOrFactory.wrap(factory); + boolean modified = false; + MathTransform step1 = this; + MathTransform move = movable(middle, mf); + if (move != null) { + /* + * Update the middle matrix with everything that we could not put in `move`. + * If `applyOtherFirst` is false: + * + * [move] → [redim] → [remaining] → [other] + * + * If `applyOtherFirst` is true: + * + * [other] → [remaining] → [redim] → [move] + * + * Usually the matrix is square before the multiplication. But if it was not the case, + * the new matrix will have the same number of columns (source coordinates) but a new + * number of rows (target coordinates). The result should be a square matrix. + */ + final Matrix remaining = remaining(applyOtherFirst, move, middle); + final WraparoundTransform redim = redim(applyOtherFirst, remaining); + if (redim != null) { + step1 = mf.concatenate(applyOtherFirst, move, redim); + middle = remaining; + modified = true; + } + } + /* + * Now look at the non-linear transform. If it is another instance of `WraparoundTransform`, + * then we may move the calculation of some coordinates before it. This is the same algorithm + * than above but with `applyOtherFirst` reversed. + */ + MathTransform step2 = steps.get(applyOtherFirst ? count - 2 : 1); + if (step2 instanceof WraparoundTransform) { + WraparoundTransform redim = (WraparoundTransform) step2; + move = redim.movable(middle, mf); + if (move != null) { + final Matrix remaining = remaining(!applyOtherFirst, move, middle); + redim = redim.redim(!applyOtherFirst, remaining); if (redim != null) { - step2 = mf.concatenate(redim, after); - matrix = remaining; + step2 = mf.concatenate(applyOtherFirst, redim, move); + middle = remaining; modified = true; } } - /* - * Now look at the non-linear transform. If it is another instance of `WraparoundTransform`, - * then we may move the calculation of some coordinates before it. - */ - MathTransform step1 = steps.get(0); - if (step1 instanceof WraparoundTransform) { - WraparoundTransform wb = (WraparoundTransform) step1; - final MathTransform before = wb.movable(matrix, mf); - if (before != null) { - final Matrix remaining = Matrices.multiply(matrix, MathTransforms.getMatrix(before.inverse())); - wb = wb.redim(remaining.getNumCol() - 1); - if (wb != null) { - step1 = mf.concatenate(before, wb); - matrix = remaining; - modified = true; - } + } + /* + * Done moving the linear operations that we can move. Put everything together. + */ + if (modified) { + MathTransform tr = mf.linear(middle); + tr = mf.concatenate(applyOtherFirst, tr, step2); + tr = mf.concatenate(applyOtherFirst, step1, tr); + if (applyOtherFirst) { + for (int i = count-2; --i >= 0;) { + tr = mf.concatenate(steps.get(i), tr); + } + } else { + for (int i = 2; i < count; i++) { + tr = mf.concatenate(tr, steps.get(i)); } } - /* - * Done moving the linear operations that we can move. Put everything together. - * Do not use the 3-arguments method because we need to control the order in which - * concatenation is applied. - */ - if (modified) { - return mf.concatenate(mf.concatenate(step1, mf.linear(matrix)), step2); - } - } catch (NoninvertibleTransformException e) { - // Should not happen. But if it is the case, just abandon the optimization effort. - Logging.recoverableException(Logging.getLogger(Modules.REFERENCING), getClass(), "tryConcatenate", e); + return tr; } + } catch (NoninvertibleTransformException e) { + // Should not happen. But if it is the case, just abandon the optimization effort. + Logging.recoverableException(Logging.getLogger(Modules.REFERENCING), getClass(), "tryConcatenate", e); } } return null; @@ -483,6 +523,21 @@ public final class WraparoundTransform extends AbstractMathTransform { } /** + * Returns a matrix equivalent to applying the inverse of {@code move} first, followed by {@code middle}. + * This is appropriate if the {@code move} transform is going to be applied first, followed by {@code remaining}. + * + *

If {@code reverse} is {@code true} then this method performs the concatenation in reverse order. + * This is appropriate if the {@code remaining} matrix is going to be applied first, followed by {@code move}.

+ */ + private static Matrix remaining(final boolean reverse, final MathTransform move, final Matrix middle) + throws NoninvertibleTransformException + { + final Matrix change = MathTransforms.getMatrix(move.inverse()); + return reverse ? Matrices.multiply(change, middle) + : Matrices.multiply(middle, change); + } + + /** * Formats this transform as a pseudo-WKT element. * * @param formatter the formatter to use. @@ -490,7 +545,9 @@ public final class WraparoundTransform extends AbstractMathTransform { */ @Override protected String formatTo(final Formatter formatter) { + formatter.append(dimension); formatter.append(wraparoundDimension); + formatter.append(period); formatter.setInvalidWKT(WraparoundTransform.class, null); return "Wraparound_MT"; } @@ -505,7 +562,8 @@ public final class WraparoundTransform extends AbstractMathTransform { public boolean equals(final Object object, final ComparisonMode mode) { if (object instanceof WraparoundTransform) { final WraparoundTransform other = (WraparoundTransform) object; - return other.dimension == dimension && other.wraparoundDimension == wraparoundDimension; + return other.dimension == dimension && other.wraparoundDimension == wraparoundDimension + && Numerics.equals(period, other.period); } return false; } @@ -515,6 +573,6 @@ public final class WraparoundTransform extends AbstractMathTransform { */ @Override protected int computeHashCode() { - return dimension * 31 + wraparoundDimension; + return dimension * 31 + wraparoundDimension + Double.hashCode(period); } } diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java index 0aad7e3..dcdf148 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java @@ -47,11 +47,12 @@ public final strictfp class WraparoundTransformTest extends TestCase { */ @Test public void testCache() { + final double period = 360; final WraparoundTransform t1, t2, t3, t4; - assertSame (WraparoundTransform.create(3, 0), t1 = WraparoundTransform.create(3, 0)); - assertNotSame(WraparoundTransform.create(3, 0), t2 = WraparoundTransform.create(3, 1)); - assertNotSame(WraparoundTransform.create(3, 0), t3 = WraparoundTransform.create(2, 0)); - assertNotSame(WraparoundTransform.create(3, 0), t4 = WraparoundTransform.create(3, 2)); + assertSame (WraparoundTransform.create(3, 0, period), t1 = WraparoundTransform.create(3, 0, period)); + assertNotSame(WraparoundTransform.create(3, 0, period), t2 = WraparoundTransform.create(3, 1, period)); + assertNotSame(WraparoundTransform.create(3, 0, period), t3 = WraparoundTransform.create(2, 0, period)); + assertNotSame(WraparoundTransform.create(3, 0, period), t4 = WraparoundTransform.create(3, 2, period)); assertEquals(3, t1.getSourceDimensions()); assertEquals(3, t2.getSourceDimensions()); assertEquals(2, t3.getSourceDimensions()); @@ -71,8 +72,9 @@ public final strictfp class WraparoundTransformTest extends TestCase { public void testOneAxis() throws TransformException { final AbstractCoordinateOperation op = new AbstractCoordinateOperation( Collections.singletonMap(AbstractCoordinateOperation.NAME_KEY, "Wrapper"), + HardCodedCRS.WGS84_φλ, HardCodedCRS.WGS84_φλ.forConvention(AxesConvention.POSITIVE_RANGE), - HardCodedCRS.WGS84_φλ, null, MathTransforms.scale(3, 5)); + null, MathTransforms.scale(3, 5)); /* * Transform should be [scale & normalization] → [wraparound] → [denormalization]. * The wraparound is applied on target coordinates, which is why it appears after [scale]. @@ -84,27 +86,41 @@ public final strictfp class WraparoundTransformTest extends TestCase { assertEquals(3, steps.size()); assertEquals(1, ((WraparoundTransform) steps.get(1)).wraparoundDimension); /* - * Wraparound outputs are in [0 … 1) range (0 inclusive and 1 exclusive), so we expect a - * multiplication by the span of each axis for getting the final range. + * WraparoundTransform outputs are in [−180 … 180] range, so we expect + * a 180° shift for getting results in the [0 … 360]° range. */ assertMatrixEquals("denormalize", new Matrix3( 1, 0, 0, // Latitude (no wrap around) - 0, 360, -180, // Longitude in [-180 … 180) range. + 0, 1, 180, // Longitude in [0 … 360] range. 0, 0, 1), MathTransforms.getMatrix(steps.get(2)), STRICT); /* - * The normalization is the inverse of above matrix (when source and target axes have the same span). - * But we expect the normalization matrix to be concatenated with the (3, 2, 5) scale operation. + * The normalization is the inverse of above matrix. + * But we expect the normalization matrix to be concatenated with the (3, 5) scale operation. */ assertMatrixEquals("normalize", new Matrix3( - 3, 0, 0, // 3 is a factor in MathTransforms.scale(…). - 0, 5./360, -(-180./360), // 5 is (idem). - 0, 0, 1), - MathTransforms.getMatrix(steps.get(0)), 1E-15); + 3, 0, 0, // 3 is a factor in MathTransforms.scale(…). + 0, 5, -180, // 5 is (idem). + 0, 0, 1), + MathTransforms.getMatrix(steps.get(0)), STRICT); + /* + * Test transforming some points. + */ + final double[] pts = { + 2, -100/5, + 6, -200/5, + 9, 200/5, + 3, 400/5}; + wt.transform(pts, 0, pts, 0, 4); + assertArrayEquals(new double[] { + 6, 260, + 18, 160, + 27, 200, + 9, 40}, pts, STRICT); } /** - * Tests wraparound on two axes. We expects two instances of {@link WraparoundTransform} without linear + * Tests wraparound on two axes. We expect two instances of {@link WraparoundTransform} without linear * transform between them. The absence of separation between the two {@link WraparoundTransform}s is an * indirect test of {@link WraparoundTransform#tryConcatenate(boolean, MathTransform, MathTransformFactory)}. * @@ -120,7 +136,7 @@ public final strictfp class WraparoundTransformTest extends TestCase { * Transform should be [scale & normalization] → [wraparound 1] → [wraparound 2] → [denormalization]. * At first an affine transform existed between the two [wraparound] operations, but that affine transform * should have been moved by `WraparoundTransform.tryConcatenate(…)` in order to combine them with initial - * [normalization} and final {denormalization]. + * [normalization] and final [denormalization]. */ final MathTransform wt = WraparoundTransform.forTargetCRS(op); final List steps = MathTransforms.getSteps(wt); @@ -128,24 +144,24 @@ public final strictfp class WraparoundTransformTest extends TestCase { assertEquals(0, ((WraparoundTransform) steps.get(1)).wraparoundDimension); assertEquals(2, ((WraparoundTransform) steps.get(2)).wraparoundDimension); /* - * Wraparound outputs are in [0 … 1) range (0 inclusive and 1 exclusive), so we expect a - * multiplication by the span of each axis for getting the final range. + * WraparoundTransform outputs are in [−180 … 180] range in longitude case, + * so we expect a 180° shift for getting results in the [0 … 360]° range. */ assertMatrixEquals("denormalize", new Matrix4( - 360, 0, 0, -180, // Longitude in [-180 … 180) range. - 0, 1, 0, 0, // Latitude (no wrap around) - 0, 0, 365, 1, // Day of year in [1 … 366) range. - 0, 0, 0, 1), + 1, 0, 0, 0, // Longitude in [-180 … 180] range. + 0, 1, 0, 0, // Latitude (no wrap around) + 0, 0, 1, 183.5, // Day of year in [1 … 366] range. + 0, 0, 0, 1), MathTransforms.getMatrix(steps.get(3)), STRICT); /* * The normalization is the inverse of above matrix (when source and target axes have the same span). * But we expect the normalization matrix to be concatenated with the (3, 2, 5) scale operation. */ assertMatrixEquals("normalize", new Matrix4( - 3./360, 0, 0, -(-180./360), // 3 is a factor in MathTransforms.scale(…). - 0, 2, 0, 0, // 2 is (idem). - 0, 0, 5./365, -(1./365), // 5 is (idem). - 0, 0, 0, 1), + 3, 0, 0, 0, // 3 is a factor in MathTransforms.scale(…). + 0, 2, 0, 0, // 2 is (idem). + 0, 0, 5, -183.5, // 5 is (idem). + 0, 0, 0, 1), MathTransforms.getMatrix(steps.get(0)), 1E-15); } } diff --git a/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java b/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java index fe049be..3f6be3f 100644 --- a/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java +++ b/core/sis-utility/src/main/java/org/apache/sis/measure/Longitude.java @@ -151,7 +151,7 @@ public final class Longitude extends Angle { public static double normalize(double λ) { /* * Following should be simplified as only one branch by javac since - * the values used in the 'if' statement are compile-time constants. + * the values used in the `if` statement are compile-time constants. * For verifying: javap -c org.apache.sis.measure.Longitude */ if (MIN_VALUE == -MAX_VALUE) {