Examples with Matrix org.opengis.referencing.operation.Matrix used on opensource projects

Example 26 with Matrix

use of org.opengis.referencing.operation.Matrix in project sis by apache.

the class ConcatenatedTransformDirect2D method derivative.

/**
 * Gets the derivative of this transform at a point.
 *
 * @param  point  the coordinate point where to evaluate the derivative.
 * @return the derivative at the specified point (never {@code null}).
 * @throws TransformException if the derivative can't be evaluated at the specified point.
 */
@Override
public Matrix derivative(final Point2D point) throws TransformException {
    final MathTransform2D transform1 = (MathTransform2D) this.transform1;
    final MathTransform2D transform2 = (MathTransform2D) this.transform2;
    final Matrix matrix1 = transform1.derivative(point);
    final Matrix matrix2 = transform2.derivative(transform1.transform(point, null));
    return Matrices.multiply(matrix2, matrix1);
}

Example 27 with Matrix

use of org.opengis.referencing.operation.Matrix 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 28 with Matrix

use of org.opengis.referencing.operation.Matrix in project sis by apache.

the class ContextualParameters method getMatrix.

/**
 * Returns the affine transforms to be applied before or after the non-linear kernel operation.
 * Immediately after {@linkplain #ContextualParameters(OperationMethod) construction}, those matrices
 * are modifiable identity matrices. Callers can modify the matrix element values, typically by calls to
 * the {@link MatrixSIS#convertBefore(int, Number, Number) MatrixSIS.convertBefore(…)} method.
 * Alternatively, the following methods can be invoked for applying some frequently used configurations:
 *
 * <ul>
 *   <li>{@link #normalizeGeographicInputs(double)}</li>
 *   <li>{@link #denormalizeGeographicOutputs(double)}</li>
 * </ul>
 *
 * After the {@link #completeTransform(MathTransformFactory, MathTransform) completeTransform(…)} method has been
 * invoked, the matrices returned by this method are {@linkplain Matrices#unmodifiable(Matrix) unmodifiable}.
 *
 * <div class="note"><b>Application to map projections:</b>
 * after {@link org.apache.sis.referencing.operation.projection.NormalizedProjection} construction, the matrices
 * returned by {@code projection.getContextualParameters().getMatrix(…)} are initialized to the values shown below.
 * Note that some {@code NormalizedProjection} subclasses apply further modifications to those matrices.
 *
 * <table class="sis">
 *   <caption>Initial matrix coefficients after {@code NormalizedProjection} construction</caption>
 *   <tr>
 *     <th>{@code getMatrix(NORMALIZATION)}</th>
 *     <th class="sep">{@code getMatrix(DENORMALIZATION)}</th>
 *   </tr><tr>
 *     <td>{@include formulas.html#NormalizeGeographic}</td>
 *     <td class="sep">{@include formulas.html#DenormalizeCartesian}</td>
 *   </tr>
 * </table>
 * </div>
 *
 * @param  role  {@code NORMALIZATION} for fetching the <cite>normalization</cite> transform to apply before the kernel,
 *               {@code DENORMALIZATION} for the <cite>denormalization</cite> transform to apply after the kernel, or
 *               {@code INVERSE_*} for the inverse of the above-cited matrices.
 * @return the matrix for the requested normalization or denormalization affine transform.
 *
 * @since 0.7
 */
public final MatrixSIS getMatrix(MatrixRole role) {
    final Matrix fallback;
    final ContextualParameters inverse;
    synchronized (this) {
        switch(role) {
            default:
                throw new AssertionError(role);
            case NORMALIZATION:
                return toMatrixSIS(normalize);
            case DENORMALIZATION:
                return toMatrixSIS(denormalize);
            case INVERSE_NORMALIZATION:
                role = MatrixRole.DENORMALIZATION;
                fallback = normalize;
                break;
            case INVERSE_DENORMALIZATION:
                role = MatrixRole.NORMALIZATION;
                fallback = denormalize;
                break;
        }
        // Copy the reference while we are inside the synchronized block.
        inverse = this.inverse;
    }
    /*
         * Following must be outside the synchronized block in order to avoid potential deadlock while invoking
         * inverse.getMatrix(role). We do not cache the matrix here, but 'inverse' is likely to have cached it.
         */
    final Matrix m;
    if (inverse != null) {
        m = inverse.getMatrix(role);
    } else
        try {
            m = Matrices.inverse(fallback);
        } catch (NoninvertibleMatrixException e) {
            throw new IllegalStateException(Errors.format(Errors.Keys.CanNotCompute_1, role), e);
        }
    return Matrices.unmodifiable(m);
}

Example 29 with Matrix

use of org.opengis.referencing.operation.Matrix in project sis by apache.

the class DefaultMathTransformFactory method swapAndScaleAxes.

/**
 * Given a transform between normalized spaces,
 * creates a transform taking in account axis directions, units of measurement and longitude rotation.
 * This method {@linkplain #createConcatenatedTransform concatenates} the given parameterized transform
 * with any other transform required for performing units changes and ordinates swapping.
 *
 * <p>The given {@code parameterized} transform shall expect
 * {@linkplain org.apache.sis.referencing.cs.AxesConvention#NORMALIZED normalized} input coordinates and
 * produce normalized output coordinates. See {@link org.apache.sis.referencing.cs.AxesConvention} for more
 * information about what Apache SIS means by "normalized".</p>
 *
 * <div class="note"><b>Example:</b>
 * The most typical examples of transforms with normalized inputs/outputs are normalized
 * map projections expecting (<cite>longitude</cite>, <cite>latitude</cite>) inputs in degrees
 * and calculating (<cite>x</cite>, <cite>y</cite>) coordinates in metres,
 * both of them with ({@linkplain org.opengis.referencing.cs.AxisDirection#EAST East},
 * {@linkplain org.opengis.referencing.cs.AxisDirection#NORTH North}) axis orientations.</div>
 *
 * <div class="section">Controlling the normalization process</div>
 * Users who need a different normalized space than the default one way find more convenient to
 * override the {@link Context#getMatrix Context.getMatrix(ContextualParameters.MatrixRole)} method.
 *
 * @param  parameterized  a transform for normalized input and output coordinates.
 * @param  context        source and target coordinate systems in which the transform is going to be used.
 * @return a transform taking in account unit conversions and axis swapping.
 * @throws FactoryException if the object creation failed.
 *
 * @see org.apache.sis.referencing.cs.AxesConvention#NORMALIZED
 * @see org.apache.sis.referencing.operation.DefaultConversion#DefaultConversion(Map, OperationMethod, MathTransform, ParameterValueGroup)
 *
 * @since 0.7
 */
public MathTransform swapAndScaleAxes(final MathTransform parameterized, final Context context) throws FactoryException {
    ArgumentChecks.ensureNonNull("parameterized", parameterized);
    ArgumentChecks.ensureNonNull("context", context);
    /*
         * Computes matrix for swapping axis and performing units conversion.
         * There is one matrix to apply before projection on (longitude,latitude)
         * coordinates, and one matrix to apply after projection on (easting,northing)
         * coordinates.
         */
    final Matrix swap1 = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
    final Matrix swap3 = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
    /*
         * Prepares the concatenation of the matrices computed above and the projection.
         * Note that at this stage, the dimensions between each step may not be compatible.
         * For example the projection (step2) is usually two-dimensional while the source
         * coordinate system (step1) may be three-dimensional if it has a height.
         */
    MathTransform step1 = (swap1 != null) ? createAffineTransform(swap1) : MathTransforms.identity(parameterized.getSourceDimensions());
    MathTransform step3 = (swap3 != null) ? createAffineTransform(swap3) : MathTransforms.identity(parameterized.getTargetDimensions());
    MathTransform step2 = parameterized;
    /*
         * Special case for the way EPSG handles reversal of axis direction. For now the "Vertical Offset" (EPSG:9616)
         * method is the only one for which we found a need for special case. But if more special cases are added in a
         * future SIS version, then we should replace the static method by a non-static one defined in AbstractProvider.
         */
    if (context.provider instanceof VerticalOffset) {
        step2 = VerticalOffset.postCreate(step2, swap3);
    }
    /*
         * If the target coordinate system has a height, instructs the projection to pass
         * the height unchanged from the base CRS to the target CRS. After this block, the
         * dimensions of 'step2' and 'step3' should match.
         */
    final int numTrailingOrdinates = step3.getSourceDimensions() - step2.getTargetDimensions();
    if (numTrailingOrdinates > 0) {
        step2 = createPassThroughTransform(0, step2, numTrailingOrdinates);
    }
    /*
         * If the source CS has a height but the target CS doesn't, drops the extra coordinates.
         * After this block, the dimensions of 'step1' and 'step2' should match.
         */
    final int sourceDim = step1.getTargetDimensions();
    final int targetDim = step2.getSourceDimensions();
    if (sourceDim > targetDim) {
        final Matrix drop = Matrices.createDiagonal(targetDim + 1, sourceDim + 1);
        // Element in the lower-right corner.
        drop.setElement(targetDim, sourceDim, 1);
        step1 = createConcatenatedTransform(createAffineTransform(drop), step1);
    }
    MathTransform mt = createConcatenatedTransform(createConcatenatedTransform(step1, step2), step3);
    /*
         * At this point we finished to create the transform.  But before to return it, verify if the
         * parameterized transform given in argument had some custom parameters. This happen with the
         * Equirectangular projection, which can be simplified as an AffineTransform while we want to
         * continue to describe it with the "semi_major", "semi_minor", etc. parameters  instead than
         * "elt_0_0", "elt_0_1", etc.  The following code just forwards those parameters to the newly
         * created transform; it does not change the operation.
         */
    if (parameterized instanceof ParameterizedAffine && !(mt instanceof ParameterizedAffine)) {
        mt = ((ParameterizedAffine) parameterized).newTransform(mt);
    }
    return mt;
}

Example 30 with Matrix

use of org.opengis.referencing.operation.Matrix in project sis by apache.

the class TransferFunction method setLinearTerms.

/**
 * Sets the {@link #scale} and {@link #offset} terms from the given function.
 *
 * @param  function  the transform to set.
 * @throws IllegalArgumentException if this method does not recognize the given transform.
 */
private void setLinearTerms(final LinearTransform function) throws IllegalArgumentException {
    final Matrix m = function.getMatrix();
    final int numRow = m.getNumRow();
    final int numCol = m.getNumCol();
    if (numRow != 2 || numCol != 2) {
        final Integer two = 2;
        throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedMatrixSize_4, two, two, numRow, numCol));
    }
    scale = m.getElement(0, 0);
    offset = m.getElement(0, 1);
}