Examples with LinearTransform org.apache.sis.referencing.operation.transform.LinearTransform used on opensource projects

Example 1 with LinearTransform

use of org.apache.sis.referencing.operation.transform.LinearTransform 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);
    }
}

Example 2 with LinearTransform

use of org.apache.sis.referencing.operation.transform.LinearTransform in project sis by apache.

the class GeocentricAffine method asDatumShift.

/**
 * Given a transformation chain, conditionally replaces the affine transform elements by an alternative object
 * showing the Bursa-Wolf parameters. The replacement is applied if and only if the affine transform is a scale,
 * translation or rotation in the geocentric domain.
 *
 * <p>This method is invoked only by {@code ConcatenatedTransform.getPseudoSteps()} for the need of WKT formatting.
 * The purpose of this method is very similar to the purpose of {@code AbstractMathTransform.beforeFormat(List, int,
 * boolean)} except that we need to perform the {@code forDatumShift(…)} work only after {@code beforeFormat(…)}
 * finished its work for all {@code ContextualParameters}, including the {@code EllipsoidToCentricTransform}'s one.</p>
 *
 * @param  transforms  the full chain of concatenated transforms.
 */
public static void asDatumShift(final List<Object> transforms) {
    for (int i = transforms.size() - 2; --i >= 0; ) {
        if (isOperation(GeographicToGeocentric.NAME, transforms.get(i)) && isOperation(GeocentricToGeographic.NAME, transforms.get(i + 2))) {
            final Object step = transforms.get(i + 1);
            if (step instanceof LinearTransform) {
                final BursaWolfParameters parameters = new BursaWolfParameters(null, null);
                try {
                    parameters.setPositionVectorTransformation(((LinearTransform) step).getMatrix(), BURSAWOLF_TOLERANCE);
                } catch (IllegalArgumentException e) {
                    /*
                         * Should not occur, except sometime on inverse transform of relatively complex datum shifts
                         * (more than just translation terms). We can fallback on formatting the full matrix.
                         */
                    log(Loggers.WKT, "asDatumShift", e);
                    continue;
                }
                final boolean isTranslation = parameters.isTranslation();
                final Parameters values = createParameters(isTranslation ? GeocentricTranslation.PARAMETERS : PositionVector7Param.PARAMETERS, parameters, isTranslation);
                transforms.set(i + 1, new FormattableObject() {

                    @Override
                    protected String formatTo(final Formatter formatter) {
                        WKTUtilities.appendParamMT(values, formatter);
                        return WKTKeywords.Param_MT;
                    }
                });
            }
        }
    }
}

Example 3 with LinearTransform

use of org.apache.sis.referencing.operation.transform.LinearTransform in project sis by apache.

the class LocalizationGridBuilder method create.

/**
 * 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.
 *
 * @param  factory  the factory to use for creating the transform, or {@code null} for the default factory.
 *                  The {@link MathTransformFactory#createAffineTransform(Matrix)} method of that factory
 *                  shall return {@link LinearTransform} instances.
 * @return the transform from source to target points.
 * @throws FactoryException if the transform can not be created,
 *         for example because the target points have not be specified.
 */
@Override
public MathTransform create(final MathTransformFactory factory) throws FactoryException {
    final LinearTransform gridToCoord = linear.create(factory);
    /*
         * Make a first check about whether the result of above LinearTransformBuilder.create() call
         * can be considered a good fit. If true, then we may return the linear transform directly.
         */
    boolean isExact = true;
    boolean isLinear = true;
    for (final double c : linear.correlation()) {
        isExact &= (c == 1);
        if (c < 0.9999) {
            // Empirical threshold (may need to be revisited).
            isLinear = false;
            break;
        }
    }
    if (isExact) {
        return MathTransforms.concatenate(sourceToGrid, gridToCoord);
    }
    final int width = linear.gridSize(0);
    final int height = linear.gridSize(1);
    final int tgtDim = gridToCoord.getTargetDimensions();
    final double[] residual = new double[tgtDim * linear.gridLength];
    final double[] point = new double[tgtDim + 1];
    double gridPrecision = precision;
    try {
        /*
             * If the user specified a precision, we need to convert it from source units to grid units.
             * We convert each dimension separately, then retain the largest magnitude of vector results.
             */
        if (gridPrecision > 0 && !sourceToGrid.isIdentity()) {
            final double[] vector = new double[sourceToGrid.getSourceDimensions()];
            final double[] offset = new double[sourceToGrid.getTargetDimensions()];
            double converted = 0;
            for (int i = 0; i < vector.length; i++) {
                vector[i] = precision;
                sourceToGrid.deltaTransform(vector, 0, offset, 0, 1);
                final double length = MathFunctions.magnitude(offset);
                if (length > converted)
                    converted = length;
                vector[i] = 0;
            }
            gridPrecision = converted;
        }
        /*
             * Compute the residuals, i.e. the differences between the coordinates that we get by a linear
             * transformation and the coordinates that we want to get. If at least one residual is greater
             * than the desired precision,  then the returned MathTransform will need to apply corrections
             * after linear transforms. Those corrections will be done by InterpolatedTransform.
             */
        final MatrixSIS coordToGrid = MatrixSIS.castOrCopy(gridToCoord.inverse().getMatrix());
        final DirectPosition2D src = new DirectPosition2D();
        point[tgtDim] = 1;
        for (int k = 0, y = 0; y < height; y++) {
            src.y = y;
            tmp[1] = y;
            for (int x = 0; x < width; x++) {
                src.x = x;
                tmp[0] = x;
                // Expected position.
                linear.getControlPoint2D(tmp, point);
                // As grid coordinate.
                double[] grid = coordToGrid.multiply(point);
                isLinear &= (residual[k++] = grid[0] - x) <= gridPrecision;
                isLinear &= (residual[k++] = grid[1] - y) <= gridPrecision;
            }
        }
    } catch (TransformException e) {
        // Should never happen.
        throw new FactoryException(e);
    }
    if (isLinear) {
        return MathTransforms.concatenate(sourceToGrid, gridToCoord);
    }
    return InterpolatedTransform.createGeodeticTransformation(nonNull(factory), new ResidualGrid(sourceToGrid, gridToCoord, width, height, tgtDim, residual, (gridPrecision > 0) ? gridPrecision : DEFAULT_PRECISION));
}

Example 4 with LinearTransform

use of org.apache.sis.referencing.operation.transform.LinearTransform in project sis by apache.

the class DatumShiftGrid method interpolateAt.

/**
 * Interpolates the translation to apply for the given coordinate.
 * The input values are in the unit given by {@link #getCoordinateUnit()}.
 * The output values are in the unit given by {@link #getTranslationUnit()}.
 * The length of the returned array is given by {@link #getTranslationDimensions()}.
 *
 * <div class="section">Default implementation</div>
 * The default implementation performs the following steps:
 * <ol>
 *   <li>Convert the given coordinate into grid indices using the transform given by {@link #getCoordinateToGrid()}.</li>
 *   <li>Interpolate the translation vector at the above grid indices with a call to {@link #interpolateInCell}.</li>
 *   <li>If {@link #isCellValueRatio()} returns {@code true}, {@linkplain LinearTransform#deltaTransform delta transform}
 *       the translation vector by the inverse of the conversion given at step 1.</li>
 * </ol>
 *
 * @param ordinates  the "real world" ordinate (often longitude and latitude, but not necessarily)
 *                   of the point for which to get the translation.
 * @return the translation vector at the given position.
 * @throws TransformException if an error occurred while computing the translation vector.
 */
public double[] interpolateAt(final double... ordinates) throws TransformException {
    final LinearTransform c = getCoordinateToGrid();
    ArgumentChecks.ensureDimensionMatches("ordinates", c.getSourceDimensions(), ordinates);
    final int dim = getTranslationDimensions();
    double[] vector = new double[Math.max(dim, c.getTargetDimensions())];
    c.transform(ordinates, 0, vector, 0, 1);
    interpolateInCell(vector[0], vector[1], vector);
    if (isCellValueRatio()) {
        c.inverse().deltaTransform(vector, 0, vector, 0, 1);
    }
    if (vector.length != dim) {
        vector = Arrays.copyOf(vector, dim);
    }
    return vector;
}

Example 5 with LinearTransform

use of org.apache.sis.referencing.operation.transform.LinearTransform in project sis by apache.

the class GeodeticObjectParserTest method conversion.

/**
 * Returns the conversion from {@code north} to {@code south}.
 */
private static Matrix conversion(final ProjectedCRS north, final ProjectedCRS south) throws NoninvertibleTransformException {
    final MathTransform transform = MathTransforms.concatenate(north.getConversionFromBase().getMathTransform().inverse(), south.getConversionFromBase().getMathTransform());
    assertInstanceOf("North to South", LinearTransform.class, transform);
    return ((LinearTransform) transform).getMatrix();
}