Skip to content

Fix GeoTIFF localization grid when tie points are rectilinear but unevenly spaced #48

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
JF3Env wants to merge 1 commit into
base: main
Choose a base branch
from
Open

Fix GeoTIFF localization grid when tie points are rectilinear but unevenly spaced #48

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
126 changes: 126 additions & 0 deletions ...g.apache.sis.storage.geotiff/main/org/apache/sis/storage/geotiff/reader/Localization.java
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down Expand Up @@ -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
Expand All @@ -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<Double> uniques = new HashSet<>(100);
final double splitX = threshold(x, uniques);
Expand Down Expand Up @@ -182,6 +204,7 @@ private static MathTransform localizationGrid(final Vector modelTiePoints, final
MathTransform global = null;
final Map<Envelope,MathTransform> specialization = new LinkedHashMap<>(4);
for (int i=0; i<indices.length; i++) {
if (indices[i].length == 0) continue; // Part without points (see above comment).
final Vector sub = modelTiePoints.pick(indices[i]);
if (i == largestPart) {
global = localizationGrid(sub, null);
Expand All @@ -193,6 +216,109 @@ private static MathTransform localizationGrid(final Vector modelTiePoints, final
}
}

/**
* Builds a single localization grid when the given tie points are on a rectilinear grid.
* The tie points are on a rectilinear grid if every combination of the distinct <var>x</var>
* and <var>y</var> 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)}.
*
* <p>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.</p>
*
* @param modelTiePoints the model tie points read from GeoTIFF file.
* @param x the <var>x</var> pixel coordinates of the tie points.
* @param y the <var>y</var> 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<Envelope,MathTransform> 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; i<size; i += RECORD_LENGTH) {
final int gx = Arrays.binarySearch(ux, modelTiePoints.doubleValue(i ));
final int gy = Arrays.binarySearch(uy, modelTiePoints.doubleValue(i+1));
if (gx < 0 || gy < 0) {
return null; // Paranoiac safety; should never happen.
}
grid.setControlPoint(gx, gy, modelTiePoints.doubleValue(i + (RECORD_LENGTH/2)),
modelTiePoints.doubleValue(i + (RECORD_LENGTH/2 + 1)));
}
grid.setDesiredPrecision(PRECISION);
final MathTransform tr = MathTransforms.concatenate(
MathTransforms.linear(new Matrix3(1/sx, 0, -ux[0]/sx,
0, 1/sy, -uy[0]/sy,
0, 0, 1)),
grid.create(null));
if (addTo != null) {
final Envelope domain = new GeneralEnvelope(new double[] {ux[0], uy[0]},
new double[] {ux[nx-1], uy[ny-1]});
if (addTo.put(domain, tr) != null) {
throw new InternalFactoryException(); // Should never happen. If it does, we have a bug in our algorithm.
}
}
return tr;
}

/**
* Returns the sorted distinct values of the given vector.
*
* @param values the <var>x</var> or <var>y</var> 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<Double> uniques = new HashSet<>(100);
for (int i=0; i<n; i++) {
uniques.add(values.doubleValue(i));
}
final double[] array = new double[uniques.size()];
int i = 0;
for (final Double value : uniques) {
array[i++] = value;
}
Arrays.sort(array);
return array;
}

/**
* Returns whether the given sorted distinct coordinates are evenly spaced by the given step,
* up to a rounding to integers. This is the case of ICEYE images, where the tie points are at
* {@code k × (size-1) / (n-1)} pixels, sometimes rounded to integers: the coordinates deviate
* from an evenly spaced sequence by less than one pixel, and the resulting error on the
* interpolated coordinates is a small fraction of a pixel. This is not the case of Sentinel 1
* images, where the last step is genuinely 18 pixels shorter than the others and where
* splitting the grid gives a better result.
*
* @param values the distinct pixel coordinates, in increasing order.
* @param step the average step between two consecutive values.
* @return whether the values are evenly spaced up to a rounding to integers.
*/
private static boolean isUniformAfterRounding(final double[] values, final double step) {
for (int i = values.length - 1; --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
Expand Down
Loading