From c3d4ab228ee12f806f77e6919f7f42beda5f27af Mon Sep 17 00:00:00 2001 From: JF3Env Date: Mon, 24 Aug 2026 17:48:40 +0200 Subject: [PATCH] fix(geotiff): build localization grid for rectilinear tie points MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit `Localization.nonLinear(…)` cannot build a localization grid when the model tie points are on a rectilinear grid, i.e. when every combination of the distinct x and y pixel coordinates is present exactly once, but the spacing between those coordinates is not an integer amount of pixels repeated exactly. `LocalizationGridBuilder` infers the grid size from the greatest common divisor of the tie point coordinates. That divisor is much smaller than the actual step when the step is fractional, and it collapses to 1 when a rounding to integers makes one step differ from the others by one pixel. Either way the inferred grid becomes as large as the image and the constructor throws `ArithmeticException`. The recovery heuristic then splits the tie points in four parts around a threshold computed per axis. That heuristic was written for Sentinel 1 images, whose last step is genuinely 18 pixels shorter than the others, and it behaves badly on a grid that is rectilinear: * If the spacing is irregular on a single axis, `threshold(…)` returns NaN for the other axis, the comparisons against it are always false, two of the four parts receive no point, and `Vector.pick(int[])` throws `IndexOutOfBoundsException: Range [0 … -1] is not valid`. * If a part is itself irregular it is split again, but the transform computed by the recursive call is added to the map of specializations only when it could be built without splitting. Otherwise the sub-area silently falls back to an extrapolation of the largest part, far from the tie points it was supposed to honor. This affects every ICEYE product. They carry no `ModelPixelScaleTag` and no `ModelTransformationTag`, so they always take this code path, and their tie points are at k × (size-1) / (n-1) pixels, sometimes rounded to integers. Such grids are now built without splitting: the ranks of the distinct pixel coordinates are used as grid indices, and the linear relationship between pixel coordinates and ranks is applied before the localization grid. The resulting transform is continuous and honors every tie point. Grids whose spacing is irregular for a real reason, such as Sentinel 1, are detected by `isUniformAfterRounding(…)` and are still handled by splitting, unchanged. The empty parts of the split are skipped, so the exception above cannot happen anymore even for the grids that are still split. Measured on five ICEYE products, the worst distance between a tie point as declared in the file and as computed by the transform. "Before" is the same code with only the empty parts skipped, since without that the products cannot be opened at all. About 0.2 m of the remaining error on the GRD and SLC products, and about 4 m on the ScanSAR product, is the half pixel between the tie point convention of those files and `PixelInCell.CELL_CORNER`; it is unchanged by this commit. product size spacing before after ScanSAR 19250 × 19510 fractional 2161.51 m 4.55 m GRD 20000 × 20000 rounded, both 2.09 m 0.33 m GRD 20000 × 20000 rounded, both 2.64 m 0.33 m SLC 34484 × 15342 rounded, alt. 17.09 m 0.34 m SLC 114644 × 16714 rounded, one 2.23 m 0.16 m The ScanSAR product is the sample published by ESA at https://earth.esa.int/eogateway/ftp/missions/sample-data/third-party-missions/iceye/ICEYE-Scan-mode.zip so this can be reproduced without an ICEYE licence. `LocalizationTest` covers the spacings of those products, the Sentinel 1 spacing, and a regular grid as a control. The tie points are built from a non-linear model, otherwise every candidate transform would reproduce them exactly and the tests would only verify the absence of exception. Co-Authored-By: Claude Opus 5 Claude-Session: https://claude.ai/code/session_01AG48J4hdPvWYWd3hFrhZGz --- .../storage/geotiff/reader/Localization.java | 126 +++++++++ .../geotiff/reader/LocalizationTest.java | 264 ++++++++++++++++++ 2 files changed, 390 insertions(+) create mode 100644 endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/reader/LocalizationTest.java 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