diff --git a/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/reader/Localization.java b/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/reader/Localization.java index 78cbc7e196..f0cf26bd49 100644 --- a/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/reader/Localization.java +++ b/endorsed/src/org.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/reader/Localization.java @@ -25,6 +25,8 @@ import org.opengis.util.FactoryException; import org.opengis.referencing.operation.MathTransform; import org.opengis.referencing.operation.TransformException; +import org.apache.sis.geometry.GeneralEnvelope; +import org.apache.sis.referencing.operation.matrix.Matrix3; import org.apache.sis.referencing.operation.transform.MathTransforms; import org.apache.sis.referencing.operation.transform.LinearTransform; import org.apache.sis.referencing.operation.builder.LocalizationGridBuilder; @@ -110,6 +112,22 @@ private static MathTransform localizationGrid(final Vector modelTiePoints, final } catch (ArithmeticException | FactoryException e) { /* * May happen when the model tie points are not distributed on a regular grid. + * The tie points may nevertheless be on a rectilinear grid, i.e. every combination + * of the distinct x and y pixel coordinates may be present exactly once, with only + * the spacing between those coordinates being unsuitable for the above inference. + * This is the case of ICEYE images, where the tie points are at k × (size-1) / (n-1) + * pixels, sometimes rounded to integers. The greatest common divisor of those + * coordinates is then much smaller than the actual step, either because the step is + * fractional or because the rounding makes one step differ from the others by one + * pixel. Such grids are handled without splitting them, by using the ranks of the + * distinct coordinates as grid indices. + */ + final MathTransform rectilinear = rectilinearGrid(modelTiePoints, x, y, addTo); + if (rectilinear != null) { + return rectilinear; + } + /* + * Otherwise the tie points are really irregular. * For example, Sentinel 1 images may have tie points spaced by 1320 pixels on the X axis, * except the very last point which is only 1302 pixels after the previous one. We try to * handle such grids by splitting them in two parts: one grid for the columns where points @@ -124,6 +142,10 @@ private static MathTransform localizationGrid(final Vector modelTiePoints, final * │ 2 │ 3 │ * └──────────────────┴───┘ * splitX + * + * If the irregular spacing is on a single axis, then the threshold of the other axis is NaN, + * the comparisons against it are always false and only two of the four parts receive points. + * The empty parts are skipped. */ final Set uniques = new HashSet<>(100); final double splitX = threshold(x, uniques); @@ -182,6 +204,7 @@ private static MathTransform localizationGrid(final Vector modelTiePoints, final MathTransform global = null; final Map specialization = new LinkedHashMap<>(4); for (int i=0; ix + * and y pixel coordinates is present exactly once, and if those coordinates are + * evenly spaced up to a rounding to integers. The latter condition is verified by + * {@link #isUniformAfterRounding(double[], double)}. + * + *

Contrarily to the {@code localizationGrid(…)} fallback, this method does not split the + * tie points: the grid indices are the ranks of the distinct pixel coordinates, and the linear + * relationship between pixel coordinates and ranks is applied before the localization grid. + * Consequently the transform has no discontinuity and honors all tie points.

+ * + * @param modelTiePoints the model tie points read from GeoTIFF file. + * @param x the x pixel coordinates of the tie points. + * @param y the y pixel coordinates of the tie points. + * @param addTo if non-null, add the transform result to this map. + * @return the "grid to CRS" transform, or {@code null} if the tie points are not on a rectilinear grid. + */ + private static MathTransform rectilinearGrid(final Vector modelTiePoints, final Vector x, final Vector y, + final Map addTo) throws FactoryException, TransformException + { + final int size = modelTiePoints.size(); + final double[] ux = distinctSorted(x); + final double[] uy = distinctSorted(y); + final int nx = ux.length; + final int ny = uy.length; + if (nx < 2 || ny < 2 || ((long) nx) * ny != size / RECORD_LENGTH) { + return null; // Not a complete rectilinear grid. + } + final double sx = (ux[nx-1] - ux[0]) / (nx - 1); + final double sy = (uy[ny-1] - uy[0]) / (ny - 1); + if (!isUniformAfterRounding(ux, sx) || !isUniformAfterRounding(uy, sy)) { + return null; // Spacing is irregular for a real reason. + } + final LocalizationGridBuilder grid = new LocalizationGridBuilder(nx, ny); + for (int i=0; ix or y vector of tie points pixel coordinates. + * @return the distinct values, in increasing order. + */ + private static double[] distinctSorted(final Vector values) { + final int n = values.size(); + final Set uniques = new HashSet<>(100); + for (int i=0; i= 1;) { + if (!(Math.abs(values[i] - (values[0] + i*step)) < 1)) { + return false; // Use `!` for catching NaN. + } + } + return true; + } + /** * Finds the value at which the increment in localization grid seems to change. * This is used when not all tie points in a GeoTIFF images are distributed on diff --git a/endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/reader/LocalizationTest.java b/endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/reader/LocalizationTest.java new file mode 100644 index 0000000000..66f407f9dd --- /dev/null +++ b/endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/reader/LocalizationTest.java @@ -0,0 +1,264 @@ +/* + * 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.storage.geotiff.reader; + +import org.opengis.referencing.operation.MathTransform; +import org.apache.sis.math.Vector; + +// Test dependencies +import org.junit.jupiter.api.Test; +import static org.junit.jupiter.api.Assertions.*; +import org.apache.sis.test.TestCase; + + +/** + * Tests the construction of a localization grid from GeoTIFF tie points. + * The tie point spacings tested here are the spacings observed in real products: + * Sentinel 1 images, and ICEYE images in Ground Range Detected, Single Look Complex + * and ScanSAR flavors. + * + *

All tests build the tie points from a slightly non-linear model, then verify that the + * resulting transform maps each tie point to the model coordinates declared in the same record. + * The model has to be non-linear, otherwise every candidate transform would reproduce the tie + * points exactly and the tests would only verify the absence of exception.

+ * + * @author Jonatas Fischer + */ +public final class LocalizationTest extends TestCase { + /** + * Tolerance threshold, in degrees, when comparing tie point coordinates. + * This is about one metre, while the tie points of the grids tested here + * are kilometres apart. + */ + private static final double TOLERANCE = 1E-5; + + /** + * Number of tie points on each axis of most grids tested here. + * This is the number of tie points in ICEYE products. + */ + private static final int GRID_SIZE = 10; + + /** + * Creates a new test case. + */ + public LocalizationTest() { + } + + /** + * Returns the pixel coordinates of {@code count} evenly spaced tie points, the first one + * on the first pixel and the last one on the given last pixel. The coordinates are exact, + * i.e. they are not necessarily integers. + * + * @param count number of tie points. + * @param last pixel coordinate of the last tie point, usually the image size minus 1. + * @return pixel coordinates of the tie points. + */ + private static double[] evenSpacing(final int count, final double last) { + final double[] coordinates = new double[count]; + for (int i=0; ix axis. + * @param rows pixel coordinates of the tie points along the y axis. + * @return the (I,J,K,X,Y,Z) records of the tie points. + */ + private static Vector tiePoints(final double[] columns, final double[] rows) { + final double[] records = new double[columns.length * rows.length * Localization.RECORD_LENGTH]; + int p = 0; + for (final double y : rows) { + for (final double x : columns) { + records[p++] = x; + records[p++] = y; + records[p++] = 0; + records[p++] = -66 + x*1E-6 + y*3E-8 + (x*y)*2E-13; + records[p++] = 45 - y*1E-6 + x*5E-8 - (x*x)*1E-13; + records[p++] = 0; + } + } + return Vector.create(records, false); + } + + /** + * Builds the localization grid for the given tie points, then verifies that the resulting + * transform maps the pixel coordinates of each tie point to the model coordinates declared + * in the same record. + * + * @param columns pixel coordinates of the tie points along the x axis. + * @param rows pixel coordinates of the tie points along the y axis. + * @throws Exception if the transform cannot be created or used. + */ + private static void verify(final double[] columns, final double[] rows) throws Exception { + final Vector tiePoints = tiePoints(columns, rows); + final MathTransform gridToCRS = Localization.nonLinear(tiePoints); + assertNotNull(gridToCRS); + final double[] source = new double[2]; + final double[] target = new double[2]; + for (int i=0; iThis is the case of ICEYE ScanSAR images of 19250 × 19510 pixels, which have 39 × 40 + * tie points spaced by 506.55 and 500.23 pixels respectively.

+ * + * @throws Exception if the transform cannot be created or used. + */ + @Test + public void testGridWithFractionalSpacing() throws Exception { + verify(evenSpacing(39, 19249), evenSpacing(40, 19509)); + } + + /** + * Tests a grid where the steps differ by one pixel on a single axis. Only two of the four parts + * in which {@code Localization} would split such a grid receive points; the empty parts shall + * not cause an {@link IndexOutOfBoundsException}. + * + *

This is the case of ICEYE Single Look Complex images of 114644 × 16714 pixels: + * the tie points are spaced by 12738 pixels along x except one step of 12739 pixels, + * and evenly spaced by 1857 pixels along y.

+ * + * @throws Exception if the transform cannot be created or used. + */ + @Test + public void testGridWithIrregularStepOnOneAxis() throws Exception { + verify(roundedSpacing(114643), roundedSpacing(16713)); + } + + /** + * Tests a grid where the steps differ by one pixel on both axes. + * This is the case of ICEYE Ground Range Detected images of 20000 × 20000 pixels: + * the tie points are spaced by 2222 pixels except one step of 2223 pixels on each axis. + * + * @throws Exception if the transform cannot be created or used. + */ + @Test + public void testGridWithIrregularStepOnBothAxes() throws Exception { + verify(roundedSpacing(19999), roundedSpacing(19999)); + } + + /** + * Tests a grid where the steps alternate between two values differing by one pixel. + * This is the case of ICEYE Single Look Complex images of 34484 × 15342 pixels: + * the tie points are spaced by 3831 and 3832 pixels alternately along x, + * and by 1705 and 1704 pixels alternately along y. Splitting such a grid + * gives parts that are still irregular, so the split has to recurse. + * + * @throws Exception if the transform cannot be created or used. + */ + @Test + public void testGridWithAlternatingSteps() throws Exception { + verify(roundedSpacing(34483), roundedSpacing(15341)); + } + + /** + * Tests a grid where the last step is genuinely shorter than the other steps, + * as in Sentinel 1 images where the tie points are spaced by 1320 pixels except + * the last two which are 1302 pixels apart. Contrarily to the ICEYE grids, the + * spacing of this grid is not uniform up to a rounding to integers, so it has to + * be handled by splitting the grid in parts. + * + * @throws Exception if the transform cannot be created or used. + */ + @Test + public void testGridWithShorterLastStep() throws Exception { + final double[] coordinates = shorterLastStep(1320, 1302); + verify(coordinates, coordinates); + } + + /** + * Tests a grid where the last step is genuinely shorter on a single axis. + * This combines the Sentinel 1 spacing with the empty parts of + * {@link #testGridWithIrregularStepOnOneAxis()}. + * + * @throws Exception if the transform cannot be created or used. + */ + @Test + public void testGridWithShorterLastStepOnOneAxis() throws Exception { + final double[] rows = new double[GRID_SIZE]; + for (int i=1; i