Author: desruisseaux Date: Tue Mar 6 15:34:47 2018 New Revision: 1825995 URL: http://svn.apache.org/viewvc?rev=1825995&view=rev Log: Enable the use of Sentinel 1 localization grids in GeoTIFF reader. https://issues.apache.org/jira/browse/SIS-407
Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/geometry/AbstractEnvelope.java sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java sis/branches/JDK8/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/geometry/AbstractEnvelope.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/geometry/AbstractEnvelope.java?rev=1825995&r1=1825994&r2=1825995&view=diff ============================================================================== --- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/geometry/AbstractEnvelope.java [UTF-8] (original) +++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/geometry/AbstractEnvelope.java [UTF-8] Tue Mar 6 15:34:47 2018 @@ -695,6 +695,7 @@ public abstract class AbstractEnvelope e /** * Tests if a specified coordinate is inside the boundary of this envelope. + * Both lower and upper values of this envelope are considered inclusive. * If it least one ordinate value in the given point is {@link Double#NaN NaN}, * then this method returns {@code false}. * @@ -710,7 +711,7 @@ public abstract class AbstractEnvelope e * * @param position the point to text. * @return {@code true} if the specified coordinate is inside the boundary of this envelope; {@code false} otherwise. - * @throws MismatchedDimensionException if the specified point doesn't have the expected dimension. + * @throws MismatchedDimensionException if the specified point does not have the expected number of dimensions. * @throws AssertionError if assertions are enabled and the envelopes have mismatched CRS. */ public boolean contains(final DirectPosition position) throws MismatchedDimensionException { Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java?rev=1825995&r1=1825994&r2=1825995&view=diff ============================================================================== --- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java [UTF-8] (original) +++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java [UTF-8] Tue Mar 6 15:34:47 2018 @@ -20,11 +20,13 @@ import java.util.Map; import java.util.Arrays; import java.io.IOException; import org.opengis.util.FactoryException; +import org.opengis.geometry.Envelope; import org.opengis.geometry.DirectPosition; import org.opengis.geometry.MismatchedDimensionException; import org.opengis.geometry.coordinate.Position; import org.opengis.referencing.operation.Matrix; import org.opengis.referencing.operation.MathTransformFactory; +import org.apache.sis.geometry.GeneralEnvelope; import org.apache.sis.io.TableAppender; import org.apache.sis.math.Line; import org.apache.sis.math.Plane; @@ -339,6 +341,64 @@ search: for (int j=0; j<numPoints; j++) } /** + * Returns the envelope of source points. The lower and upper values are inclusive. + * + * @return the envelope of source points. + * @throws IllegalStateException if the source points are not yet known. + * + * @since 1.0 + */ + public Envelope getSourceEnvelope() { + if (gridSize != null) { + final int dim = gridSize.length; + final GeneralEnvelope envelope = new GeneralEnvelope(dim); + for (int i=0; i <dim; i++) { + envelope.setRange(i, 0, gridSize[i] - 1); + } + return envelope; + } else { + return envelope(sources); + } + } + + /** + * Returns the envelope of target points. The lower and upper values are inclusive. + * + * @return the envelope of target points. + * @throws IllegalStateException if the target points are not yet known. + * + * @since 1.0 + */ + public Envelope getTargetEnvelope() { + return envelope(targets); + } + + /** + * Implementation of {@link #getSourceEnvelope()} and {@link #getTargetEnvelope()}. + */ + private static Envelope envelope(final double[][] points) { + if (points == null) { + throw new IllegalStateException(noData()); + } + final int dim = points.length; + final GeneralEnvelope envelope = new GeneralEnvelope(dim); + for (int i=0; i <dim; i++) { + final double[] data = points[i]; + double lower = Double.POSITIVE_INFINITY; + double upper = Double.NEGATIVE_INFINITY; + for (final double value : data) { + if (value < lower) lower = value; + if (value > upper) upper = value; + } + if (lower > upper) { + lower = upper = Double.NaN; + } + envelope.setRange(i, lower, upper); + } + return envelope; + } + + /** * Returns the direct position of the given position, or {@code null} if none. */ private static DirectPosition position(final Position p) { @@ -569,7 +629,7 @@ search: for (int j=0; j<numPoints; j++) final double[][] sources = this.sources; // Protect from changes. final double[][] targets = this.targets; if (targets == null) { - throw new FactoryException(noData()); + throw new InvalidGeodeticParameterException(noData()); } final int sourceDim = (sources != null) ? sources.length : gridSize.length; final int targetDim = targets.length; Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java?rev=1825995&r1=1825994&r2=1825995&view=diff ============================================================================== --- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java [UTF-8] (original) +++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java [UTF-8] Tue Mar 6 15:34:47 2018 @@ -17,6 +17,7 @@ package org.apache.sis.referencing.operation.builder; import org.opengis.util.FactoryException; +import org.opengis.geometry.Envelope; import org.opengis.geometry.MismatchedDimensionException; import org.opengis.referencing.operation.Matrix; import org.opengis.referencing.operation.MathTransform; @@ -31,6 +32,7 @@ import org.apache.sis.referencing.operat import org.apache.sis.referencing.datum.DatumShiftGrid; import org.apache.sis.internal.referencing.Resources; import org.apache.sis.geometry.DirectPosition2D; +import org.apache.sis.geometry.Envelopes; import org.apache.sis.internal.util.Numerics; import org.apache.sis.util.ArgumentChecks; import org.apache.sis.measure.NumberRange; @@ -314,6 +316,25 @@ public class LocalizationGridBuilder ext } /** + * Returns the envelope of source coordinates. This is the envelope of the grid domain + * (i.e. the ranges of valid {@code gridX} and {@code gridY} argument values in calls + * to {@code get/setControlPoint(…)} methods) transformed by the inverse of + * {@linkplain #getSourceToGrid() source to grid} transform. + * The lower and upper values are inclusive. + * + * @return the envelope of grid points. + * @throws IllegalStateException if the grid points are not yet known. + * @throws TransformException if the envelope can not be calculated. + * + * @see LinearTransformBuilder#getSourceEnvelope() + * + * @since 1.0 + */ + public Envelope getSourceEnvelope() throws TransformException { + return Envelopes.transform(sourceToGrid.inverse(), linear.getSourceEnvelope()); + } + + /** * Creates a transform from the source points to the target points. * This method assumes that source points are precise and all uncertainty is in the target points. * If this transform is close enough to an affine transform, then an instance of {@link LinearTransform} is returned. Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java?rev=1825995&r1=1825994&r2=1825995&view=diff ============================================================================== --- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java [UTF-8] (original) +++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/MathTransforms.java [UTF-8] Tue Mar 6 15:34:47 2018 @@ -184,9 +184,9 @@ public final class MathTransforms extend * Creates a transform defined as one transform applied globally except in sub-areas where more accurate * transforms are available. Such constructs appear in some datum shift files. The result of transforming * a point by the returned {@code MathTransform} is as if iterating over all given {@link Envelope}s in - * no particular order, find the smallest one containing the point to transform, then use the associated - * {@link MathTransform} for transforming the point. If the point is not found in any envelope, - * then the global transform is applied. + * no particular order, find the smallest one containing the point to transform (envelope border considered + * inclusive), then use the associated {@link MathTransform} for transforming the point. + * If the point is not found in any envelope, then the global transform is applied. * * <p>The following constraints apply:</p> * <ul> Modified: sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java?rev=1825995&r1=1825994&r2=1825995&view=diff ============================================================================== --- sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java [UTF-8] (original) +++ sis/branches/JDK8/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/SpecializableTransform.java [UTF-8] Tue Mar 6 15:34:47 2018 @@ -42,6 +42,7 @@ import org.apache.sis.util.ArraysExt; /** * A transform having sub-areas where more accurate transforms can be used. * The global transform must be a reasonable approximation of the specialized transforms. + * The lower and upper values of given envelopes are inclusive. * * @author Martin Desruisseaux (Geomatys) * @version 1.0 @@ -277,7 +278,7 @@ next: for (final Map.Entry<Envelope,Ma } } for (int i=0; i<n; i++) { - if (area.intersects(areas[n])) { + if (area.intersects(areas[i])) { // Pending implementation of R-Tree in Apache SIS. throw new IllegalArgumentException("Current implementation does not accept overlapping envelopes."); } Modified: sis/branches/JDK8/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java URL: http://svn.apache.org/viewvc/sis/branches/JDK8/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java?rev=1825995&r1=1825994&r2=1825995&view=diff ============================================================================== --- sis/branches/JDK8/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java [UTF-8] (original) +++ sis/branches/JDK8/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java [UTF-8] Tue Mar 6 15:34:47 2018 @@ -18,9 +18,12 @@ package org.apache.sis.storage.geotiff; import java.util.Arrays; import java.util.Set; +import java.util.Map; import java.util.HashSet; +import java.util.LinkedHashMap; import java.util.Collection; import java.util.Collections; +import org.opengis.geometry.Envelope; import org.opengis.util.FactoryException; import org.opengis.metadata.quality.DataQuality; import org.opengis.metadata.spatial.GeolocationInformation; @@ -28,6 +31,7 @@ import org.opengis.referencing.operation import org.opengis.referencing.operation.TransformException; import org.opengis.referencing.crs.CoordinateReferenceSystem; import org.apache.sis.referencing.operation.matrix.MatrixSIS; +import org.apache.sis.referencing.operation.transform.MathTransforms; import org.apache.sis.referencing.operation.transform.LinearTransform; import org.apache.sis.referencing.operation.builder.LocalizationGridBuilder; import org.apache.sis.referencing.operation.AbstractCoordinateOperation; @@ -93,16 +97,22 @@ final class GridGeometry extends Abstrac GridGeometry(final String name, final CoordinateReferenceSystem crs, final Vector modelTiePoints) throws FactoryException, TransformException { - super(Collections.singletonMap(NAME_KEY, name), null, crs, null, localizationGrid(modelTiePoints)); + super(Collections.singletonMap(NAME_KEY, name), null, crs, null, localizationGrid(modelTiePoints, null)); } /** * Builds a localization grid from the given GeoTIFF tie points. * This method may invoke itself recursively. + * + * @param modelTiePoints the model tie points read from GeoTIFF file. + * @param addTo if non-null, add the transform result to this map. */ - private static MathTransform localizationGrid(final Vector modelTiePoints) throws FactoryException, TransformException { + private static MathTransform localizationGrid(final Vector modelTiePoints, final Map<Envelope,MathTransform> addTo) + throws FactoryException, TransformException + { final int size = modelTiePoints.size(); final int n = size / RECORD_LENGTH; + if (n == 0) return null; final Vector x = modelTiePoints.subSampling(0, RECORD_LENGTH, n); final Vector y = modelTiePoints.subSampling(1, RECORD_LENGTH, n); try { @@ -119,47 +129,94 @@ final class GridGeometry extends Abstrac modelTiePoints.doubleValue(i+4)); } grid.setDesiredPrecision(PRECISION); - return grid.create(null); + final MathTransform tr = grid.create(null); + if (addTo != null && addTo.put(grid.getSourceEnvelope(), tr) != null) { + throw new FactoryException(); // Should never happen. If it does, we have a bug in our algorithm. + } + return tr; } catch (ArithmeticException | FactoryException e) { /* * May happen when the model tie points are not distributed on a regular grid. * 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 - * are spaced by 1320 pixels and one grid for the last column. + * are spaced by 1320 pixels and one grid for the last column. Such splitting needs to be + * done horizontally and vertically, which result in four grids: + * + * ┌──────────────────┬───┐ + * │ │ │ + * │ 0 │ 1 │ + * │ │ │ + * ├──────────────────┼───┤ splitY + * │ 2 │ 3 │ + * └──────────────────┴───┘ + * splitX */ final Set<Double> uniques = new HashSet<>(100); - final double lastX = threshold(x, uniques); - final double lastY = threshold(y, uniques); - if (!Double.isNaN(lastX) || !Double.isNaN(lastY)) { - int[] indicesGrid1 = new int[size]; - int[] indicesGrid2 = new int[size]; - int sizeGrid1 = 0; - int sizeGrid2 = 0; - for (int i=0; i<size;) { - final int[] indicesGrid; - int sizeGrid; - final int s; - if (modelTiePoints.doubleValue(i) > lastX || modelTiePoints.doubleValue(i+1) > lastY) { - indicesGrid = indicesGrid2; - sizeGrid = sizeGrid2; - s = sizeGrid2 += RECORD_LENGTH; - } else { - indicesGrid = indicesGrid1; - sizeGrid = sizeGrid1; - s = sizeGrid1 += RECORD_LENGTH; + final double splitX = threshold(x, uniques); + final double splitY = threshold(y, uniques); + if (Double.isNaN(splitX) && Double.isNaN(splitY)) { + throw e; // Can not do better. Report the failure. + } + final int[][] indices = new int[4][size]; + final int[] lengths = new int[4]; + for (int i=0; i<size;) { + final double px = modelTiePoints.doubleValue(i ); + final double py = modelTiePoints.doubleValue(i+1); + int part = 0; // Number of the part where to add current point. + if (px > splitX) part = 1; // Point will be added to part #1 or #3. + if (py > splitY) part |= 2; // Point will be added to part #2 or #3. + int parts = 1 << part; // Bitmask of the parts where to add the point. + if (px == splitX) parts |= 1 << (part | 1); // Add also the point to part #1 or #3. + if (py == splitY) parts |= 1 << (part | 2); // Add also the point to part #2 or #3. + if (parts == 0b0111) { + parts = 0b1111; // Add also the point to part #3. + assert px == splitX && py == splitY; + } + final int upper = i + RECORD_LENGTH; + do { + part = Integer.numberOfTrailingZeros(parts); + @SuppressWarnings("MismatchedReadAndWriteOfArray") + final int[] tileIndices = indices[part]; + int k = lengths[part]; + for (int j=i; j<upper; j++) { + tileIndices[k++] = j; } - while (sizeGrid < s) indicesGrid[sizeGrid++] = i++; + lengths[part] = k; + } while ((parts &= ~(1 << part)) != 0); // Clear the bit of the part we processed. + i = upper; + } + /* + * At this point, we finished to collect indices of the points to use for parts #0, 1, 2 and 3. + * Verify that each part has less points than the initial vector (otherwise it would be a bug), + * and identify which part is the biggest one. This is usually part #0. + */ + int maxLength = 0; + int largestPart = 0; + for (int i=0; i<indices.length; i++) { + final int length = lengths[i]; + if (length >= size) throw e; // Safety against infinite recursivity. + indices[i] = Arrays.copyOf(indices[i], length); + if (length > maxLength) { + maxLength = length; + largestPart = i; } - if (sizeGrid1 < size && sizeGrid2 < size) { // Paranoiac check against never-ending recursivity. - indicesGrid1 = Arrays.copyOf(indicesGrid1, sizeGrid1); - indicesGrid2 = Arrays.copyOf(indicesGrid2, sizeGrid2); - final MathTransform grid1 = localizationGrid(modelTiePoints.pick(indicesGrid1)); - final MathTransform grid2 = localizationGrid(modelTiePoints.pick(indicesGrid2)); - // TODO: pending completing of SIS-408. + } + /* + * The biggest part will define the global transform. All other parts will define a specialization + * 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); + 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); } } - throw e; + return MathTransforms.specialize(global, specialization); } }