You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2018/03/08 15:44:56 UTC
svn commit: r1826248 [3/3] - in /sis/trunk: ./
application/sis-console/src/main/java/org/apache/sis/console/
core/sis-metadata/src/main/java/org/apache/sis/internal/metadata/
core/sis-metadata/src/main/java/org/apache/sis/io/wkt/
core/sis-referencing/s...
Modified: sis/trunk/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java
URL: http://svn.apache.org/viewvc/sis/trunk/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java?rev=1826248&r1=1826247&r2=1826248&view=diff
==============================================================================
--- sis/trunk/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java [UTF-8] (original)
+++ sis/trunk/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometry.java [UTF-8] Thu Mar 8 15:44:56 2018
@@ -16,8 +16,14 @@
*/
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;
@@ -25,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;
@@ -67,7 +74,7 @@ import org.apache.sis.math.Vector;
* are considered equivalent.
*
* @author Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.0
* @since 0.8
* @module
*/
@@ -90,34 +97,161 @@ 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 is a workaround for RFE #4093999 in Sun's bug database
- * ("Relax constraint on placement of this()/super() call in constructors").
+ * 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;
- final LocalizationGridBuilder grid = new LocalizationGridBuilder(
- modelTiePoints.subSampling(0, RECORD_LENGTH, n),
- modelTiePoints.subSampling(1, RECORD_LENGTH, n));
-
- final LinearTransform sourceToGrid = grid.getSourceToGrid();
- final double[] ordinates = new double[2];
- for (int i=0; i<size; i += RECORD_LENGTH) {
- ordinates[0] = modelTiePoints.doubleValue(i);
- ordinates[1] = modelTiePoints.doubleValue(i+1);
- sourceToGrid.transform(ordinates, 0, ordinates, 0, 1);
- grid.setControlPoint((int) Math.round(ordinates[0]),
- (int) Math.round(ordinates[1]),
- modelTiePoints.doubleValue(i+3),
- modelTiePoints.doubleValue(i+4));
+ if (n == 0) return null;
+ 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 LinearTransform sourceToGrid = grid.getSourceToGrid();
+ final double[] ordinates = new double[2];
+ for (int i=0; i<size; i += RECORD_LENGTH) {
+ ordinates[0] = modelTiePoints.doubleValue(i);
+ ordinates[1] = modelTiePoints.doubleValue(i+1);
+ sourceToGrid.transform(ordinates, 0, ordinates, 0, 1);
+ grid.setControlPoint(Math.toIntExact(Math.round(ordinates[0])),
+ Math.toIntExact(Math.round(ordinates[1])),
+ modelTiePoints.doubleValue(i+3),
+ modelTiePoints.doubleValue(i+4));
+ }
+ grid.setDesiredPrecision(PRECISION);
+ 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. 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 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;
+ }
+ 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;
+ }
+ }
+ /*
+ * 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);
+ }
+ }
+ return MathTransforms.specialize(global, specialization);
+ }
+ }
+
+ /**
+ * 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
+ * a regular grid (e.g. Sentinel 1 image). This method tells where to split in
+ * two grids.
+ *
+ * @param values the x or y vector of tie points pixel coordinates.
+ * @param uniques an initially empty set to be used for this method internal working.
+ * @return value after which a different step is used, or {@code NaN} if none.
+ * The value returned by this method should be included in the first grid.
+ */
+ private static double threshold(final Vector values, final Set<Double> uniques) {
+ final int n = values.size();
+ for (int i=0; i<n; i++) {
+ uniques.add(values.doubleValue(i));
+ }
+ final Double[] array = uniques.toArray(new Double[uniques.size()]);
+ uniques.clear();
+ int i = array.length;
+ if (i >= 3) {
+ Arrays.sort(array);
+ double value;
+ final double inc = array[--i] - (value = array[--i]);
+ do {
+ final double lower = array[--i];
+ if (value - lower != inc) {
+ return value;
+ }
+ value = lower;
+ } while (i > 0);
}
- grid.setDesiredPrecision(PRECISION);
- return grid.create(null);
+ return Double.NaN;
}
/**