This is an automated email from the ASF dual-hosted git repository. jsorel pushed a commit to branch feat/matrixutils in repository https://gitbox.apache.org/repos/asf/sis.git
commit d63bee034a26b9c21211e1af3068e671e4b0e935 Author: jsorel <[email protected]> AuthorDate: Mon May 23 11:02:54 2022 +0200 TileMatrix : add utility method for tiling scheme detection --- .../apache/sis/storage/tiling/TileMatrices.java | 188 +++++++++++++++++++++ .../sis/storage/tiling/TileMatricesTest.java | 152 +++++++++++++++++ .../apache/sis/test/suite/StorageTestSuite.java | 1 + 3 files changed, 341 insertions(+) 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 0000000000..a2cd272153 --- /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 0000000000..bdbd90390d --- /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 2392c0f1de..c24f31381c 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 @@ import org.junit.BeforeClass; 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,
