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 f1c09cd373da5f39ae22d8413b6a0181228a1e28 Author: Martin Desruisseaux <martin.desruisse...@geomatys.com> AuthorDate: Thu Sep 8 17:37:23 2022 +0200 Review the new `upsampling(int...)` methods with the following changes: In GridExtent: * Operate directly on high grid coordinates instead of (low + grid size). The operations on size in `subsample` were because of rounding problem. There is no such rounding with `upsample`. * Check for integer overflow in multiplications (ArithmeticException). * Opportunistic documentation fixes in various places (not only upsample) In particular the discussion about envelope was unclear in `subsample`. In GridGeometry: * Make the code safe to missing extent, gridToCRS, envelope, or CRS. * Matrix may be non-square (source and target dimensions not equal). * Avoid `IndexOutOfBoundsException` if `periods` array is truncated. * Use `DoubleDouble` arithmetic for reducing rounding errors. * Share the same `Envelope` instance. * Moved the method closer to related methods (subgrid). --- .../apache/sis/coverage/grid/GridDerivation.java | 2 +- .../org/apache/sis/coverage/grid/GridExtent.java | 41 +++--- .../org/apache/sis/coverage/grid/GridGeometry.java | 160 ++++++++++++--------- .../apache/sis/coverage/grid/GridExtentTest.java | 5 +- .../apache/sis/coverage/grid/GridGeometryTest.java | 65 ++++----- 5 files changed, 147 insertions(+), 126 deletions(-) diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java index 6a5ef345ef..6b82fbe6e7 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridDerivation.java @@ -807,7 +807,7 @@ public class GridDerivation { * They are not in units of cells of the size that we get after subsampling. * * @param indices the envelopes to intersect in units of {@link #base} grid coordinates. - * Shall contains at least one element. + * Shall contain at least one element. * @throws DisjointExtentException if the given envelope does not intersect the grid extent. * * @see #getBaseExtentExpanded(boolean) diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java index d0488546cb..ce9832b8fa 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridExtent.java @@ -91,6 +91,7 @@ import org.opengis.coverage.PointOutsideCoverageException; * * @author Martin Desruisseaux (IRD, Geomatys) * @author Alexis Manin (Geomatys) + * @author Johann Sorel (Geomatys) * @version 1.3 * @since 1.0 * @module @@ -1466,7 +1467,10 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable * of one cell may exist). * * <div class="note"><b>Note:</b> - * The envelope computed from a grid extent may become <em>larger</em> after subsampling, not smaller. + * If the "real world" envelope computed from grid extent needs to stay approximately the same, then the + * {@linkplain GridGeometry#getGridToCRS grid to CRS} transform needs to compensate the subsampling with + * a pre-multiplication of each grid coordinates by {@code periods}. + * However the envelope computed that way may become <em>larger</em> after subsampling, not smaller. * This effect can be understood intuitively if we consider that cells become larger after subsampling, * which implies that accurate representation of the same envelope may require fractional cells on some * grid borders.</div> @@ -1480,7 +1484,7 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable * If the array is longer, extraneous values are ignored. * * @param periods the subsampling. Length shall be equal to the number of dimension and all values shall be greater than zero. - * @return the subsampled extent, or {@code this} is subsampling results in the same extent. + * @return the subsampled extent, or {@code this} if subsampling results in the same extent. * @throws IllegalArgumentException if a period is not greater than zero. * * @see GridDerivation#subgrid(GridExtent, int...) @@ -1500,11 +1504,12 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable throw new ArithmeticException(Errors.format(Errors.Keys.IntegerOverflow_1, Long.SIZE)); } long r = Long.divideUnsigned(size, s); - if (r*s == size) r--; // Make inclusive if the division did not already rounded toward 0. + if (r*s == size) r--; // Make inclusive if the division did not already rounded toward 0. sub.coordinates[i] = low /= s; sub.coordinates[j] = low + r; } else if (s <= 0) { - throw new IllegalArgumentException(Errors.format(Errors.Keys.ValueNotGreaterThanZero_2, Strings.toIndexed("periods", i), s)); + throw new IllegalArgumentException(Errors.format( + Errors.Keys.ValueNotGreaterThanZero_2, Strings.toIndexed("periods", i), s)); } } return Arrays.equals(coordinates, sub.coordinates) ? this : sub; @@ -1512,13 +1517,9 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable /** * Creates a new grid extent upsampled by the given amount of cells along each grid dimensions. - * This method multiplies {@linkplain #getLow(int) low coordinates} and {@linkplain #getSize(int) grid sizes} + * This method multiplies {@linkplain #getLow(int) low} and {@linkplain #getHigh(int) high} coordinates * by the given periods. * - * <div class="note"><b>Note:</b> - * The envelope computed from a grid extent is preserved after upsampling. - * </div> - * * This method does not change the number of dimensions of the grid extent. * * <h4>Number of arguments</h4> @@ -1527,27 +1528,27 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable * If the array is longer, extraneous values are ignored. * * @param periods the upsampling. Length shall be equal to the number of dimension and all values shall be greater than zero. - * @return the upsampled extent, or {@code this} is upsampling results in the same extent. + * @return the upsampled extent, or {@code this} if upsampling results in the same extent. * @throws IllegalArgumentException if a period is not greater than zero. + * @throws ArithmeticException if the upsampled extent overflows the {@code long} capacity. + * + * @see GridGeometry#upsample(int...) + * @since 1.3 */ - public GridExtent upsample(int... periods) { + public GridExtent upsample(final int... periods) { ArgumentChecks.ensureNonNull("periods", periods); final int m = getDimension(); final int length = Math.min(m, periods.length); final GridExtent sub = new GridExtent(this); for (int i=0; i<length; i++) { - final long s = periods[i]; + final int s = periods[i]; if (s > 1) { final int j = i + m; - long low = coordinates[i]; - long size = coordinates[j] - low + 1; // Result is an unsigned number. - if (size == 0) { - throw new ArithmeticException(Errors.format(Errors.Keys.IntegerOverflow_1, Long.SIZE)); - } - sub.coordinates[i] = low * s; - sub.coordinates[j] = sub.coordinates[i] + size * s -1; + sub.coordinates[i] = Math.multiplyExact(coordinates[i], s); + sub.coordinates[j] = Math.addExact(Math.multiplyExact(coordinates[j], s), s-1); } else if (s <= 0) { - throw new IllegalArgumentException(Errors.format(Errors.Keys.ValueNotGreaterThanZero_2, Strings.toIndexed("periods", i), s)); + throw new IllegalArgumentException(Errors.format( + Errors.Keys.ValueNotGreaterThanZero_2, Strings.toIndexed("periods", i), s)); } } return Arrays.equals(coordinates, sub.coordinates) ? this : sub; diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java index 8ebb22fb70..1c0d5e8f5c 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java @@ -55,6 +55,7 @@ import org.apache.sis.referencing.operation.matrix.MatrixSIS; import org.apache.sis.referencing.operation.transform.MathTransforms; import org.apache.sis.referencing.operation.transform.LinearTransform; import org.apache.sis.referencing.operation.transform.PassThroughTransform; +import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix; import org.apache.sis.internal.referencing.DirectPositionView; import org.apache.sis.internal.referencing.TemporalAccessor; import org.apache.sis.internal.referencing.AxisDirections; @@ -132,6 +133,7 @@ import static org.apache.sis.referencing.CRS.findOperation; * The same instance can be shared by different {@link GridCoverage} instances. * * @author Martin Desruisseaux (IRD, Geomatys) + * @author Johann Sorel (Geomatys) * @version 1.3 * @since 1.0 * @module @@ -729,6 +731,18 @@ public class GridGeometry implements LenientComparable, Serializable { } } + /** + * Converts a "grid to CRS" transform from "cell corner" convention to "cell center" convention. + * This is a helper method for use of {@link #GridGeometry(GridExtent, MathTransform, MathTransform, + * ImmutableEnvelope, double[], long)} constructor. + * + * @param gridToCRS the transform to convert, or {@code null} if none. + * @return the converted transform, or {@code null} if the given transform was null. + */ + private static MathTransform cornerToCenter(final MathTransform gridToCRS) { + return PixelTranslation.translate(gridToCRS, PixelInCell.CELL_CORNER, PixelInCell.CELL_CENTER); + } + /** * Returns the number of dimensions of the <em>grid</em>. This is typically the same * than the number of {@linkplain #getEnvelope() envelope} dimensions or the number of @@ -853,7 +867,7 @@ public class GridGeometry implements LenientComparable, Serializable { * @return linear approximation of the conversion from grid coordinates to "real world" coordinates. * @throws IllegalArgumentException if the given {@code anchor} is not a known code list value. * @throws IncompleteGridGeometryException if this grid geometry has no transform, - * of if the transform is non-linear but this grid geometry has no extent. + * or if the transform is non-linear but this grid geometry has no extent. * @throws TransformException if an error occurred while computing the tangent. * * @since 1.3 @@ -1364,6 +1378,70 @@ public class GridGeometry implements LenientComparable, Serializable { return new GridDerivation(this); } + /** + * Creates a new grid geometry upsampled by the given amount of cells along each grid dimensions. + * This method multiplies {@linkplain GridExtent#getLow(int) low} and {@linkplain GridExtent#getHigh(int) high} + * coordinates by the given periods, then scales the {@link #getGridToCRS(PixelInCell) grid to CRS} transform + * for compensating the grid change. + * The grid geometry {@link #getEnvelope() envelope} is preserved after upsampling. + * + * <h4>Number of arguments</h4> + * The {@code periods} array length should be equal to the {@linkplain #getDimension() number of grid dimensions}. + * If the array is shorter, missing values default to 1 (i.e. samplings in unspecified dimensions are unchanged). + * If the array is longer, extraneous values are ignored. + * + * @param periods the upsampling. Length shall be equal to the number of dimension and all values shall be greater than zero. + * @return the upsampled grid geometry, or {@code this} is upsampling results in the same extent. + * @throws IllegalArgumentException if a period is not greater than zero. + * + * @see GridExtent#upsample(int...) + * @since 1.3 + */ + public GridGeometry upsample(final int... periods) { + GridExtent newExtent = extent; + if (newExtent != null) { + newExtent = newExtent.upsample(periods); + if (newExtent == extent) return this; // Unchanged. + } + boolean changed = false; // Additional check in case `extent` was null. + MathTransform newGridToCRS = null; + double[] newResolution = null; + if (cornerToCRS != null) { + final int tgtDim = cornerToCRS.getTargetDimensions(); + final int srcDim = cornerToCRS.getSourceDimensions(); + MatrixSIS matrix = Matrices.copy(MathTransforms.getMatrix(cornerToCRS)); + final boolean isNonLinear = (matrix == null); + if (isNonLinear) { + matrix = Matrices.create(tgtDim+1, srcDim+1, ExtendedPrecisionMatrix.IDENTITY); + } + /* + * By dividing the matrix elements directly, we avoid some numerical errors. + * This is because 100*0.1 is not exactly equal to 100/10. + * Each period value must be multiplied by a full column. + */ + newResolution = resolution.clone(); + final DoubleDouble div = new DoubleDouble(); + for (int i = Math.min(srcDim, periods.length); --i >= 0;) { + for (int j=0; j<tgtDim; j++) { + div.error = 0; + div.value = periods[i]; + if (div.value != 1) { + newResolution[i] /= div.value; + div.inverseDivide(DoubleDouble.castOrCopy(matrix.getNumber(j, i))); + matrix.setNumber(j, i, div); + changed = true; + } + } + } + newGridToCRS = MathTransforms.linear(matrix); + if (isNonLinear) { + newGridToCRS = MathTransforms.concatenate(newGridToCRS, cornerToCRS); + } + } + if (!changed) return this; + return new GridGeometry(newExtent, cornerToCenter(newGridToCRS), newGridToCRS, envelope, newResolution, nonLinears); + } + /** * Returns a grid geometry translated by the given amount of cells compared to this grid. * The returned grid has the same {@linkplain GridExtent#getSize(int) size} than this grid, @@ -1388,10 +1466,10 @@ public class GridGeometry implements LenientComparable, Serializable { */ public GridGeometry translate(final long... translation) { ArgumentChecks.ensureNonNull("translation", translation); - GridExtent te = extent; - if (te != null) { - te = te.translate(translation); - if (te == extent) return this; + GridExtent newExtent = extent; + if (newExtent != null) { + newExtent = newExtent.translate(translation); + if (newExtent == extent) return this; } MathTransform t1 = gridToCRS; MathTransform t2 = cornerToCRS; @@ -1407,7 +1485,7 @@ public class GridGeometry implements LenientComparable, Serializable { t1 = MathTransforms.concatenate(t, t1); t2 = MathTransforms.concatenate(t, t2); } - return new GridGeometry(te, t1, t2, envelope, resolution, nonLinears); + return new GridGeometry(newExtent, t1, t2, envelope, resolution, nonLinears); } /** @@ -1418,27 +1496,27 @@ public class GridGeometry implements LenientComparable, Serializable { * <p>The given extent is taken verbatim; this method does no clipping. * The given extent does not need to intersect the extent of this grid geometry.</p> * - * @param extent extent of the grid geometry to return. + * @param newExtent extent of the grid geometry to return. * @return grid geometry with the given extent. May be {@code this} if there is no change. * @throws TransformException if the geospatial envelope can not be recomputed with the new grid extent. * * @since 1.3 */ - public GridGeometry relocate(final GridExtent extent) throws TransformException { - ArgumentChecks.ensureNonNull("size", extent); - if (extent.equals(this.extent)) { + public GridGeometry relocate(final GridExtent newExtent) throws TransformException { + ArgumentChecks.ensureNonNull("newExtent", newExtent); + if (newExtent.equals(extent)) { return this; } - ensureDimensionMatches(getDimension(), extent); + ensureDimensionMatches(getDimension(), newExtent); final ImmutableEnvelope relocated; if (cornerToCRS != null) { - final GeneralEnvelope env = extent.toEnvelope(cornerToCRS, gridToCRS, null); + final GeneralEnvelope env = newExtent.toEnvelope(cornerToCRS, gridToCRS, null); env.setCoordinateReferenceSystem(getCoordinateReferenceSystem(envelope)); relocated = new ImmutableEnvelope(env); } else { relocated = envelope; // Either null or contains only the CRS. } - return new GridGeometry(extent, gridToCRS, cornerToCRS, relocated, resolution, nonLinears); + return new GridGeometry(newExtent, gridToCRS, cornerToCRS, relocated, resolution, nonLinears); } /** @@ -1470,62 +1548,6 @@ public class GridGeometry implements LenientComparable, Serializable { return this; } - /** - * Creates a new grid geometry upsampling the GridExtent by the given amount of cells along each grid dimensions. - * This method multiplies {@linkplain GridExtent#getLow(int) low coordinates} and {@linkplain GridExtent#getSize(int) grid sizes} - * by the given periods. - * - * <div class="note"><b>Note:</b> - * The envelope of the new grid geometry is preserved after upsampling. - * </div> - * - * This method does not change the number of dimensions of the grid geometry. - * - * <h4>Number of arguments</h4> - * The {@code periods} array length should be equal to the {@linkplain #getDimension() number of dimensions}. - * If the array is shorter, missing values default to 1 (i.e. samplings in unspecified dimensions are unchanged). - * If the array is longer, extraneous values are ignored. - * - * @param periods the upsampling. Length shall be equal to the number of dimension and all values shall be greater than zero. - * @return the upsampled grid geometry, or {@code this} is upsampling results in the same extent. - * @throws IllegalArgumentException if a period is not greater than zero. - * - * @see GridExtent#upsample(int...) - */ - public GridGeometry upsample(int... periods) { - - final GridExtent extent = getExtent(); - final GridExtent upExtent = extent.upsample(periods); - if (upExtent == extent) { - //unchanged - return this; - } - final int dimension = upExtent.getDimension(); - - final MathTransform gridToCrs = getGridToCRS(PixelInCell.CELL_CORNER); - final MathTransform upGridToCrs; - if (gridToCrs instanceof LinearTransform) { - /* - By dividing the matrix elements directly we avoid some numeric errors. - */ - final LinearTransform lnt = (LinearTransform) gridToCrs; - final MatrixSIS matrix = Matrices.copy(lnt.getMatrix()); - for (int i = 0; i < dimension; i++) { - for (int k = 0; k < dimension; k++) { - matrix.setElement(k, i, matrix.getElement(k, i) / periods[i]); - } - } - upGridToCrs = MathTransforms.linear(matrix); - } else { - final double[] scaling = new double[dimension]; - for (int i = 0; i < dimension; i++) { - scaling[i] = 1.0 / periods[i]; - } - upGridToCrs = MathTransforms.concatenate(MathTransforms.scale(scaling), gridToCrs); - } - return new GridGeometry(upExtent, PixelInCell.CELL_CORNER, upGridToCrs, getCoordinateReferenceSystem()); - } - /** * Creates a one-, two- or three-dimensional coordinate reference system for cell indices in the grid. * This method returns a CRS which is derived from the "real world" CRS or a subset of it. diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java index 430d76647d..4a7b022bd5 100644 --- a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java +++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridExtentTest.java @@ -47,6 +47,7 @@ import static org.apache.sis.test.ReferencingAssert.*; * * @author Martin Desruisseaux (IRD, Geomatys) * @author Alexis Manin (Geomatys) + * @author Johann Sorel (Geomatys) * @version 1.3 * @since 1.0 * @module @@ -72,7 +73,7 @@ public final strictfp class GridExtentTest extends TestCase { } /** - * Tests the {@link GridExtent#subsample(int...)}. + * Tests {@link GridExtent#subsample(int...)}. */ @Test public void testSubsample() { @@ -84,7 +85,7 @@ public final strictfp class GridExtentTest extends TestCase { } /** - * Tests the {@link GridExtent#upsample(int...)}. + * Tests {@link GridExtent#upsample(int...)}. */ @Test public void testUpsample() { diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java index 02e8b6a86b..4654589a6b 100644 --- a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java +++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/GridGeometryTest.java @@ -44,6 +44,7 @@ import static org.apache.sis.test.ReferencingAssert.*; * Tests the {@link GridGeometry} implementation. * * @author Martin Desruisseaux (IRD, Geomatys) + * @author Johann Sorel (Geomatys) * @version 1.3 * @since 1.0 * @module @@ -467,6 +468,36 @@ public final strictfp class GridGeometryTest extends TestCase { new double[] { 111319.49, -331876.53}), envelope, 0.01); } + /** + * Tests {@link GridGeometry#upsample(int...)}. + */ + @Test + public void testUpsample() { + final GridGeometry grid; + { // Source grid + final GridExtent extent = new GridExtent(null, new long[] {10,-8}, new long[] {100, 50}, false); + final Matrix3 mat = new Matrix3( + 1, 0, 10, + 0, -2, 50, + 0, 0, 1); + final MathTransform gridToCRS = MathTransforms.linear(mat); + grid = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, HardCodedCRS.CARTESIAN_2D); + } + final GridGeometry upsampled = grid.upsample(4, 4); + final GridGeometry expected; + { // Expected grid + GridExtent extent = new GridExtent(null, new long[] {40,-32}, new long[] {400, 200}, false); + final Matrix3 mat = new Matrix3( + 0.25, 0, 9.625, + 0, -0.5, 50.750, + 0, 0, 1); + final MathTransform gridToCRS = MathTransforms.linear(mat); + expected = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, HardCodedCRS.CARTESIAN_2D); + } + assertSame(grid.getEnvelope(), upsampled.getEnvelope()); + assertEquals(expected, upsampled); + } + /** * Tests {@link GridGeometry#translate(long...)}. */ @@ -561,40 +592,6 @@ public final strictfp class GridGeometryTest extends TestCase { verifyGridToCRS(grid); } - /** - * Tests {@link GridGeometry#upsample(int...)}. - */ - @Test - public void testUpsample() { - - final GridGeometry grid; - { //grid - final GridExtent extent = new GridExtent(null, new long[]{10,-8}, new long[]{100, 50}, false); - final Matrix3 mat = new Matrix3( - 1, 0, 10, - 0, -2, 50, - 0, 0, 1); - final MathTransform gridToCRS = MathTransforms.linear(mat); - grid = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, HardCodedCRS.CARTESIAN_2D); - } - - final GridGeometry surSampling = grid.upsample(4, 4); - - final GridGeometry expected; - { //expected grid - GridExtent extent = new GridExtent(null, new long[]{40,-32}, new long[]{400, 200}, false); - final Matrix3 mat = new Matrix3( - 0.25, 0, 9.625, - 0, -0.5, 50.750, - 0, 0, 1); - final MathTransform gridToCRS = MathTransforms.linear(mat); - expected = new GridGeometry(extent, PixelInCell.CELL_CENTER, gridToCRS, HardCodedCRS.CARTESIAN_2D); - } - - assertEquals(grid.getEnvelope(), surSampling.getEnvelope()); - assertEquals(expected, surSampling); - } - /** * Tests {@link GridGeometry#reduce(int...)} with a {@code gridToCRS} transform having a constant value * in one dimension. This method tests indirectly {@link SliceGeometry#findTargetDimensions(MathTransform,