Examples with MatrixSIS org.apache.sis.referencing.operation.matrix.MatrixSIS used on opensource projects

Example 6 with MatrixSIS

use of org.apache.sis.referencing.operation.matrix.MatrixSIS in project sis by apache.

the class ContextualParameters method beforeFormat.

/**
 * Given a transformation chain, replaces the elements around {@code transforms.get(index)} transform by
 * alternative objects to use when formatting WKT. The replacement is performed in-place in the given list.
 *
 * <p>This method shall replace only the previous element and the few next elements that need
 * to be changed as a result of the previous change. This method is not expected to continue
 * the iteration after the changes that are of direct concern to this object.</p>
 *
 * <p>This method is invoked (indirectly) only by {@link ConcatenatedTransform#getPseudoSteps()} in order
 * to get the {@link ParameterValueGroup} of a map projection, or to format a {@code ProjectedCRS} WKT.</p>
 *
 * @param  transforms  the full chain of concatenated transforms.
 * @param  index       the index of this transform in the {@code transforms} chain.
 * @param  inverse     always {@code false}, except if we are formatting the inverse transform.
 * @return index of this transform in the {@code transforms} chain after processing.
 *
 * @see ConcatenatedTransform#getPseudoSteps()
 * @see AbstractMathTransform#beforeFormat(List, int, boolean)
 */
final int beforeFormat(final List<Object> transforms, int index, final boolean inverse) {
    /*
         * We expect affine transforms before and after the normalized projection. Extracts those
         * affine transforms now. If one or both are missing, we will treat null as an identity
         * transform. We will not replace the elements in the list before new values for those
         * affine transforms have been fully calculated.
         */
    Matrix before = null;
    Matrix after = null;
    if (index != 0) {
        final Object candidate = transforms.get(index - 1);
        if (candidate instanceof MathTransform) {
            before = MathTransforms.getMatrix((MathTransform) candidate);
        }
    }
    if (index + 1 < transforms.size()) {
        final Object candidate = transforms.get(index + 1);
        if (candidate instanceof MathTransform) {
            after = MathTransforms.getMatrix((MathTransform) candidate);
        }
    }
    final boolean hasBefore = (before != null);
    final boolean hasAfter = (after != null);
    /*
         * We assume that the "before" affine contains the normalize operation to be applied
         * before the projection. However it may contains more than just this normalization,
         * because it may have been concatenated with any user-defined transform (for example
         * in order to apply a change of axis order). We need to separate the "user-defined"
         * step from the "normalize" step.
         */
    MatrixSIS userDefined;
    try {
        userDefined = getMatrix(inverse ? MatrixRole.DENORMALIZATION : MatrixRole.INVERSE_NORMALIZATION);
    } catch (IllegalStateException e) {
        /*
             * Should never happen. But if it does, we abandon the attempt to change
             * the list elements and will format the objects in their "raw" format.
             */
        unexpectedException(e);
        return index;
    }
    if (hasBefore) {
        userDefined = userDefined.multiply(before);
    }
    /*
         * At this point "userDefined" is the affine transform to show to user instead of the
         * "before" affine transform. Replaces "before" by "userDefined" locally (but not yet
         * in the list), or set it to null (meaning that it will be removed from the list) if
         * it is identity, which happen quite often.
         *
         * Note on rounding error: the coefficients are often either 0 or 1 since the transform
         * is often for changing axis order. Thanks to double-double arithmetic in SIS matrices,
         * the non-zero values are usually accurate. But the values that should be zero are much
         * harder to get right. Sometime we see small values (around 1E-12) in the last column of
         * the 'before' matrix below. Since this column contains translation terms, those numbers
         * are in the unit of measurement of input values of the MathTransform after the matrix.
         *
         *   - For forward map projections, those values are conceptually in decimal degrees
         *     (in fact the values are converted to radians but not by this 'before' matrix).
         *
         *   - For inverse map projections, those values are conceptually in metres (in fact
         *     converted to distances on a unitary ellipsoid but not by this 'before' matrix).
         *
         *   - Geographic/Geocentric transformations behave like map projections in regard to units.
         *     Molodensky transformations conceptually use always decimal degrees. There is not much
         *     other cases since this mechanism is internal to SIS (not in public API).
         *
         * Consequently we set the tolerance threshold to ANGULAR_TOLERANCE. We do not bother (at least
         * for now) to identify the cases where we could use LINEAR_TOLERANCE because just checking the
         * 'inverse' flag is not sufficient (e.g. the Molodensky case). Since the angular tolerance is
         * smaller than the linear one, unconditional usage of ANGULAR_TOLERANCE is more conservative.
         */
    before = Matrices.isIdentity(userDefined, Formulas.ANGULAR_TOLERANCE) ? null : userDefined;
    /*
         * Compute the "after" affine transform in a way similar than the "before" affine.
         * Note that if this operation fails, we will cancel everything we would have done
         * in this method (i.e. we do not touch the transforms list at all).
         */
    try {
        userDefined = getMatrix(inverse ? MatrixRole.NORMALIZATION : MatrixRole.INVERSE_DENORMALIZATION);
    } catch (IllegalStateException e) {
        unexpectedException(e);
        return index;
    }
    if (hasAfter) {
        userDefined = Matrices.multiply(after, userDefined);
    }
    /*
         * Note on rounding error: same discussion than the "note on rounding error" of the 'before' matrix,
         * with the following differences:
         *
         *   - For forward map projections, unit of measurements of translation terms are conceptually
         *     metres (instead than degrees) multiplied by the scale factors in the 'after' matrix.
         *
         *   - For inverse map projections, unit of measurements of translation terms are conceptually
         *     degrees (instead than metres) multiplied by the scale factors in the 'after' matrix.
         *
         *   - And so on for all cases: swap the units of the forward and inverse cases, then multiply
         *     by the scale factor. Note that the multiplication step does not exist in the 'before' case.
         *
         * Since we are seeking for the identity matrix, the scale factor is 1. We do not bother to distinguish
         * the ANGULAR_TOLERANCE and LINEAR_TOLERANCE cases for the same reasons than for the 'before' matrix.
         */
    after = Matrices.isIdentity(userDefined, Formulas.ANGULAR_TOLERANCE) ? null : userDefined;
    /*
         * At this point we have computed all the affine transforms to show to the user.
         * We can replace the elements in the list. The transform referenced by transforms.get(index)
         * is usually a NormalizedProjection, to be replaced by a ContextualParameters instance in order
         * to format real parameter values (semi-major axis, scale factor, etc.)
         * instead than a semi-major axis length of 1.
         */
    if (before == null) {
        if (hasBefore) {
            final Object old = transforms.remove(--index);
            assert (old instanceof LinearTransform);
        }
    } else {
        if (hasBefore) {
            final Object old = transforms.set(index - 1, before);
            assert (old instanceof LinearTransform);
        } else {
            transforms.add(index++, before);
        }
    }
    transforms.set(index, new WKT(inverse));
    if (after == null) {
        if (hasAfter) {
            final Object old = transforms.remove(index + 1);
            assert (old instanceof LinearTransform);
        }
    } else {
        if (hasAfter) {
            final Object old = transforms.set(index + 1, after);
            assert (old instanceof LinearTransform);
        } else {
            transforms.add(index + 1, after);
        }
    }
    return index;
}

Example 7 with MatrixSIS

use of org.apache.sis.referencing.operation.matrix.MatrixSIS in project sis by apache.

the class PassThroughTransform method toSubMatrix.

/**
 * If the given matrix to be concatenated to this transform, can be concatenated to the
 * sub-transform instead, returns the matrix to be concatenated to the sub-transform.
 * Otherwise returns {@code null}.
 *
 * <p>This method assumes that the matrix size is compatible with this transform source dimension.
 * It is caller responsibility to verify this condition.</p>
 */
final Matrix toSubMatrix(final Matrix matrix) {
    final int numRow = matrix.getNumRow();
    final int numCol = matrix.getNumCol();
    if (numRow != numCol) {
        // Current implementation requires a square matrix.
        return null;
    }
    final int subDim = subTransform.getSourceDimensions();
    final MatrixSIS sub = Matrices.createIdentity(subDim + 1);
    /*
         * Ensure that every dimensions which are scaled by the affine transform are one
         * of the dimensions modified by the sub-transform, and not any other dimension.
         */
    for (int j = numRow; --j >= 0; ) {
        final int sj = j - firstAffectedOrdinate;
        for (int i = numCol; --i >= 0; ) {
            final double element = matrix.getElement(j, i);
            if (sj >= 0 && sj < subDim) {
                final int si;
                final boolean pass;
                if (i == numCol - 1) {
                    // Translation term (last column)
                    si = subDim;
                    pass = true;
                } else {
                    // Any term other than translation.
                    si = i - firstAffectedOrdinate;
                    pass = (si >= 0 && si < subDim);
                }
                if (pass) {
                    sub.setElement(sj, si, element);
                    continue;
                }
            }
            if (element != (i == j ? 1 : 0)) {
                // Found a dimension which perform some scaling or translation.
                return null;
            }
        }
    }
    return sub;
}

Example 8 with MatrixSIS

use of org.apache.sis.referencing.operation.matrix.MatrixSIS in project sis by apache.

the class ScaleTransform method derivative.

/**
 * Gets the derivative of this transform at a point.
 * For a matrix transform, the derivative is the same everywhere.
 *
 * @param  point  ignored (can be {@code null}).
 */
@Override
public Matrix derivative(final DirectPosition point) {
    final int n = factors.length;
    final MatrixSIS matrix = Matrices.createZero(n, n + numDroppedDimensions);
    for (int i = 0; i < n; i++) {
        matrix.setElement(i, i, factors[i]);
    }
    return matrix;
}

Example 9 with MatrixSIS

use of org.apache.sis.referencing.operation.matrix.MatrixSIS in project sis by apache.

the class MathTransformContext method getMatrix.

/**
 * Returns the normalization or denormalization matrix.
 */
@Override
@SuppressWarnings("fallthrough")
public Matrix getMatrix(final MatrixRole role) throws FactoryException {
    final CoordinateSystem cs;
    boolean inverse = false;
    double rotation;
    switch(role) {
        default:
            throw new IllegalArgumentException(Errors.format(Errors.Keys.IllegalArgumentValue_2, "role", role));
        // Fall through
        case INVERSE_NORMALIZATION:
            inverse = true;
        case NORMALIZATION:
            rotation = sourceMeridian;
            cs = getSourceCS();
            break;
        // Fall through
        case INVERSE_DENORMALIZATION:
            inverse = true;
        case DENORMALIZATION:
            inverse = !inverse;
            rotation = targetMeridian;
            cs = getTargetCS();
            break;
    }
    Matrix matrix = super.getMatrix(role);
    if (rotation != 0) {
        if (inverse)
            rotation = -rotation;
        MatrixSIS cm = MatrixSIS.castOrCopy(matrix);
        if (cs instanceof CartesianCS) {
            rotation = Math.toRadians(rotation);
            final Matrix4 rot = new Matrix4();
            rot.m00 = rot.m11 = Math.cos(rotation);
            rot.m01 = -(rot.m10 = Math.sin(rotation));
            if (inverse) {
                // Apply the rotation after denormalization.
                matrix = Matrices.multiply(rot, cm);
            } else {
                // Apply the rotation before normalization.
                matrix = cm.multiply(rot);
            }
        } else if (cs == null || cs instanceof EllipsoidalCS || cs instanceof SphericalCS) {
            final Double value = rotation;
            if (inverse) {
                // Longitude is the first axis in normalized CS.
                cm.convertBefore(0, null, value);
            } else {
                cm.convertAfter(0, null, value);
            }
            matrix = cm;
        } else {
            throw new FactoryException(Errors.format(Errors.Keys.UnsupportedCoordinateSystem_1, cs.getName()));
        }
    }
    return matrix;
}

Example 10 with MatrixSIS

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