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 9f8a4020152877c9a3dc82a162080a6c4d28f894 Author: Martin Desruisseaux <martin.desruisse...@geomatys.com> AuthorDate: Fri Apr 14 20:27:09 2023 +0200 Relax the restriction that all inputs to `BandAggregateGridCoverage` have the same grid geometry. With this commit, they are allowed to have different translations but not yet more complex changes. --- .../coverage/grid/BandAggregateGridCoverage.java | 9 +- .../coverage/grid/CoordinateOperationFinder.java | 5 +- .../apache/sis/coverage/grid/DefaultEvaluator.java | 2 +- .../sis/coverage/grid/GridCoverageProcessor.java | 2 + .../org/apache/sis/coverage/grid/GridExtent.java | 22 +- .../org/apache/sis/image/MultiSourceLayout.java | 17 +- .../sis/internal/coverage/CommonDomainFinder.java | 372 +++++++++++++++++++++ .../sis/internal/coverage/MultiSourceArgument.java | 39 ++- .../org/apache/sis/internal/feature/Resources.java | 5 + .../sis/internal/feature/Resources.properties | 1 + .../sis/internal/feature/Resources_fr.properties | 1 + .../grid/BandAggregateGridCoverageTest.java | 137 ++++++++ .../apache/sis/test/suite/FeatureTestSuite.java | 1 + .../apache/sis/storage/aggregate/GridSlice.java | 64 +--- .../sis/storage/aggregate/GridSliceLocator.java | 3 +- .../sis/storage/aggregate/GroupByTransform.java | 4 +- 16 files changed, 592 insertions(+), 92 deletions(-) diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java index 9e671073a9..11a98fe7a4 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/BandAggregateGridCoverage.java @@ -64,6 +64,12 @@ final class BandAggregateGridCoverage extends GridCoverage { */ private final int numBands; + /** + * Translations in units of grid cells to apply on an extent of this grid coverage + * for getting a "grid to CRS" transform compatible with each source. + */ + private final long[][] gridTranslations; + /** * Index of a sources having the same "grid to CRS" than this grid coverage, or -1 if none. */ @@ -94,6 +100,7 @@ final class BandAggregateGridCoverage extends GridCoverage { this.sources = aggregate.sources(); this.bandsPerSource = aggregate.bandsPerSource(true); this.numBands = aggregate.numBands(); + this.gridTranslations = aggregate.gridTranslations(); this.sourceOfGridToCRS = aggregate.sourceOfGridToCRS(); this.processor = processor; this.dataType = sources[0].getBandType(); @@ -138,7 +145,7 @@ final class BandAggregateGridCoverage extends GridCoverage { } final RenderedImage[] images = new RenderedImage[sources.length]; for (int i=0; i<images.length; i++) { - images[i] = sources[i].render(sliceExtent); + images[i] = sources[i].render(sliceExtent.translate(gridTranslations[i])); } return processor.aggregateBands(images, bandsPerSource); } diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java index e5f653f6d5..ee1ab98896 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java @@ -233,10 +233,11 @@ final class CoordinateOperationFinder implements Supplier<double[]> { } /** - * Verifies whether the presence of a CRS considered mandatory, unless the CRS of opposite grid - * is also missing. + * Verifies the presence of a CRS considered mandatory, + * unless the CRS of the opposite grid is also missing. * * @param rs {@code true} is source CRS is mandatory, {@code false} if target CRS is mandatory. + * @throws IncompleteGridGeometryException if a mandatory CRS is missing. */ final void verifyPresenceOfCRS(final boolean rs) { if ((rs ? target : source).isDefined(GridGeometry.CRS)) { diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java index 66a7572630..65467914e6 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DefaultEvaluator.java @@ -639,7 +639,7 @@ class DefaultEvaluator implements GridCoverage.Evaluator { } } // Modify fields only after everything else succeeded. - position = new FractionalGridCoordinates.Position(crsToGrid.getTargetDimensions()); + position = new FractionalGridCoordinates.Position(crsToGrid.getTargetDimensions()); inputCRS = crs; inputToGrid = crsToGrid; } diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java index 9394a7c603..d04f92c113 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverageProcessor.java @@ -773,6 +773,8 @@ public class GridCoverageProcessor implements Cloneable { * <li>All coverages shall use the same data type in their rendered image.</li> * </ul> * + * Some of those restrictions may be relaxed in future versions. + * * @param sources coverages whose bands shall be aggregated, in order. At least one coverage must be provided. * @param bandsPerSource bands to use for each source coverage, in order. May contain {@code null} elements. * @return the aggregated coverage, or one of the sources if it can be used directly. 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 de45ba8a24..99da3df96e 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 @@ -1724,18 +1724,18 @@ public class GridExtent implements GridEnvelope, LenientComparable, Serializable ArgumentChecks.ensureNonNull("translation", translation); final int m = getDimension(); final int length = Math.min(m, translation.length); - if (!isZero(translation, length)) { - final GridExtent translated = new GridExtent(this); - final long[] c = translated.coordinates; - for (int i=0; i < length; i++) { - final int j = i + m; - final long t = translation[i]; - c[i] = Math.addExact(c[i], t); - c[j] = Math.addExact(c[j], t); - } - return translated; + if (isZero(translation, length)) { + return this; + } + final GridExtent translated = new GridExtent(this); + final long[] c = translated.coordinates; + for (int i=0; i < length; i++) { + final int j = i + m; + final long t = translation[i]; + c[i] = Math.addExact(c[i], t); + c[j] = Math.addExact(c[j], t); } - return this; + return translated; } /** diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java b/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java index a49f466d96..d4eae377bc 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java +++ b/core/sis-feature/src/main/java/org/apache/sis/image/MultiSourceLayout.java @@ -34,6 +34,7 @@ import org.apache.sis.internal.coverage.j2d.ImageLayout; import org.apache.sis.internal.coverage.j2d.ImageUtilities; import org.apache.sis.internal.coverage.j2d.ColorModelFactory; import org.apache.sis.internal.coverage.MultiSourceArgument; +import org.apache.sis.internal.coverage.CommonDomainFinder; import org.apache.sis.coverage.grid.DisjointExtentException; @@ -174,12 +175,18 @@ final class MultiSourceLayout extends ImageLayout { } /* * Get the domain of the combined image to create. - * TODO: current implementation computes the intersection of all sources. - * But a future version should allow users to specify if they want intersection, - * union or strict mode instead. A "strict" mode would prevent the combination of - * images using different domains (i.e. raise an error if domains are not the same). + * Current implementation computes the intersection of all sources. */ - ImageUtilities.clipBounds(source, domain); + if (CommonDomainFinder.INTERSECTION) { + ImageUtilities.clipBounds(source, domain); + } else if (!domain.equals(ImageUtilities.getBounds(source))) { + /* + * TODO: a future version should allow users to specify if they want intersection, + * union or strict mode instead. A "strict" mode would prevent the combination of + * images using different domains (i.e. raise an error if domains are not the same). + */ + throw new IllegalArgumentException(); + } if (domain.isEmpty()) { throw new DisjointExtentException(Resources.format(Resources.Keys.SourceImagesDoNotIntersect)); } diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/CommonDomainFinder.java b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/CommonDomainFinder.java new file mode 100644 index 0000000000..f6d28d8698 --- /dev/null +++ b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/CommonDomainFinder.java @@ -0,0 +1,372 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.sis.internal.coverage; + +import java.util.Map; +import java.util.LinkedHashMap; +import java.util.NoSuchElementException; +import org.opengis.util.FactoryException; +import org.opengis.geometry.MismatchedDimensionException; +import org.opengis.referencing.datum.PixelInCell; +import org.opengis.referencing.operation.Matrix; +import org.opengis.referencing.operation.MathTransform; +import org.opengis.referencing.operation.TransformException; +import org.opengis.referencing.crs.CoordinateReferenceSystem; +import org.opengis.referencing.operation.NoninvertibleTransformException; +import org.apache.sis.referencing.CRS; +import org.apache.sis.referencing.operation.transform.MathTransforms; +import org.apache.sis.coverage.grid.GridExtent; +import org.apache.sis.coverage.grid.GridGeometry; +import org.apache.sis.coverage.grid.IllegalGridGeometryException; +import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix; +import org.apache.sis.internal.feature.Resources; +import org.apache.sis.internal.util.Numerics; +import org.apache.sis.util.Numbers; + + +/** + * Helper class for building a combined domain from a list of grid geometries. + * After construction, one of the following methods shall be invoked exactly once. + * + * <ul> + * <li>{@link #setFromGridAligned(GridGeometry...)}</li> + * </ul> + * + * Then, the result can be obtained by the given methods: + * + * <ul> + * <li>{@link #result()}</li> + * <li>{@link #gridTranslations()}</li> + * <li>{@link #sourceOfGridToCRS()}</li> + * </ul> + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.4 + * @since 1.4 + */ +public final class CommonDomainFinder { + /** + * Whether to compute intersection (true) or union (false) of all grid coverages. + * This is not yet configurable in current version. + * + * <p>We use this constant for tracking the code to update when we will want to provide an option for using + * union or strict equality instead (the latter would be a mode that fails if all extents are not identical). + * Note that in the case of unions, it would be possible to specify coverages with no intersection. + * Whether we should accept that or raise an exception is still an open question.</p> + */ + public static final boolean INTERSECTION = true; + + /** + * The grid geometry taken as the reference for computing the translations of all grid geometry items. + * Values of {@link #gridTranslations} and {@link #itemTranslations} are relative to that reference. + * At first this is an arbitrary grid geometry, for example the first encountered one. + * After {@code CommonDomainFinder} completed its task, this is updated to the final result. + * + * @see #result() + */ + private GridGeometry reference; + + /** + * Coordinate reference system of the reference, or {@code null} if not yet known. + * It will be set to the CRS of the first grid geometry where the CRS is defined. + */ + private CoordinateReferenceSystem crs; + + /** + * The inverse of the "grid to CRS" transform of the grid geometry taken as a reference. + */ + private MathTransform crsToGrid; + + /** + * The convention to use for fetching the "grid to CRS" transforms. + */ + private final PixelInCell anchor; + + /** + * The combined extent, as the union or intersection of all grid extents. + */ + private GridExtent extent; + + /** + * The translation in units of grid cells from the {@linkplain #reference} grid geometry + * to the grid geometry in the key. + */ + private final Map<GridGeometry,long[]> gridTranslations; + + /** + * Translations in units of grid cells from the {@linkplain #reference} grid geometry to each item. + * For each index <var>i</var>, {@code itemTranslations[i]} is a value from {@link #gridTranslations} + * map and may be reused at more than one index <var>i</var>. + * + * @see #gridTranslations() + */ + private long[][] itemTranslations; + + /** + * If one of the grid geometries has the same "grid to CRS" than the common grid geometry, the index. + * Otherwise -1. + */ + private int sourceOfGridToCRS; + + /** + * Creates a new common domain finder. + * + * @param anchor the convention to use for fetching the "grid to CRS" transforms. + */ + CommonDomainFinder(final PixelInCell anchor) { + this.anchor = anchor; + gridTranslations = new LinkedHashMap<>(); + } + + /** + * Computes a common grid geometry from the given items. + * All items shall share be aligned on the same grid. + * Items may be translated relative to each other, + * but the translations shall be an integer number of grid cells. + * + * <h4>Coordinate reference system</h4> + * If the CRS of a grid geometry is undefined, it is assumed the same than other grid geometries. + * + * @param items the grid geometries for which to compute a common grid geometry. + * @throws IllegalGridGeometryException if the specified item is not compatible with the reference grid geometry. + */ + final void setFromGridAligned(final GridGeometry... items) { + itemTranslations = new long[items.length][]; + for (int i=0; i<items.length; i++) { + itemTranslations[i] = gridTranslations.computeIfAbsent(items[i], this::itemToCommon); + } + /* + * Change the reference grid geometry for matching more closely the desired grid extent. + * If one item has exactly the desired grid extent, use it. Otherwise search for an item + * having the same origin. This criterion is arbitrary and may change in future version. + */ + GridGeometry fallback = null; + for (final Map.Entry<GridGeometry,long[]> entry : gridTranslations.entrySet()) { + final GridGeometry item = entry.getKey(); + final GridExtent actual = item.getExtent(); + final GridExtent expected = extent.translate(entry.getValue()); + if (actual.equals(expected)) { + setGridToCRS(items, item); // Must be before `reference` assignation. + reference = item; + return; + } + if (fallback == null && extent.getLow().equals(item.getExtent().getLow())) { + fallback = item; + } + } + if (fallback == null) { + fallback = reference; + } + setGridToCRS(items, fallback); + try { + reference = fallback.relocate(extent); + } catch (TransformException e) { + throw new IllegalGridGeometryException(Resources.format(Resources.Keys.IncompatibleGridGeometries), e); + } + } + + /** + * Given a grid geometry with the desired "grid to CRS", saves its index in {@link #sourceOfGridToCRS}. + * This method updates all previously computed translations for making them relative to the new reference. + * + * Note: updating values of the {@link #gridTranslations} map indirectly update all + * {@link #itemTranslations} array elements. + */ + private void setGridToCRS(final GridGeometry[] items, final GridGeometry item) { + sourceOfGridToCRS = indexOf(items, item); + if (item == reference) { // Quick check for a common case. + return; + } + final long[] oldReference = itemTranslations[indexOf(items, reference)]; + final long[] newReference = itemTranslations[sourceOfGridToCRS]; + final long[] change = new long[newReference.length]; + for (int i=0; i < newReference.length; i++) { + change[i] = Math.subtractExact(newReference[i], oldReference[i]); + } + for (final long[] offset : gridTranslations.values()) { + for (int i=0; i < offset.length; i++) { + offset[i] = Math.subtractExact(offset[i], change[i]); + } + } + } + + /** + * Returns the index of the given grid geometry. + * This method is invoked when the instance should always exist in the array. + */ + private static int indexOf(final GridGeometry[] items, final GridGeometry item) { + for (int i=0; i<items.length; i++) { + if (items[i] == item) { + return i; + } + } + throw new NoSuchElementException(); // Should never happen. + } + + /** + * Returns the common grid geometry computed from all specified items. + * + * @return a grid geometry which is the union or intersection of all specified items. + */ + final GridGeometry result() { + if (crs != null && !reference.isDefined(GridGeometry.CRS)) { + reference = new GridGeometry(extent, anchor, reference.getGridToCRS(anchor), crs); + } + return reference; + } + + /** + * Returns the translations (in units of grid cells) from the common grid geometry to all items. + * The items are the arguments given to {@link #setFromGridAligned(GridGeometry...)}, in order. + * The common grid geometry is the value returned by {@link #result()}. + * + * <p>The returned array should not be modified because it is not cloned.</p> + * + * @return translation from the common grid geometry to all items. This array is not cloned. + */ + @SuppressWarnings("ReturnOfCollectionOrArrayField") + final long[][] gridTranslations() { + return itemTranslations; + } + + /** + * If one items has the same "grid to CRS" than the common grid geometry, returns its index. + * + * @return index of an items having the desired "grid to CRS", or -1 if none. + */ + final int sourceOfGridToCRS() { + return sourceOfGridToCRS; + } + + /** + * Computes the translation in units of grid cells from the common grid geometry to the given item. + * This method also opportunistically computes the union or intersection of all grid extents. + * + * @param item the grid geometry for which to get a translation from the common grid geometry. + * @return translation in unis of grid cells. Note that the caller may reuse this array for many grid geometries. + * @throws IllegalGridGeometryException if the specified item is not compatible with the reference grid geometry. + */ + private long[] itemToCommon(final GridGeometry item) { + /* + * Compute the change ourselves instead of invoking `GridGeometry.createTransformTo(…)` + * because we do not want wraparound handling when we search for a simple translation. + */ + MathTransform change = item.getGridToCRS(anchor); + try { + if (crsToGrid == null) { + crsToGrid = change.inverse(); + reference = item; + } + if (item.isDefined(GridGeometry.CRS)) { + final CoordinateReferenceSystem src = item.getCoordinateReferenceSystem(); + if (crs == null) { + crs = src; + } else { + /* + * Ask for a change of CRS without specifying an area of interest (AOI) on the assumption + * that if the transformation is only a translation, the AOI would not make a difference. + * It save not only the AOI computation cost, but also make easier for `findOperation(…)` + * to use its cache. + */ + change = MathTransforms.concatenate(change, CRS.findOperation(src, crs, null).getMathTransform()); + } + } + change = MathTransforms.concatenate(change, crsToGrid); + } catch (FactoryException | NoninvertibleTransformException | MismatchedDimensionException e) { + throw new IllegalGridGeometryException(Resources.format(Resources.Keys.IncompatibleGridGeometries), e); + } + final long[] offset = integerTranslation(MathTransforms.getMatrix(change), null); + if (offset == null) { + throw new IllegalGridGeometryException(Resources.format(Resources.Keys.IncompatibleGridGeometries)); + } + /* + * The grid geometry has been accepted as valid. Now compute the combined extent, + * taking offset in account. At this point this is an offset TO the common grid geometry. + * It will be converted to an offset FROM the common grid geometry after the extent update. + */ + final GridExtent e = item.getExtent().translate(offset); + if (extent == null) { + extent = e; + } else if (INTERSECTION) { + extent = extent.intersect(e); + } else if (!extent.equals(e)) { + throw new IllegalGridGeometryException(); + } + for (int i=0; i<offset.length; i++) { + offset[i] = Math.negateExact(offset[i]); + } + return offset; + } + + /** + * If the given matrix is the identity matrix except for translation terms, returns the translation. + * Translation terms must be integer values and will be stored in the given {@code offset} array. + * If the matrix is not an integer translation, this method return {@code null} without modifying + * the given {@code offset} array. + * + * @param change conversion between two grid geometries, or {@code null}. + * @param offset where to store translation terms if the change is an integer translation, or {@code null}. + * @return the translation terms, or {@code null} if the given matrix does not met the conditions. + * + * @see org.apache.sis.referencing.operation.matrix.Matrices#isTranslation(Matrix) + */ + public static long[] integerTranslation(final Matrix change, long[] offset) { + if (change == null) { + return null; + } + final int numRows = change.getNumRow(); + final int numCols = change.getNumCol(); + for (int j=0; j<numRows; j++) { + for (int i=0; i<numCols; i++) { + double tolerance = Numerics.COMPARISON_THRESHOLD; + double e = change.getElement(j, i); + if (i == j) { + e--; + } else if (i == numCols - 1) { + final double a = Math.abs(e); + if (a > 1) { + tolerance = Math.min(tolerance*a, 0.125); + } + e -= Math.rint(e); + } + if (!(Math.abs(e) <= tolerance)) { // Use `!` for catching NaN. + return null; + } + } + } + /* + * Store the translation terms after we have determined that they are integers. + * It must be an "all or nothing" operation (unless an exception is thrown). + */ + if (offset == null) { + offset = new long[numRows - 1]; + } + final int i = numCols - 1; + if (change instanceof ExtendedPrecisionMatrix) { + final var epm = (ExtendedPrecisionMatrix) change; + for (int j=0; j<offset.length; j++) { + final Number e = epm.getElementOrNull(j, i); + offset[j] = (e != null) ? Numbers.round(e) : 0; + } + } else { + for (int j=0; j<offset.length; j++) { + offset[j] = Math.round(change.getElement(j, i)); + } + } + return offset; + } +} diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java index 520ab3479a..3f3c1fe100 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java +++ b/core/sis-feature/src/main/java/org/apache/sis/internal/coverage/MultiSourceArgument.java @@ -25,6 +25,7 @@ import java.util.Objects; import java.util.function.Consumer; import java.util.function.Function; import java.util.function.ToIntFunction; +import org.opengis.referencing.datum.PixelInCell; import org.apache.sis.coverage.SampleDimension; import org.apache.sis.coverage.grid.GridGeometry; import org.apache.sis.coverage.grid.IllegalGridGeometryException; @@ -32,11 +33,10 @@ import org.apache.sis.internal.util.Numerics; import org.apache.sis.util.resources.Errors; import org.apache.sis.util.ArgumentChecks; import org.apache.sis.util.ArraysExt; -import org.apache.sis.util.ComparisonMode; /** - * Helper class for building a list of sample dimensions from aggregated sources. + * Helper class for building a combined domain or range from aggregated sources. * This helper class is shared for aggregation operations on different sources: * rendered images, grid coverages and resources. * @@ -112,6 +112,12 @@ public final class MultiSourceArgument<S> { */ private List<SampleDimension> ranges; + /** + * Translations in units of grid cells to apply for obtaining a grid geometry + * compatible with the "grid to CRS" transform of a source. + */ + private long[][] gridTranslations; + /** * Index of a source having the same "grid to CRS" transform than the grid geometry * returned by {@link #domain(Function)}. If there is none, then this value is -1. @@ -510,19 +516,28 @@ next: for (int i=0; i<sources.length; i++) { // `sources.length` may * @param getter the method to invoke for getting grid geometry from a source. * @return intersection of all grid geometries. * @throws IllegalGridGeometryException if a grid geometry is not compatible with the others. - * - * @todo Current implementation requires that all grid geometry are equal. We need to relax that. */ public GridGeometry domain(final Function<S, GridGeometry> getter) { checkValidationState(true); - GridGeometry intersection = getter.apply(sources[0]); - for (int i=1; i < validatedSourceCount; i++) { - if (!intersection.equals(getter.apply(sources[i]), ComparisonMode.IGNORE_METADATA)) { - throw new IllegalGridGeometryException("Not yet supported on coverages with different grid geometries."); - } - } - sourceOfGridToCRS = 0; // TODO: to be computed when different grid geometries will be allowed. Prefer widest extent. - return intersection; + final var finder = new CommonDomainFinder(PixelInCell.CELL_CORNER); + finder.setFromGridAligned(Arrays.stream(sources).map(getter).toArray(GridGeometry[]::new)); + sourceOfGridToCRS = finder.sourceOfGridToCRS(); + gridTranslations = finder.gridTranslations(); + return finder.result(); + } + + /** + * Returns the translations in units of grid cells to apply for obtaining a grid geometry + * compatible with the "grid to CRS" transform of a source. + * + * <p>The returned array should not be modified because it is not cloned.</p> + * + * @return translations from the common grid geometry to all items. This array is not cloned. + */ + @SuppressWarnings("ReturnOfCollectionOrArrayField") + public long[][] gridTranslations() { + checkValidationState(true); + return gridTranslations; } /** diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java index caf2a3ec8e..1c89d38a7d 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java +++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.java @@ -256,6 +256,11 @@ public final class Resources extends IndexedResourceBundle { */ public static final short IncompatibleColorModel = 82; + /** + * At least two coverages have mutually incompatible grid geometries. + */ + public static final short IncompatibleGridGeometries = 87; + /** * The ({0}, {1}) tile has an unexpected size, number of bands or sample layout. */ diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties index 0a1c6c2018..853ab78e0a 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties +++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources.properties @@ -58,6 +58,7 @@ ImageAllowsTransparency = Image allows transparency. ImageHasAlphaChannel = Image has alpha channel. ImageIsOpaque = Image is opaque. IncompatibleColorModel = Color model is incompatible with sample model. +IncompatibleGridGeometries = At least two coverages have mutually incompatible grid geometries. IncompatibleTile_2 = The ({0}, {1}) tile has an unexpected size, number of bands or sample layout. InsufficientBufferCapacity_3 = Data buffer capacity is insufficient for a grid of {0} cells \u00d7 {1} bands. Missing {2} elements. IterationIsFinished = Iteration is finished. diff --git a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties index 356dc1cecf..0ae9c2ce3b 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties +++ b/core/sis-feature/src/main/java/org/apache/sis/internal/feature/Resources_fr.properties @@ -63,6 +63,7 @@ ImageAllowsTransparency = L\u2019image permet la transparence. ImageHasAlphaChannel = L\u2019image a un canal alpha. ImageIsOpaque = L\u2019image est opaque. IncompatibleColorModel = Le mod\u00e8le de couleurs est incompatible. +IncompatibleGridGeometries = Au moins deux couvertures de donn\u00e9es ont des g\u00e9om\u00e9tries de grilles mutuellement incompatibles. IncompatibleTile_2 = La tuile ({0}, {1}) a une taille, un nombre de bandes ou une disposition des valeurs inattendu. InsufficientBufferCapacity_3 = La capacit\u00e9 du buffer est insuffisante pour une grille de {0} cellules \u00d7 {1} bandes. Il lui manque {2} \u00e9l\u00e9ments. IterationIsFinished = L\u2019it\u00e9ration est termin\u00e9e. diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/BandAggregateGridCoverageTest.java b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/BandAggregateGridCoverageTest.java new file mode 100644 index 0000000000..ade6c2caa5 --- /dev/null +++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/BandAggregateGridCoverageTest.java @@ -0,0 +1,137 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.sis.coverage.grid; + +import java.util.Arrays; +import java.awt.image.Raster; +import java.awt.image.DataBufferInt; +import org.opengis.referencing.datum.PixelInCell; +import org.apache.sis.coverage.SampleDimension; +import org.apache.sis.referencing.crs.HardCodedCRS; +import org.apache.sis.referencing.operation.transform.MathTransforms; +import org.apache.sis.test.TestCase; +import org.junit.Test; + +import static org.junit.Assert.*; + + +/** + * Tests {@link BandAggregateGridCoverage}. + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.4 + * @since 1.4 + */ +public final class BandAggregateGridCoverageTest extends TestCase { + /** + * Width and height of images created for tests. + */ + private static final int WIDTH = 3, HEIGHT = 2; + + /** + * The processor to use for creating the aggregated grid coverage. + */ + private final GridCoverageProcessor processor; + + /** + * Creates a new test case. + */ + public BandAggregateGridCoverageTest() { + processor = new GridCoverageProcessor(); + } + + /** + * Tests aggregation with two coverages having the same grid geometry. + */ + @Test + public void testSameGridGeometry() { + final GridCoverage c1 = createCoverage(-2, 4, 3, -1, 100, 200); + final GridCoverage c2 = createCoverage(-2, 4, 3, -1, 300); + final GridCoverage cr = processor.aggregateRanges(c1, c2); + assertEquals(c1.getGridGeometry(), cr.getGridGeometry()); + assertEquals(c2.getGridGeometry(), cr.getGridGeometry()); + assertEquals(3, cr.getSampleDimensions().size()); + assertPixelsEqual(cr, 100, 200, 300, 101, 201, 301, 102, 202, 302, + 103, 203, 303, 104, 204, 304, 105, 205, 305); + } + + /** + * Tests aggregation with two coverages having a translation in their grid extents. + * Their "grid to CRS" transforms are the same, which implies that the "real world" + * coordinates are different. The intersection is not equal to any source extent. + */ + @Test + public void testNoCommonGridGeometry() { + final GridCoverage c1 = createCoverage(-2, 4, 3, -1, 100, 200); + final GridCoverage c2 = createCoverage(-1, 3, 3, -1, 300); + final GridCoverage cr = processor.aggregateRanges(c1, c2); + final GridExtent extent = cr.getGridGeometry().getExtent(); + assertNotEquals(c1.getGridGeometry().getExtent(), extent); + assertNotEquals(c2.getGridGeometry().getExtent(), extent); + assertEquals(extent(-1, 4, 1, 5), extent); + assertEquals(3, cr.getSampleDimensions().size()); + assertPixelsEqual(cr, 101, 201, 303, 102, 202, 304); + } + + /** + * Returns a two-dimensional grid extents with the given bounding box. + * The maximal coordinates are exclusive. + */ + private static GridExtent extent(final int minX, final int minY, final int maxX, final int maxY) { + return new GridExtent(null, new long[] {minX, minY}, new long[] {maxX, maxY}, false); + } + + /** + * Creates a new grid coverage with bands starting with the given sample values. + * The length of the {@code bandValues} array is the number of bands to create. + * In a given band <var>b</var>, all pixels have the {@code bandValues[b]}. + * + * @param minX minimal <var>x</var> coordinate value of the grid extent. + * @param minY minimal <var>y</var> coordinate value of the grid extent. + * @param translateX <var>x</var> component of the "grid to CRS" translation. + * @param translateY <var>y</var> component of the "grid to CRS" translation. + * @param bandValues sample values for the first pixel. + * @return a coverage with an image where all pixels have the specified sample values. + */ + private static GridCoverage createCoverage(final int minX, final int minY, + final int translateX, final int translateY, final int... bandValues) + { + final int numBands = bandValues.length; + final var samples = new SampleDimension[numBands]; + final var builder = new SampleDimension.Builder(); + for (int i=0; i<numBands; i++) { + samples[i] = builder.setName("Band with value " + bandValues[i]).build(); + } + final int[] data = new int[WIDTH * HEIGHT * numBands]; + for (int i=0; i<data.length; i++) { + data[i] = bandValues[i % numBands] + i / numBands; + } + final var values = new DataBufferInt(data, data.length); + final var domain = new GridGeometry(extent(minX, minY, minX+WIDTH, minY+HEIGHT), + PixelInCell.CELL_CORNER, MathTransforms.translation(translateX, translateY), HardCodedCRS.WGS84); + return new BufferedGridCoverage(domain, Arrays.asList(samples), values); + } + + /** + * Asserts that all pixels from the given coverage have the expected values. + */ + private static void assertPixelsEqual(final GridCoverage cr, final int... expected) { + final Raster r = cr.render(null).getData(); + final int[] data = r.getPixels(r.getMinX(), r.getMinY(), r.getWidth(), r.getHeight(), (int[]) null); + assertArrayEquals(expected, data); + } +} diff --git a/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java b/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java index 3d09e913e5..4f233d0744 100644 --- a/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java +++ b/core/sis-feature/src/test/java/org/apache/sis/test/suite/FeatureTestSuite.java @@ -124,6 +124,7 @@ import org.junit.runners.Suite; org.apache.sis.coverage.grid.TranslatedGridCoverageTest.class, org.apache.sis.coverage.grid.ResampledGridCoverageTest.class, org.apache.sis.coverage.grid.DimensionalityReductionTest.class, + org.apache.sis.coverage.grid.BandAggregateGridCoverageTest.class, // Index and processing org.apache.sis.index.tree.PointTreeNodeTest.class, diff --git a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java index cc8bff034f..502cf6878a 100644 --- a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java +++ b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSlice.java @@ -29,10 +29,9 @@ import org.apache.sis.storage.DataStoreException; import org.apache.sis.coverage.grid.GridCoverage; import org.apache.sis.coverage.grid.GridGeometry; import org.apache.sis.coverage.grid.GridExtent; +import org.apache.sis.internal.coverage.CommonDomainFinder; import org.apache.sis.internal.storage.MemoryGridResource; -import org.apache.sis.internal.util.Numerics; import org.apache.sis.internal.util.Strings; -import org.apache.sis.util.Numbers; /** @@ -48,6 +47,11 @@ import org.apache.sis.util.Numbers; * @since 1.3 */ final class GridSlice { + /** + * The "pixel in cell" value used for all "grid to CRS" computations. + */ + static final PixelInCell CELL_ANCHOR = PixelInCell.CELL_CORNER; + /** * The resource associated to this slice. */ @@ -87,57 +91,6 @@ final class GridSlice { offset = new long[geometry.getDimension()]; } - /** - * Sets the {@link #offset} terms to the values of the translation columns of the given matrix. - * This method shall be invoked if and only if {@link #isIntegerTranslation(Matrix)} returned {@code true}. - * - * @param groupToSlice conversion from source coordinates of {@link GroupByTransform#gridToCRS} - * to grid coordinates of {@link #geometry}. - * @throws ArithmeticException if a translation term is NaN or overflows {@code long} integer capacity. - * - * @see #getOffset(Map) - */ - private void setOffset(final MatrixSIS groupToSlice) { - final int i = groupToSlice.getNumCol() - 1; - for (int j=0; j<offset.length; j++) { - offset[j] = Numbers.round(groupToSlice.getNumber(j, i)); - } - } - - /** - * Returns {@code true} if the given matrix is the identity matrix except for translation terms. - * Translation terms must be integer values. - * - * @param groupToSlice conversion from {@link GroupByTransform#gridToCRS} - * source coordinates to {@link #gridToCRS} source coordinates. - * @return whether the matrix is identity, ignoring integer translation. - * - * @see org.apache.sis.referencing.operation.matrix.Matrices#isTranslation(Matrix) - */ - private static boolean isIntegerTranslation(final Matrix groupToSlice) { - final int numRows = groupToSlice.getNumRow(); - final int numCols = groupToSlice.getNumCol(); - for (int j=0; j<numRows; j++) { - for (int i=0; i<numCols; i++) { - double tolerance = Numerics.COMPARISON_THRESHOLD; - double e = groupToSlice.getElement(j, i); - if (i == j) { - e--; - } else if (i == numCols - 1) { - final double a = Math.abs(e); - if (a > 1) { - tolerance = Math.min(tolerance*a, 0.125); - } - e -= Math.rint(e); - } - if (!(Math.abs(e) <= tolerance)) { // Use `!` for catching NaN. - return false; - } - } - } - return true; - } - /** * Returns the group of objects associated to the CRS and "grid to CRS" transform. * The CRS comparisons ignore metadata and transform comparisons ignore integer translations. @@ -152,14 +105,13 @@ final class GridSlice { final GroupByTransform getList(final List<GroupByCRS<GroupByTransform>> groups, final MergeStrategy strategy) throws NoninvertibleTransformException { - final MathTransform gridToCRS = geometry.getGridToCRS(PixelInCell.CELL_CORNER); + final MathTransform gridToCRS = geometry.getGridToCRS(CELL_ANCHOR); final MathTransform crsToGrid = gridToCRS.inverse(); final List<GroupByTransform> transforms = GroupByCRS.getOrAdd(groups, geometry).members; synchronized (transforms) { for (final GroupByTransform c : transforms) { final Matrix groupToSlice = c.linearTransform(crsToGrid); - if (groupToSlice != null && isIntegerTranslation(groupToSlice)) { - setOffset(MatrixSIS.castOrCopy(groupToSlice)); + if (CommonDomainFinder.integerTranslation(groupToSlice, offset) != null) { c.strategy = strategy; return c; } diff --git a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java index 6cecc715f3..137809f741 100644 --- a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java +++ b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GridSliceLocator.java @@ -22,7 +22,6 @@ import java.util.Arrays; import java.util.HashMap; import java.util.function.Function; import org.opengis.util.InternationalString; -import org.opengis.referencing.datum.PixelInCell; import org.opengis.metadata.spatial.DimensionNameType; import org.apache.sis.coverage.grid.GridExtent; import org.apache.sis.coverage.grid.GridGeometry; @@ -137,7 +136,7 @@ final class GridSliceLocator { return base; } extent = new GridExtent(axes, low, high, true); - return new GridGeometry(extent, PixelInCell.CELL_CORNER, base.getGridToCRS(PixelInCell.CELL_CORNER), + return new GridGeometry(extent, GridSlice.CELL_ANCHOR, base.getGridToCRS(GridSlice.CELL_ANCHOR), base.isDefined(GridGeometry.CRS) ? base.getCoordinateReferenceSystem() : null); } diff --git a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java index 540483160e..d8884efe15 100644 --- a/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java +++ b/storage/sis-storage/src/main/java/org/apache/sis/storage/aggregate/GroupByTransform.java @@ -51,7 +51,7 @@ final class GroupByTransform extends Group<GridSlice> { private final GridGeometry geometry; /** - * Value or {@code geometry.getGridToCRS(PixelInCell.CELL_CORNER)}. + * Value of {@code geometry.getGridToCRS(GridSlice.CELL_ANCHOR)}. * All {@linkplain #members} of this group use this transform, * possibly with integer differences in translation terms. */ @@ -68,7 +68,7 @@ final class GroupByTransform extends Group<GridSlice> { * Creates a new group of objects associated to the given transform. * * @param geometry geometry of the grid coverage or resource. - * @param gridToCRS value or {@code geometry.getGridToCRS(PixelInCell.CELL_CORNER)}. + * @param gridToCRS value of {@code geometry.getGridToCRS(GridSlice.CELL_ANCHOR)}. * @param strategy algorithm to apply when more than one grid coverage can be found at the same grid index. */ GroupByTransform(final GridGeometry geometry, final MathTransform gridToCRS, final MergeStrategy strategy) {