Examples with LocalizationGridBuilder org.apache.sis.referencing.operation.builder.LocalizationGridBuilder used on opensource projects

Example 1 with LocalizationGridBuilder

use of org.apache.sis.referencing.operation.builder.LocalizationGridBuilder in project sis by apache.

the class GridGeometry method localizationGrid.

/**
 * 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, 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 {
        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) {
            // Should never happen. If it does, we have a bug in our algorithm.
            throw new FactoryException();
        }
        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)) {
            // Can not do better. Report the failure.
            throw e;
        }
        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);
            // Number of the part where to add current point.
            int part = 0;
            // Point will be added to part #1 or #3.
            if (px > splitX)
                part = 1;
            // Point will be added to part #2 or #3.
            if (py > splitY)
                part |= 2;
            // Bitmask of the parts where to add the point.
            int parts = 1 << part;
            // Add also the point to part #1 or #3.
            if (px == splitX)
                parts |= 1 << (part | 1);
            // Add also the point to part #2 or #3.
            if (py == splitY)
                parts |= 1 << (part | 2);
            if (parts == 0b0111) {
                // Add also the point to part #3.
                parts = 0b1111;
                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 (// Clear the bit of the part we processed.
            (parts &= ~(1 << part)) != 0);
            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];
            // Safety against infinite recursivity.
            if (length >= size)
                throw e;
            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);
    }
}