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);
}
}
