TileMatrix : add utility method for tiling scheme detection
diff --git a/storage/sis-storage/src/main/java/org/apache/sis/storage/tiling/TileMatrices.java b/storage/sis-storage/src/main/java/org/apache/sis/storage/tiling/TileMatrices.java
new file mode 100644
index 0000000..a2cd272
--- /dev/null
+++ b/storage/sis-storage/src/main/java/org/apache/sis/storage/tiling/TileMatrices.java
@@ -0,0 +1,188 @@
+/*
+ * (C) 2022, Geomatys
+ */
+package org.apache.sis.storage.tiling;
+
+import java.util.AbstractMap;
+import java.util.Arrays;
+import java.util.HashMap;
+import java.util.Map;
+import java.util.Map.Entry;
+import java.util.Optional;
+import java.util.stream.LongStream;
+import org.apache.sis.coverage.grid.GridClippingMode;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.coverage.grid.GridRoundingMode;
+import org.apache.sis.referencing.operation.transform.LinearTransform;
+import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.util.Static;
+import org.apache.sis.util.Utilities;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.datum.PixelInCell;
+import org.opengis.referencing.operation.MathTransform;
+import org.opengis.referencing.operation.Matrix;
+
+/**
+ * TileMatrix utilities.
+ *
+ * @author Johann Sorel (Geomatys)
+ * @version 1.3
+ * @since   1.3
+ */
+public final class TileMatrices extends Static {
+
+    private TileMatrices(){}
+
+    /**
+     * Compute tiling scheme of a list of tile geometries.
+     *
+     * @param grids tile grid geometries
+     * @return tiling scheme and map of each tile indices or empty if any tile do not fit in the scheme.
+     */
+    public static Optional<Entry<GridGeometry, Map<GridGeometry,long[]>>> toTilingScheme(GridGeometry... grids) {
+
+        final GridGeometry first = grids[0];
+        GridGeometry tilingScheme;
+
+        {   // Create a first tiling scheme, pick the first grid as a tile reference
+
+            if (!first.isDefined(GridGeometry.EXTENT | GridGeometry.CRS | GridGeometry.GRID_TO_CRS)) {
+                //we don't have enough informations to compute the tiling scheme
+                return Optional.empty();
+            }
+            final MathTransform firstGridToCRS = first.getGridToCRS(PixelInCell.CELL_CENTER);
+            if (!(firstGridToCRS instanceof LinearTransform) || !((LinearTransform) firstGridToCRS).isAffine()) {
+                //detection works only for affine transforms
+                return Optional.empty();
+            }
+
+            //create a first tiling scheme
+            final GridGeometry forceLowerToZero = forceLowerToZero(first);
+            final int[] subsampling = LongStream.of(forceLowerToZero.getExtent().getHigh().getCoordinateValues())
+                    .mapToInt(Math::toIntExact)
+                    .map((int operand) -> operand+1) //high values are inclusive
+                    .toArray();
+            tilingScheme = forceLowerToZero.derive().subgrid(null, subsampling).build();
+        }
+
+        //compute all tile indices
+        final Map<GridGeometry,long[]> tiles = new HashMap<>();
+        final int dimension = tilingScheme.getExtent().getDimension();
+        tiles.put(first, tilingScheme.getExtent().getLow().getCoordinateValues());
+        final long[] min = new long[dimension];
+        final long[] max = new long[dimension];
+        for (int i = 1; i < grids.length; i++) {
+            final Optional<long[]> indice = getTileIndices(first, tilingScheme, grids[i]);
+            if (!indice.isPresent()) return Optional.empty();
+            long[] r = indice.get();
+            tiles.put(grids[i], r);
+
+            //keep track of min/max range
+            for (int k = 0; k < dimension; k++) {
+                min[k] = Math.min(min[k], r[k]);
+                max[k] = Math.max(max[k], r[k]);
+            }
+        }
+
+        //rebuild the tiling scheme extent to contain all tiles
+        tilingScheme = new GridGeometry(
+                new GridExtent(null, min, max, true),
+                PixelInCell.CELL_CENTER,
+                tilingScheme.getGridToCRS(PixelInCell.CELL_CENTER),
+                tilingScheme.getCoordinateReferenceSystem());
+        tilingScheme = forceLowerToZero(tilingScheme);
+
+        //offset all indices
+        for (Entry<GridGeometry,long[]> entry : tiles.entrySet()) {
+            final long[] indices = entry.getValue();
+            for (int i = 0; i < dimension; i++) {
+                indices[i] -= min[i];
+            }
+        }
+
+        return Optional.of(new AbstractMap.SimpleImmutableEntry<>(tilingScheme, tiles));
+    }
+
+    /**
+     * Find tile indice in given tiling scheme.
+     *
+     * @param referenceTile a valid tile used a reference.
+     * @param tilingScheme the tiling scheme geometry.
+     * @param tileGrid searched tile grid geometry.
+     * @return tile index or empty if tile do not fit in the scheme.
+     */
+    private static Optional<long[]> getTileIndices(GridGeometry referenceTile, GridGeometry tilingScheme, GridGeometry tileGrid) {
+        if (!tileGrid.isDefined(GridGeometry.EXTENT | GridGeometry.CRS | GridGeometry.GRID_TO_CRS)) {
+            //we don't have enough informations to compute the tile indices
+            return Optional.empty();
+        }
+        if (!Utilities.equalsIgnoreMetadata(referenceTile.getCoordinateReferenceSystem(), tileGrid.getCoordinateReferenceSystem())) {
+            //tile candidate has different CRS
+            return Optional.empty();
+        }
+        final MathTransform gridToCRS = tileGrid.getGridToCRS(PixelInCell.CELL_CENTER);
+        if (!(gridToCRS instanceof LinearTransform) || !((LinearTransform) gridToCRS).isAffine()) {
+            //indice computation works only for affine transforms
+            return Optional.empty();
+        }
+
+        //matrices must differ only by the last column (translation)
+        final LinearTransform firstLinear = (LinearTransform) referenceTile.getGridToCRS(PixelInCell.CELL_CENTER);
+        final Matrix matrix1 = firstLinear.getMatrix();
+        final LinearTransform linear2 = (LinearTransform) gridToCRS;
+        final Matrix matrix2 = linear2.getMatrix();
+        for (int x = 0, xn = matrix1.getNumCol() - 1, yn = matrix1.getNumRow(); x < xn; x++) {
+            for (int y = 0; y < yn; y++) {
+                if (matrix1.getElement(y, x) != matrix2.getElement(y, x)) {
+                    return Optional.empty();
+                }
+            }
+        }
+
+        //tiles must have the same extent size
+        final GridExtent referenceExtent = referenceTile.getExtent();
+        final GridExtent candidateExtent = tileGrid.getExtent();
+        for (int i = 0, n = referenceExtent.getDimension(); i < n; i++) {
+            if (referenceExtent.getSize(i) != candidateExtent.getSize(i)) {
+                return Optional.empty();
+            }
+        }
+
+        //compute the tile indice
+        final GridExtent intersection = tilingScheme.derive().clipping(GridClippingMode.NONE).rounding(GridRoundingMode.ENCLOSING).subgrid(tileGrid).getIntersection();
+        final long[] low = intersection.getLow().getCoordinateValues();
+        final long[] high = intersection.getHigh().getCoordinateValues();
+
+        //if tile overlaps several indices then it's not part of the tiling scheme
+        if (!Arrays.equals(low, high)) {
+            return Optional.empty();
+        }
+
+        return Optional.of(low);
+    }
+
+    /**
+     * Shift lower extent to zero.
+     */
+    private static GridGeometry forceLowerToZero(final GridGeometry gg) {
+        final GridExtent extent = gg.getExtent();
+        if (!extent.startsAtZero()) {
+            CoordinateReferenceSystem crs = null;
+            if (gg.isDefined(GridGeometry.CRS)) crs = gg.getCoordinateReferenceSystem();
+            final int dimension = extent.getDimension();
+            final double[] vector = new double[dimension];
+            final long[] high = new long[dimension];
+            for (int i = 0; i < dimension; i++) {
+                final long low = extent.getLow(i);
+                high[i] = extent.getHigh(i) - low;
+                vector[i] = low;
+            }
+            MathTransform gridToCRS = gg.getGridToCRS(PixelInCell.CELL_CENTER);
+            gridToCRS = MathTransforms.concatenate(MathTransforms.translation(vector), gridToCRS);
+            return new GridGeometry(new GridExtent(null, null, high, true), PixelInCell.CELL_CENTER, gridToCRS, crs);
+        }
+        return gg;
+    }
+
+}
diff --git a/storage/sis-storage/src/test/java/org/apache/sis/storage/tiling/TileMatricesTest.java b/storage/sis-storage/src/test/java/org/apache/sis/storage/tiling/TileMatricesTest.java
new file mode 100644
index 0000000..bdbd903
--- /dev/null
+++ b/storage/sis-storage/src/test/java/org/apache/sis/storage/tiling/TileMatricesTest.java
@@ -0,0 +1,152 @@
+/*
+ * (C) 2022, Geomatys
+ */
+package org.apache.sis.storage.tiling;
+
+import java.util.Map;
+import java.util.Optional;
+import org.apache.sis.coverage.grid.GridExtent;
+import org.apache.sis.coverage.grid.GridGeometry;
+import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
+import org.apache.sis.referencing.CommonCRS;
+import org.apache.sis.test.TestCase;
+import org.junit.Test;
+import org.opengis.referencing.crs.CoordinateReferenceSystem;
+import org.opengis.referencing.datum.PixelInCell;
+
+import static org.junit.Assert.*;
+
+/**
+ * Tests for {@link TileMatrices}.
+ *
+ * @author  Johann Sorel (Geomatys)
+ * @version 1.3
+ * @since   1.3
+ *
+ * @author Johann Sorel (Geomatys)
+ */
+public final class TileMatricesTest extends TestCase {
+
+    public TileMatricesTest() {
+    }
+
+    /**
+     * Test tiling scheme detection from tiles.
+     */
+    @Test
+    public void testToTilingScheme() {
+
+        final long tileWidth = 256;
+        final long tileHeight = 256;
+        final CoordinateReferenceSystem crs = CommonCRS.WGS84.normalizedGeographic();
+        final double scaleX = 0.1;
+        final double scaleY = 0.1;
+        final double translateX = 0.0;
+        final double translateY = 0.0;
+
+        //tile {0,0}
+        final GridGeometry tile00 = new GridGeometry(new GridExtent(tileWidth, tileHeight), PixelInCell.CELL_CORNER,
+                new AffineTransform2D(scaleX, 0, 0, scaleY, translateX, translateY), crs);
+        //tile {1,0}
+        final GridGeometry tile10 = new GridGeometry(new GridExtent(tileWidth, tileHeight), PixelInCell.CELL_CORNER,
+            new AffineTransform2D(scaleX, 0, 0, scaleY, translateX+tileWidth*scaleX, translateY), crs);
+        //tile {-4,-2} using gridToCRS
+        final GridGeometry tilem42 = new GridGeometry(new GridExtent(tileWidth, tileHeight), PixelInCell.CELL_CORNER,
+            new AffineTransform2D(scaleX, 0, 0, scaleY, translateX + -4 * tileWidth * scaleX, translateY + -2 * tileHeight * scaleY), crs);
+        //tile {-4,-2} using gridExtent
+        final GridGeometry tilem42bis = new GridGeometry(new GridExtent(null, new long[]{-4*tileWidth,-2*tileHeight}, new long[]{-3*tileWidth,-1*tileHeight}, false), PixelInCell.CELL_CORNER,
+            new AffineTransform2D(scaleX, 0, 0, scaleY, translateX, translateY), crs);
+
+        { // test single tile
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00);
+            assertTrue(result.isPresent());
+            final GridGeometry tilingScheme = result.get().getKey();
+            final Map<GridGeometry, long[]> indices = result.get().getValue();
+
+            final GridGeometry expected = new GridGeometry(new GridExtent(null, null, new long[]{1, 1}, false), PixelInCell.CELL_CORNER,
+                    new AffineTransform2D(tileWidth*scaleX, 0, 0, tileHeight*scaleY, 0, 0), crs);
+            //assertEquals(expected, tilingScheme);
+            assertTrue(expected.equals(tilingScheme));
+            assertArrayEquals(new long[]{0,0}, indices.get(tile00));
+        }
+
+
+        { // test with a valid second tile
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00, tile10);
+            assertTrue(result.isPresent());
+            final GridGeometry tilingScheme = result.get().getKey();
+            final Map<GridGeometry, long[]> indices = result.get().getValue();
+
+            final GridGeometry expected = new GridGeometry(new GridExtent(null, null, new long[]{2, 1}, false), PixelInCell.CELL_CORNER,
+                    new AffineTransform2D(tileWidth*scaleX, 0, 0, tileHeight*scaleY, 0, 0), crs);
+            //assertEquals(expected, tilingScheme);
+            assertTrue(expected.equals(tilingScheme));
+            assertArrayEquals(new long[]{0,0}, indices.get(tile00));
+            assertArrayEquals(new long[]{1,0}, indices.get(tile10));
+        }
+
+        { // test with a valid third tile
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00, tile10, tilem42);
+            assertTrue(result.isPresent());
+            final GridGeometry tilingScheme = result.get().getKey();
+            final Map<GridGeometry, long[]> indices = result.get().getValue();
+
+            final GridGeometry expected = new GridGeometry(new GridExtent(null, null, new long[]{6, 3}, false), PixelInCell.CELL_CORNER,
+                    new AffineTransform2D(tileWidth*scaleX, 0, 0, tileHeight*scaleY, -4 * tileWidth * scaleX, -2 * tileHeight * scaleY), crs);
+            //assertEquals(expected, tilingScheme);
+            assertTrue(expected.equals(tilingScheme));
+            assertArrayEquals(new long[]{4,2}, indices.get(tile00));
+            assertArrayEquals(new long[]{5,2}, indices.get(tile10));
+            assertArrayEquals(new long[]{0,0}, indices.get(tilem42));
+        }
+
+        { // test with a valid third tile, same as previous test but tile translation is expressed in the grid extent
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00, tile10, tilem42bis);
+            assertTrue(result.isPresent());
+            final GridGeometry tilingScheme = result.get().getKey();
+            final Map<GridGeometry, long[]> indices = result.get().getValue();
+
+            final GridGeometry expected = new GridGeometry(new GridExtent(null, null, new long[]{6, 3}, false), PixelInCell.CELL_CORNER,
+                    new AffineTransform2D(tileWidth*scaleX, 0, 0, tileHeight*scaleY, -4 * tileWidth * scaleX, -2 * tileHeight * scaleY), crs);
+            //assertEquals(expected, tilingScheme);
+            assertTrue(expected.equals(tilingScheme));
+            assertArrayEquals(new long[]{4,2}, indices.get(tile00));
+            assertArrayEquals(new long[]{5,2}, indices.get(tile10));
+            assertArrayEquals(new long[]{0,0}, indices.get(tilem42bis));
+        }
+
+        { // test a tile with a different crs
+            final GridGeometry badTile = new GridGeometry(new GridExtent(tileWidth, tileHeight), PixelInCell.CELL_CORNER,
+                new AffineTransform2D(scaleX, 0, 0, scaleY, translateX+tileWidth*scaleX, translateY), CommonCRS.WGS84.geographic());
+
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00, badTile);
+            assertFalse(result.isPresent());
+        }
+
+        { // test a tile with a different size
+            final GridGeometry badTile = new GridGeometry(new GridExtent(tileWidth-1, tileHeight), PixelInCell.CELL_CORNER,
+                new AffineTransform2D(scaleX, 0, 0, scaleY, translateX+tileWidth*scaleX, translateY), crs);
+
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00, badTile);
+            assertFalse(result.isPresent());
+        }
+
+        { // test a tile with a different scale
+            final GridGeometry badTile = new GridGeometry(new GridExtent(tileWidth, tileHeight), PixelInCell.CELL_CORNER,
+                new AffineTransform2D(scaleX, 0, 0, scaleY+0.001, translateX, translateY), crs);
+
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00, badTile);
+            assertFalse(result.isPresent());
+        }
+
+        { // test a tile with a none aligned translation
+            final GridGeometry badTile = new GridGeometry(new GridExtent(tileWidth, tileHeight), PixelInCell.CELL_CORNER,
+                new AffineTransform2D(scaleX, 0, 0, scaleY, translateX+0.0001, translateY), crs);
+
+            final Optional<Map.Entry<GridGeometry, Map<GridGeometry, long[]>>> result = TileMatrices.toTilingScheme(tile00, badTile);
+            assertFalse(result.isPresent());
+        }
+
+    }
+
+}
diff --git a/storage/sis-storage/src/test/java/org/apache/sis/test/suite/StorageTestSuite.java b/storage/sis-storage/src/test/java/org/apache/sis/test/suite/StorageTestSuite.java
index 2392c0f..c24f313 100644
--- a/storage/sis-storage/src/test/java/org/apache/sis/test/suite/StorageTestSuite.java
+++ b/storage/sis-storage/src/test/java/org/apache/sis/test/suite/StorageTestSuite.java
@@ -50,6 +50,7 @@
     org.apache.sis.storage.event.StoreListenersTest.class,
     org.apache.sis.storage.CoverageQueryTest.class,
     org.apache.sis.storage.FeatureQueryTest.class,
+    org.apache.sis.storage.tiling.TileMatricesTest.class,
     org.apache.sis.internal.storage.xml.MimeTypeDetectorTest.class,
     org.apache.sis.internal.storage.xml.StoreProviderTest.class,
     org.apache.sis.internal.storage.xml.StoreTest.class,