This is an automated email from the ASF dual-hosted git repository.
asf-gitbox-commits pushed a commit to branch geoapi-4.0
in repository https://gitbox.apache.org/repos/asf/sis.git
The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
new c0abcda167 Make localization grid tolerant to coordinates rounded to
nearest integers. https://github.com/apache/sis/pull/48
c0abcda167 is described below
commit c0abcda16726b5215e443b1cdfeb73290d6214cf
Author: Martin Desruisseaux <[email protected]>
AuthorDate: Tue Oct 6 20:40:55 2026 +0200
Make localization grid tolerant to coordinates rounded to nearest integers.
https://github.com/apache/sis/pull/48
Co-authored-by: JF3Env <[email protected]>
---
.../operation/builder/LocalizationGridBuilder.java | 140 +++++++----
.../operation/builder/package-info.java | 2 +-
.../builder/LocalizationGridBuilderTest.java | 42 +++-
.../sis/storage/geotiff/reader/Localization.java | 23 +-
.../storage/geotiff/reader/LocalizationTest.java | 270 +++++++++++++++++++++
.../main/org/apache/sis/math/Vector.java | 15 +-
6 files changed, 429 insertions(+), 63 deletions(-)
diff --git
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
index e635f75e4f..b775adee7e 100644
---
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
+++
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
@@ -87,7 +87,7 @@ import org.opengis.coordinate.MismatchedDimensionException;
* See the <cite>Linearizers</cite> section in {@link LinearTransformBuilder}
for more discussion.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.2
+ * @version 1.7
*
* @see InterpolatedTransform
* @see LinearTransform
@@ -96,11 +96,6 @@ import org.opengis.coordinate.MismatchedDimensionException;
* @since 0.8
*/
public class LocalizationGridBuilder extends TransformBuilder {
- /**
- * Tolerance threshold for comparing pixel coordinates relative to integer
values.
- */
- private static final double EPS = Numerics.COMPARISON_THRESHOLD;
-
/**
* The transform for the linear part.
* Always created with a grid size specified to the constructor.
@@ -116,6 +111,8 @@ public class LocalizationGridBuilder extends
TransformBuilder {
/**
* Conversions from source real-world coordinates to grid indices before
interpolation.
* If there is no such conversion to apply, then this is the identity
transform.
+ *
+ * @see #getSourceToGrid()
*/
private LinearTransform sourceToGrid;
@@ -180,9 +177,9 @@ public class LocalizationGridBuilder extends
TransformBuilder {
* @throws ArithmeticException if this constructor cannot infer a
reasonable grid size from the given vectors.
*/
public LocalizationGridBuilder(final Vector sourceX, final Vector sourceY)
{
- final Matrix fromGrid = new Matrix3();
- final int width = infer(sourceX, fromGrid, 0);
- final int height = infer(sourceY, fromGrid, 1);
+ final var fromGrid = new Matrix3();
+ final int width = infer(sourceX, fromGrid, 0);
+ final int height = infer(sourceY, fromGrid, 1);
linearBuilder = new LinearTransformBuilder(width, height);
try {
sourceToGrid = MathTransforms.linear(fromGrid).inverse();
@@ -227,7 +224,7 @@ public class LocalizationGridBuilder extends
TransformBuilder {
final Vector[] sources = localizations.sources();
n = sources.length;
if (n == SOURCE_DIMENSION) {
- final Matrix fromGrid = new Matrix3();
+ final var fromGrid = new Matrix3();
final int width = infer(sources[0], fromGrid, 0);
final int height = infer(sources[1], fromGrid, 1);
linearBuilder = new LinearTransformBuilder(width, height);
@@ -247,35 +244,83 @@ public class LocalizationGridBuilder extends
TransformBuilder {
}
/**
- * Infers a grid size by searching for the greatest common divisor (GCD)
for values in the given vector.
+ * Infers a grid size by searching for a constant increment between all
values in the given vector.
* The vector values should be integers, but this method is tolerant to
constant offsets (typically 0.5).
- * The GCD is taken as a "grid to source" scale factor and the minimal
value as the translation term.
+ * The increment is taken as a "grid to source" scale factor and the
minimal value as the translation term.
* Those two values are stored in the {@code dim} row of the given matrix.
*
- * @param source the vector of values for which to get the GCD and
minimum value.
- * @param fromGrid matrix where to store the minimum value and the GCD.
+ * @param source the vector of values for which to get the increment
and minimum value.
+ * @param fromGrid matrix where to store the minimum value and the
increment.
* @param dim index of the matrix row to update.
* @return grid size.
+ * @throws ArithmeticException if the grid size cannot be computed.
*/
private static int infer(final Vector source, final Matrix fromGrid, final
int dim) {
final NumberRange<?> range = source.range();
- final double min = range.getMinDouble(true);
- final double span = range.getSpan();
- final Number increment = source.increment(EPS * span);
- double inc;
- if (increment != null) {
- inc = increment.doubleValue();
+ final double min = range.getMinDouble(true);
+ final double span = range.getSpan();
+ final int size = source.size();
+ final Number scale = source.increment(span / (size - 1) *
DEFAULT_PRECISION);
+ double increment;
+ if (scale != null) {
+ increment = Math.abs(scale.doubleValue());
} else {
- inc = span;
- final int size = source.size();
+ /*
+ * Initialize the increment to the difference between the two
values closest to zero.
+ * They are likely to be the two most accurate values, thus
reducing rounding errors.
+ * The block is for keeping variables in local scope.
+ */
+ {
+ // Find the value closest to zero.
+ double zero = Double.POSITIVE_INFINITY;
+ double abs0 = Double.POSITIVE_INFINITY;
+ for (int i=0; i<size; i++) {
+ final double value = source.doubleValue(i);
+ final double abs = Math.abs(value);
+ if (abs < abs0) {
+ abs0 = abs;
+ zero = value;
+ if (zero == 0) break; // Optimization for a
common case.
+ }
+ }
+ // Find the second value closest to zero.
+ double tolerance = abs0 + span / (Math.sqrt(size) - 1) *
DEFAULT_PRECISION; // Assuming a square grid.
+ double closeZero = Double.POSITIVE_INFINITY;
+ abs0 = Double.POSITIVE_INFINITY;
+ for (int i=0; i<size; i++) {
+ final double value = source.doubleValue(i);
+ final double abs = Math.abs(value);
+ if (abs < abs0 && abs > tolerance) {
+ abs0 = abs;
+ closeZero = value;
+ }
+ }
+ increment = Math.abs(closeZero - zero);
+ /*
+ * Use the increment for computing the expected number of
distinct values.
+ * Then, re-compute the increment as the average delta over
all the span.
+ * It will be the same increment if the delta between values
is constant,
+ * but may differ if the delta is not constant.
+ */
+ final double average = span / Math.rint(span / increment);
+ if (Math.abs(average - increment) > span * DEFAULT_PRECISION
|| Numerics.isInteger(average)) {
+ increment = average;
+ }
+ }
+ /*
+ * If the vector contains only integer values, it may be because
all coordinates have been rounded
+ * toward nearest integer. In such case, tolerate a difference
corresponding to the rounding errors.
+ * Then verifies that all values in the vector are multiples of
the increment (ignoring `min` offset).
+ */
+ double tolerance = increment * DEFAULT_PRECISION;
+ if (tolerance < 0.5 && source.isInteger()) {
+ tolerance = 0.5;
+ }
for (int i=0; i<size; i++) {
- double v = source.doubleValue(i) - min;
- if (Math.abs(v % inc) > EPS) {
- do {
- final double r = (inc % v); // Both `inc` and `v`
are positive, so `r` will be positive too.
- inc = v;
- v = r;
- } while (Math.abs(v) > EPS);
+ final double remainder = (source.doubleValue(i) - min) %
increment;
+ if (remainder > tolerance && (increment - remainder) >
tolerance) {
+ increment = 0; // Will cause an exception to be thrown
below.
+ break;
}
}
}
@@ -285,15 +330,22 @@ public class LocalizationGridBuilder extends
TransformBuilder {
* would fail anyway, because it would contain holes where no value is
defined. A limit is important
* for preventing useless allocation of large arrays -
https://issues.apache.org/jira/browse/SIS-407
*/
- fromGrid.setElement(dim, dim, inc);
+ fromGrid.setElement(dim, dim, increment);
fromGrid.setElement(dim, SOURCE_DIMENSION, min);
- final double n = span / inc;
+ final double n = span / increment;
if (n >= 0.5 && n < source.size() - 0.5) { // Compare as
`double` in case the value is large.
return ((int) Math.round(n)) + 1;
}
throw new
ArithmeticException(Resources.format(Resources.Keys.CanNotInferGridSizeFromValues_1,
range));
}
+ /**
+ * Returns the grid size for the given dimension.
+ */
+ final int gridSize(final int srcDim) {
+ return linearBuilder.gridSize(srcDim);
+ }
+
/**
* Throws {@link IllegalStateException} if this builder cannot be modified
anymore.
*/
@@ -657,11 +709,11 @@ public class LocalizationGridBuilder extends
TransformBuilder {
if (isExact) {
step = MathTransforms.concatenate(sourceToGrid, gridToCoord);
} else {
- final int width = linearBuilder.gridSize(0);
- final int height = linearBuilder.gridSize(1);
- final float[] residual = new float [SOURCE_DIMENSION *
linearBuilder.gridLength];
- final double[] grid = new double[SOURCE_DIMENSION * width];
- double gridPrecision = precision;
+ final int width = gridSize(0);
+ final int height = gridSize(1);
+ final var residual = new float [SOURCE_DIMENSION *
linearBuilder.gridLength];
+ final var grid = new double[SOURCE_DIMENSION * width];
+ double gridPrecision = precision;
try {
/*
* If the user specified a precision, we need to convert
it from source units to grid units.
@@ -761,7 +813,7 @@ public class LocalizationGridBuilder extends
TransformBuilder {
*
* @since 1.1
*/
- public Optional<Map.Entry<String,MathTransform>> linearizer(final boolean
ifNotCompensated) {
+ public Optional<Map.Entry<String, MathTransform>> linearizer(final boolean
ifNotCompensated) {
ProjectedTransformTry linearizer = linearBuilder.appliedLinearizer();
if (ifNotCompensated && linearizer != null &&
linearizer.reverseAfterLinearization) {
linearizer = null;
@@ -803,10 +855,10 @@ public class LocalizationGridBuilder extends
TransformBuilder {
if (mt == null) {
throw new
IllegalStateException(Errors.format(Errors.Keys.Uninitialized_1,
getClass().getSimpleName()));
}
- final int tgtDim = mt.getTargetDimensions();
- final double[] point = new double[Math.max(tgtDim,
SOURCE_DIMENSION)];
- final Statistics[] stats = new Statistics[tgtDim + SOURCE_DIMENSION];
- final StringBuilder buffer = new StringBuilder();
+ final int tgtDim = mt.getTargetDimensions();
+ final var point = new double[Math.max(tgtDim, SOURCE_DIMENSION)];
+ final var stats = new Statistics[tgtDim + SOURCE_DIMENSION];
+ final var buffer = new StringBuilder();
for (int i=0; i<stats.length; i++) {
buffer.setLength(0);
buffer.append('Δ');
@@ -835,8 +887,8 @@ public class LocalizationGridBuilder extends
TransformBuilder {
} catch (NoninvertibleTransformException e) {
throw new IllegalStateException(e);
}
- final int width = linearBuilder.gridSize(0);
- final int height = linearBuilder.gridSize(1);
+ final int width = gridSize(0);
+ final int height = gridSize(1);
for (int y=0; y<height; y++) {
for (int x=0; x<width; x++) {
point[0] = gridCoordinates[0] = x;
@@ -893,7 +945,7 @@ public class LocalizationGridBuilder extends
TransformBuilder {
* @since 1.1
*/
public String toString(final boolean linear, final Locale locale) {
- final StringBuilder buffer = new StringBuilder(400);
+ final var buffer = new StringBuilder(400);
String lineSeparator = null;
try {
lineSeparator = linearBuilder.appendTo(buffer, getClass(), locale,
Vocabulary.Keys.LinearTransformation);
diff --git
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/package-info.java
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/package-info.java
index a21d184cd3..39994212b8 100644
---
a/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/package-info.java
+++
b/endorsed/src/org.apache.sis.referencing/main/org/apache/sis/referencing/operation/builder/package-info.java
@@ -30,7 +30,7 @@
* convenience.</p>
*
* @author Martin Desruisseaux (Geomatys)
- * @version 1.2
+ * @version 1.7
* @since 0.5
*/
package org.apache.sis.referencing.operation.builder;
diff --git
a/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/builder/LocalizationGridBuilderTest.java
b/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/builder/LocalizationGridBuilderTest.java
index 9ffdeddf1a..3b05075caa 100644
---
a/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/builder/LocalizationGridBuilderTest.java
+++
b/endorsed/src/org.apache.sis.referencing/test/org/apache/sis/referencing/operation/builder/LocalizationGridBuilderTest.java
@@ -16,11 +16,14 @@
*/
package org.apache.sis.referencing.operation.builder;
+import java.util.Arrays;
import java.awt.geom.Point2D;
import java.awt.geom.AffineTransform;
import org.opengis.util.FactoryException;
+import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.geometry.Envelope2D;
+import org.apache.sis.math.Vector;
// Test dependencies
import org.junit.jupiter.api.Test;
@@ -39,6 +42,12 @@ import static
org.apache.sis.referencing.Assertions.assertEnvelopeEquals;
@SuppressWarnings("exports")
@ExtendWith(FailureDetailsReporter.class)
public final class LocalizationGridBuilderTest extends TransformTestCase {
+ /**
+ * Whether to generate assertion codes.
+ * If enabled, the code is sent to the standard output stream.
+ */
+ private static final boolean GENERATE_TEST_CODE = false;
+
/**
* Creates a new test case.
*/
@@ -56,7 +65,7 @@ public final class LocalizationGridBuilderTest extends
TransformTestCase {
*/
@SuppressWarnings("UseOfSystemOutOrSystemErr")
private static LocalizationGridBuilder builder(final AffineTransform
reference, final int width, final int height) {
- final LocalizationGridBuilder builder = new
LocalizationGridBuilder(width, height);
+ final var builder = new LocalizationGridBuilder(width, height);
Point2D pt = new Point2D.Double();
for (int gridY=0; gridY < height; gridY++) {
for (int gridX=0; gridX < width; gridX++) {
@@ -67,7 +76,7 @@ public final class LocalizationGridBuilderTest extends
TransformTestCase {
final double x = pt.getX() + 0.4*gx2 + 0.7*gy2;
final double y = pt.getY() + 0.3*gx2 - 0.5*gy2;
builder.setControlPoint(gridX, gridY, x, y);
- if (false) {
+ if (GENERATE_TEST_CODE) {
// For generating verification code.
System.out.printf("verifyTransform(new double[] {%d, %d},
new double[] {%f, %f});%n", gridX, gridY, x, y);
}
@@ -84,7 +93,7 @@ public final class LocalizationGridBuilderTest extends
TransformTestCase {
*/
@Test
public void testQuadratic() throws FactoryException, TransformException {
- final AffineTransform reference = new AffineTransform(20, -30, 5, -4,
-20, 8);
+ final var reference = new AffineTransform(20, -30, 5, -4, -20, 8);
final LocalizationGridBuilder builder = builder(reference, 5, 4);
builder.setDesiredPrecision(1E-6);
transform = builder.create(null);
@@ -122,14 +131,14 @@ public final class LocalizationGridBuilderTest extends
TransformTestCase {
*/
@Test
public void testCreateFromLocalizations() throws TransformException {
- final LinearTransformBuilder localizations = new
LinearTransformBuilder();
+ final var localizations = new LinearTransformBuilder();
localizations.setControlPoint(new int[] {0, 0}, new double[] {-20.0,
8.0});
localizations.setControlPoint(new int[] {1, 0}, new double[] { 0.4,
-21.7});
localizations.setControlPoint(new int[] {0, 1}, new double[] {-14.3,
3.5});
localizations.setControlPoint(new int[] {1, 1}, new double[] { 6.1,
-26.2});
localizations.setControlPoint(new int[] {0, 2}, new double[] { 1.3,
-8.5});
localizations.setControlPoint(new int[] {1, 2}, new double[] { 87.7,
-123.7});
- LocalizationGridBuilder builder = new
LocalizationGridBuilder(localizations);
+ final var builder = new LocalizationGridBuilder(localizations);
/*
* Verifies the grid size by checking the source envelope.
* Minimum and maximum values are inclusive.
@@ -148,4 +157,27 @@ public final class LocalizationGridBuilderTest extends
TransformTestCase {
assertArrayEquals(new double[] {-8.5, -123.7}, builder.getRow(1,
2).doubleValues());
assertArrayEquals(new double[] {-21.7, -26.2, -123.7},
builder.getColumn(1, 1).doubleValues());
}
+
+ /**
+ * Tests inferring the grid size from the vectors of <var>x</var> and
<var>y</var> values.
+ * This test uses a non-integer delta between grid coordinates in order to
test robustness
+ * against rounding errors.
+ */
+ @Test
+ public void testInferGridSize() {
+ final var x = new double[39];
+ final var y = new double[40];
+ final double sx = 19249d/38d;
+ final double sy = 19509/39d;
+ Arrays.setAll(x, (i) -> -1d/3d + (i % 10) * sx + 1E-12 *
StrictMath.random());
+ Arrays.setAll(y, (i) -> -1d/6d + (i / 10) * sy + 1E-12 *
StrictMath.random());
+ final var builder = new LocalizationGridBuilder(Vector.create(x),
Vector.create(y));
+ assertEquals(10, builder.gridSize(0), "width");
+ assertEquals( 4, builder.gridSize(1), "height");
+ final Matrix m = builder.getSourceToGrid().getMatrix();
+ assertEquals(3, m.getNumCol());
+ assertEquals(3, m.getNumRow());
+ assertEquals(1/sx, m.getElement(0, 0), 1E-16);
+ assertEquals(1/sy, m.getElement(1, 1), 1E-16);
+ }
}
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..731c24b302 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
@@ -40,6 +40,7 @@ import org.apache.sis.math.Vector;
* if another data store needs similar functionality in the future.
*
* @author Martin Desruisseaux (Geomatys)
+ * @author Jonatas Fischer
*/
final class Localization {
/**
@@ -80,7 +81,7 @@ final class Localization {
* @param addTo if non-null, add the transform result to this
map.
* @return the "grid to CRS" transform backed by the localization grid.
*/
- private static MathTransform localizationGrid(final Vector modelTiePoints,
final Map<Envelope,MathTransform> addTo)
+ private static MathTransform localizationGrid(final Vector modelTiePoints,
final Map<Envelope, MathTransform> addTo)
throws FactoryException, TransformException
{
final int size = modelTiePoints.size();
@@ -89,7 +90,7 @@ final class Localization {
final Vector x = modelTiePoints.subSampling(0, RECORD_LENGTH, n);
final Vector y = modelTiePoints.subSampling(1, RECORD_LENGTH, n);
try {
- final LocalizationGridBuilder grid = new
LocalizationGridBuilder(x, y);
+ final var grid = new LocalizationGridBuilder(x, y);
final LinearTransform sourceToGrid = grid.getSourceToGrid();
final double[] coordinates = new double[2];
for (int i=0; i<size; i += RECORD_LENGTH) {
@@ -124,8 +125,12 @@ final class Localization {
* │ 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<Double> uniques = new HashSet<>(100);
+ final var uniques = new HashSet<Double>(100);
final double splitX = threshold(x, uniques);
final double splitY = threshold(y, uniques);
if (Double.isNaN(splitX) && Double.isNaN(splitY)) {
@@ -180,13 +185,15 @@ final class Localization {
* valid only in a sub-area. Put those information in a map for
MathTransforms.specialize(…).
*/
MathTransform global = null;
- final Map<Envelope,MathTransform> specialization = new
LinkedHashMap<>(4);
+ final var specialization = new LinkedHashMap<Envelope,
MathTransform>(4);
for (int i=0; i<indices.length; i++) {
final Vector sub = modelTiePoints.pick(indices[i]);
- if (i == largestPart) {
- global = localizationGrid(sub, null);
- } else {
- localizationGrid(sub, specialization);
+ if (!sub.isEmpty()) {
+ if (i == largestPart) {
+ global = localizationGrid(sub, null);
+ } else {
+ localizationGrid(sub, specialization);
+ }
}
}
return MathTransforms.specialize(global, specialization);
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..111ea6a9a9
--- /dev/null
+++
b/endorsed/src/org.apache.sis.storage.geotiff/test/org/apache/sis/storage/geotiff/reader/LocalizationTest.java
@@ -0,0 +1,270 @@
+/*
+ * 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.
+ *
+ * <p>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.</p>
+ *
+ * @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 <abbr>ICEYE</abbr> 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 coordinate is zero and the last coordinate is the given value.
+ * The coordinates 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];
+ final double step = last / (count - 1.0);
+ for (int i=0; i<count; i++) {
+ coordinates[i] = i * step;
+ }
+ return coordinates;
+ }
+
+ /**
+ * Returns the pixel coordinates of {@value #GRID_SIZE} evenly spaced tie
points,
+ * rounded to integers as in the GeoTIFF files written by
<abbr>ICEYE</abbr>.
+ * Consequently the steps may differ from each other by one pixel.
+ *
+ * @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[] roundedSpacing(final int last) {
+ final double[] coordinates = evenSpacing(GRID_SIZE, last);
+ for (int i=0; i<coordinates.length; i++) {
+ coordinates[i] = Math.round(coordinates[i]);
+ }
+ return coordinates;
+ }
+
+ /**
+ * Returns the pixel coordinates of {@value #GRID_SIZE} tie points spaced
by the given step,
+ * except the last point which is closer to its predecessor. This is the
spacing of Sentinel 1
+ * images, and the reason why {@link Localization} splits irregular grids
in four parts.
+ *
+ * @param step step between two consecutive tie points, except the last
two.
+ * @param last step between the two last tie points.
+ * @return pixel coordinates of the tie points.
+ */
+ private static double[] shorterLastStep(final int step, final int last) {
+ final double[] coordinates = new double[GRID_SIZE];
+ for (int i=0; i<GRID_SIZE; i++) {
+ coordinates[i] = i * step;
+ }
+ coordinates[GRID_SIZE-1] =last;
+ return coordinates;
+ }
+
+ /**
+ * Creates tie points on the grid formed by the given pixel coordinates.
Model coordinates are computed
+ * by an arbitrary non-linear function of the pixel coordinates, in order
to give the localization grid
+ * something to interpolate. The magnitude of the non-linear terms is
about 10 metres,
+ * which is the order of magnitude of the terrain-induced distortion in a
radar image.
+ *
+ * @param columns pixel coordinates of the tie points along the
<var>x</var> axis.
+ * @param rows pixel coordinates of the tie points along the
<var>y</var> 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
<var>x</var> axis.
+ * @param rows pixel coordinates of the tie points along the
<var>y</var> 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; i<tiePoints.size(); i += Localization.RECORD_LENGTH) {
+ source[0] = tiePoints.doubleValue(i);
+ source[1] = tiePoints.doubleValue(i+1);
+ gridToCRS.transform(source, 0, target, 0, 1);
+ assertEquals(tiePoints.doubleValue(i+3), target[0], TOLERANCE,
"Longitude of tie point");
+ assertEquals(tiePoints.doubleValue(i+4), target[1], TOLERANCE,
"Latitude of tie point");
+ }
+ }
+
+ /**
+ * Tests a grid where all points are evenly spaced by an integer amount of
pixels. This is the case
+ * handled directly by {@link
org.apache.sis.referencing.operation.builder.LocalizationGridBuilder},
+ * without any of the fallbacks tested by the other methods.
+ *
+ * @throws Exception if the transform cannot be created or used.
+ */
+ @Test
+ public void testRegularGrid() throws Exception {
+ final double[] coordinates = new double[GRID_SIZE];
+ for (int i=1; i<GRID_SIZE; i++) {
+ coordinates[i] = coordinates[i-1] + 2222;
+ }
+ verify(coordinates, coordinates);
+ }
+
+ /**
+ * Tests a grid where all points are evenly spaced, but by a fractional
amount of pixels.
+ * The greatest common divisor of those coordinates is much smaller than
the actual step,
+ * so the grid size cannot be inferred from it.
+ *
+ * <p>This is the case of <abbr>ICEYE</abbr> ScanSAR images of 19250 ×
19510 pixels,
+ * which have 39 × 40 tie points spaced by 506.55 and 500.23 pixels
respectively.</p>
+ *
+ * @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}.
+ *
+ * <p>This is the case of ICEYE Single Look Complex images of 114644 ×
16714 pixels:
+ * the tie points are spaced by 12738 pixels along <var>x</var> except one
step of 12739 pixels,
+ * and evenly spaced by 1857 pixels along <var>y</var>.</p>
+ *
+ * @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
<var>x</var>,
+ * and by 1705 and 1704 pixels alternately along <var>y</var>. 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<GRID_SIZE; i++) {
+ rows[i] = rows[i-1] + 1320;
+ }
+ verify(shorterLastStep(1320, 1302), rows);
+ }
+}
diff --git
a/endorsed/src/org.apache.sis.util/main/org/apache/sis/math/Vector.java
b/endorsed/src/org.apache.sis.util/main/org/apache/sis/math/Vector.java
index 71af8410fb..517d867515 100644
--- a/endorsed/src/org.apache.sis.util/main/org/apache/sis/math/Vector.java
+++ b/endorsed/src/org.apache.sis.util/main/org/apache/sis/math/Vector.java
@@ -90,7 +90,7 @@ import org.apache.sis.system.Loggers;
* method and by accepting buffer in the {@link #create(Object, boolean)}
method.
*
* @author Martin Desruisseaux (MPO, Geomatys)
- * @version 1.6
+ * @version 1.7
*
* @see org.apache.sis.util.collection.IntegerList
*
@@ -912,12 +912,17 @@ search: for (;;) {
*/
@SuppressWarnings("ReturnOfCollectionOrArrayField")
public Vector subSampling(final int first, final int step, final int
length) {
- final int size = size();
- if (step == 1 && first == 0 && length == size) {
+ int limit = size();
+ if (step == 1 && first == 0 && length == limit) {
return this;
}
- final long last = first + step * (length - 1L);
- if (first < 0 || first >= size || last < 0 || last >= size || length <
0) {
+ ArgumentChecks.ensurePositive("length", length);
+ long last = first;
+ if (length != 0) {
+ last += step * (length - 1L);
+ limit--;
+ }
+ if ((first | last) < 0 || Math.max(first, last) > limit) {
final short key;
final Object arg1, arg2;
if (step == 1) {
