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

Example 31 with Matrix

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

the class TransformSeparator method filterTargetDimensions.

/**
 * Creates a transform for the same mathematic than the given {@code step}
 * but producing only the given dimensions as outputs.
 * This method is invoked by {@link #separate()} when user-specified target dimensions need to be taken in account.
 * The given {@code step} and {@code dimensions} are typically the values of
 * {@link #transform} and {@link #targetDimensions} fields respectively, but not necessarily.
 *
 * <p>Subclasses can override this method if they need to handle some {@code MathTransform} implementations
 * in a special way. However all implementations of this method shall obey to the following contract:</p>
 * <ul>
 *   <li>{@link #sourceDimensions} and {@link #targetDimensions} should not be assumed accurate.</li>
 *   <li>{@link #sourceDimensions} should not be modified by this method.</li>
 *   <li>{@link #targetDimensions} should not be modified by this method.</li>
 * </ul>
 *
 * The number and nature of inputs stay unchanged. For example if the supplied {@code transform}
 * has (<var>longitude</var>, <var>latitude</var>, <var>height</var>) outputs, then a filtered
 * transform may keep only the (<var>longitude</var>, <var>latitude</var>) part for the same inputs.
 * In most cases, the filtered transform is non-invertible since it loose informations.
 *
 * @param  step        the transform for which to retain only a subset of the target dimensions.
 * @param  dimensions  indices of the target dimensions of {@code step} to retain.
 * @return a transform producing only the given target dimensions.
 * @throws FactoryException if the given transform is not separable.
 */
protected MathTransform filterTargetDimensions(MathTransform step, final int[] dimensions) throws FactoryException {
    final int numSrc = step.getSourceDimensions();
    int numTgt = step.getTargetDimensions();
    final int lower = dimensions[0];
    final int upper = dimensions[dimensions.length - 1];
    if (lower == 0 && upper == numTgt && dimensions.length == numTgt) {
        return step;
    }
    /*
         * If the transform is an instance of passthrough transform but no dimension from its sub-transform
         * is requested, then ignore the sub-transform (i.e. treat the whole transform as identity, except
         * for the number of target dimension which may be different from the number of input dimension).
         */
    int removeAt = 0;
    int numRemoved = 0;
    if (step instanceof PassThroughTransform) {
        final PassThroughTransform passThrough = (PassThroughTransform) step;
        final int subLower = passThrough.firstAffectedOrdinate;
        final int numSubTgt = passThrough.subTransform.getTargetDimensions();
        if (!containsAny(dimensions, subLower, subLower + numSubTgt)) {
            step = IdentityTransform.create(numTgt = numSrc);
            removeAt = subLower;
            numRemoved = numSubTgt - passThrough.subTransform.getSourceDimensions();
        }
    }
    /*                                                  ┌  ┐     ┌          ┐ ┌ ┐
         * Create the matrix to be used as a filter         │x'│     │1  0  0  0│ │x│
         * and concatenate it to the transform. The         │z'│  =  │0  0  1  0│ │y│
         * matrix will contain 1 only in the target         │1 │     │0  0  0  1│ │z│
         * dimensions to keep, as in this example:          └  ┘     └          ┘ │1│
         *                                                                        └ ┘
         */
    final Matrix matrix = Matrices.createZero(dimensions.length + 1, numTgt + 1);
    for (int j = 0; j < dimensions.length; j++) {
        int i = dimensions[j];
        if (i >= removeAt) {
            i -= numRemoved;
        }
        matrix.setElement(j, i, 1);
    }
    matrix.setElement(dimensions.length, numTgt, 1);
    return factory.createConcatenatedTransform(step, factory.createAffineTransform(matrix));
}

Example 32 with Matrix

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

the class AbstractMathTransform2D method derivative.

/**
 * Implementation of {@link #derivative(DirectPosition)} shared by the inverse transform.
 */
static Matrix derivative(final AbstractMathTransform tr, final Point2D point) throws TransformException {
    final double[] coordinate = new double[] { point.getX(), point.getY() };
    final Matrix derivative = tr.transform(coordinate, 0, null, 0, true);
    if (derivative == null) {
        throw new TransformException(Resources.format(Resources.Keys.CanNotComputeDerivative));
    }
    return derivative;
}

Example 33 with Matrix

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

the class ConcatenatedTransform method createOptimized.

/**
 * Tries to returns an optimized concatenation, for example by merging two affine transforms
 * into a single one. If no optimized cases has been found, returns {@code null}. In the later
 * case, the caller will need to create a more heavy {@link ConcatenatedTransform} instance.
 *
 * @param  factory  the factory which is (indirectly) invoking this method, or {@code null} if none.
 */
private static MathTransform createOptimized(final MathTransform tr1, final MathTransform tr2, final MathTransformFactory factory) throws FactoryException {
    /*
         * Trivial - but actually essential!! - check for the identity cases.
         */
    if (tr1.isIdentity())
        return tr2;
    if (tr2.isIdentity())
        return tr1;
    /*
         * If both transforms use matrix, then we can create
         * a single transform using the concatenated matrix.
         */
    final Matrix matrix1 = MathTransforms.getMatrix(tr1);
    if (matrix1 != null) {
        final Matrix matrix2 = MathTransforms.getMatrix(tr2);
        if (matrix2 != null) {
            final Matrix matrix = Matrices.multiply(matrix2, matrix1);
            if (Matrices.isIdentity(matrix, IDENTITY_TOLERANCE)) {
                // Returns a cached instance.
                return MathTransforms.identity(matrix.getNumRow() - 1);
            }
            /*
                 * NOTE: It is quite tempting to "fix rounding errors" in the matrix before to create the transform.
                 * But this is often wrong for datum shift transformations (Molodensky and the like) since the datum
                 * shifts are very small. The shift may be the order of magnitude of the tolerance threshold. Intead,
                 * Apache SIS performs matrix operations using double-double arithmetic in the hope to get exact
                 * results at the 'double' accuracy, which avoid the need for a tolerance threshold.
                 */
            if (factory != null) {
                return factory.createAffineTransform(matrix);
            } else {
                return MathTransforms.linear(matrix);
            }
        }
        /*
             * If the second transform is a passthrough transform and all passthrough ordinates
             * are unchanged by the matrix, we can move the matrix inside the passthrough transform.
             */
        if (tr2 instanceof PassThroughTransform) {
            final PassThroughTransform candidate = (PassThroughTransform) tr2;
            final Matrix sub = candidate.toSubMatrix(matrix1);
            if (sub != null) {
                if (factory != null) {
                    return factory.createPassThroughTransform(candidate.firstAffectedOrdinate, factory.createConcatenatedTransform(factory.createAffineTransform(sub), candidate.subTransform), candidate.numTrailingOrdinates);
                } else {
                    return PassThroughTransform.create(candidate.firstAffectedOrdinate, create(MathTransforms.linear(sub), candidate.subTransform, factory), candidate.numTrailingOrdinates);
                }
            }
        }
    }
    /*
         * If one transform is the inverse of the other, return the identity transform.
         */
    if (areInverse(tr1, tr2) || areInverse(tr2, tr1)) {
        assert tr1.getSourceDimensions() == tr2.getTargetDimensions();
        assert tr1.getTargetDimensions() == tr2.getSourceDimensions();
        // Returns a cached instance.
        return MathTransforms.identity(tr1.getSourceDimensions());
    }
    /*
         * Give a chance to AbstractMathTransform to returns an optimized object.
         * The main use case is Logarithmic vs Exponential transforms.
         */
    if (tr1 instanceof AbstractMathTransform) {
        final MathTransform optimized = ((AbstractMathTransform) tr1).tryConcatenate(false, tr2, factory);
        if (optimized != null) {
            return optimized;
        }
    }
    if (tr2 instanceof AbstractMathTransform) {
        final MathTransform optimized = ((AbstractMathTransform) tr2).tryConcatenate(true, tr1, factory);
        if (optimized != null) {
            return optimized;
        }
    }
    // No optimized case found.
    return null;
}

Example 34 with Matrix

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

the class ConcatenatedTransform method getParameterised.

/**
 * If there is exactly one transform step which is {@linkplain Parameterized parameterized},
 * returns that transform step. Otherwise returns {@code null}.
 *
 * <p>This method normally requires that there is exactly one transform step remaining after we
 * processed map projections in the special way described in {@link #getParameterValues()},
 * because if they were more than one remaining steps, the returned parameters would not be
 * sufficient for rebuilding the full concatenated transform. Returning parameters when there
 * is more than one remaining step, even if all other transform steps are not parameterizable,
 * would be a contract violation.</p>
 *
 * <p>However in the special case where we are getting the parameters of a {@code CoordinateOperation} instance
 * through {@link org.apache.sis.referencing.operation.AbstractCoordinateOperation#getParameterValues()} method
 * (often indirectly trough WKT formatting of a {@code "ProjectedCRS"} element), then the above rule is slightly
 * relaxed: we ignore affine transforms in order to accept axis swapping or unit conversions. We do that in that
 * particular case only because the coordinate systems given with the enclosing {@code CoordinateOperation} or
 * {@code GeneralDerivedCRS} specify the axis swapping and unit conversions.
 * This special case is internal to SIS implementation and should be unknown to users.</p>
 *
 * @return the parameterizable transform step, or {@code null} if none.
 */
private Parameterized getParameterised() {
    Parameterized param = null;
    final List<Object> transforms = getPseudoSteps();
    if (transforms.size() == 1 || Semaphores.query(Semaphores.ENCLOSED_IN_OPERATION)) {
        for (final Object candidate : transforms) {
            /*
                 * Search for non-linear parameters only, ignoring affine transforms and the matrices
                 * computed by ContextualParameters. Note that the 'transforms' list is guaranteed to
                 * contains at least one non-linear parameter, otherwise we would not have created a
                 * ConcatenatedTransform instance.
                 */
            if (!(candidate instanceof Matrix) && !(candidate instanceof LinearTransform)) {
                if ((param == null) && (candidate instanceof Parameterized)) {
                    param = (Parameterized) candidate;
                } else {
                    /*
                         * Found more than one group of non-linear parameters, or found an object
                         * that do not declare its parameters.  In the later case, conservatively
                         * returns 'null' because we do not know what the real parameters are.
                         */
                    return null;
                }
            }
        }
    }
    return param;
}

Example 35 with Matrix

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

the class MolodenskyFormula method transform.

/**
 * Implementation of {@link #transform(double[], int, double[], int, boolean)} with possibility
 * to override some field values. In this method signature, parameters having the same name than
 * fields have the same value, except in some special circumstances:
 *
 * <ul>
 *   <li>{@code tX}, {@code tY} and {@code tZ} parameters always have the values of {@link #tX}, {@link #tY}
 *       and {@link #tZ} fields when this method is invoked by {@link MolodenskyTransform}. But those values
 *       may be slightly different when this method is invoked by {@link InterpolatedMolodenskyTransform}.</li>
 * </ul>
 *
 * @param  λ           longitude (radians).
 * @param  φ           latitude (radians).
 * @param  h           height above the ellipsoid in unit of semi-major axis.
 * @param  dstPts      the array into which the transformed coordinate is returned, or {@code null}.
 * @param  dstOff      the offset to the location of the transformed point that is stored in the destination array.
 * @param  tX          the {@link #tX} field value (or a slightly different value during geocentric interpolation).
 * @param  tY          the {@link #tY} field value (or a slightly different value during geocentric interpolation).
 * @param  tZ          the {@link #tZ} field value (or a slightly different value during geocentric interpolation).
 * @param  offset      an array of length 3 if this method should use the interpolation grid, or {@code null} otherwise.
 * @param  derivate    {@code true} for computing the derivative, or {@code false} if not needed.
 * @throws TransformException if a point can not be transformed.
 */
final Matrix transform(final double λ, final double φ, final double h, final double[] dstPts, int dstOff, double tX, double tY, double tZ, double[] offset, final boolean derivate) throws TransformException {
    /*
         * Abridged Molodensky formulas from EPSG guidance note:
         *
         *     ν   = a / √(1 - ℯ²⋅sin²φ)                        : radius of curvature in the prime vertical
         *     ρ   = a⋅(1 – ℯ²) / (1 – ℯ²⋅sin²φ)^(3/2)          : radius of curvature in the meridian
         *     Δλ″ = (-tX⋅sinλ + tY⋅cosλ) / (ν⋅cosφ⋅sin1″)
         *     Δφ″ = (-tX⋅sinφ⋅cosλ - tY⋅sinφ⋅sinλ + tZ⋅cosφ + [a⋅Δf + f⋅Δa]⋅sin(2φ)) / (ρ⋅sin1″)
         *     Δh  = tX⋅cosφ⋅cosλ + tY⋅cosφ⋅sinλ + tZ⋅sinφ + (a⋅Δf + f⋅Δa)⋅sin²φ - Δa
         *
         * we set:
         *
         *    dfm     = (a⋅Δf + f⋅Δa) in abridged case (b⋅Δf in non-abridged case)
         *    sin(2φ) = 2⋅sin(φ)⋅cos(φ)
         */
    final double sinλ = sin(λ);
    final double cosλ = cos(λ);
    final double sinφ = sin(φ);
    final double cosφ = cos(φ);
    final double sin2φ = sinφ * sinφ;
    // Square of the denominator of ν
    final double ν2den = 1 - eccentricitySquared * sin2φ;
    // Denominator of ν
    final double νden = sqrt(ν2den);
    // Denominator of ρ
    final double ρden = ν2den * νden;
    // Other notation: Rm = ρ
    double ρ = semiMajor * (1 - eccentricitySquared) / ρden;
    // Other notation: Rn = ν
    double ν = semiMajor / νden;
    // A term in the calculation of Δφ
    double t = Δfmod * 2;
    if (!isAbridged) {
        ρ += h;
        ν += h;
        t = // = Δf⋅[ν⋅(b/a) + ρ⋅(a/b)]     (without the +h in ν and ρ)
        t * (0.5 / νden + 0.5 / ρden) + // = Δa⋅[ℯ²⋅ν/a]
        Δa * eccentricitySquared / νden;
    }
    // Intermediate terms to be reused by the derivative
    double spcλ, cmsλ, cmsφ, scaleX, scaleY;
    // The target geographic coordinates
    double λt, φt;
    do {
        // "spc" stands for "sin plus cos"
        spcλ = tY * sinλ + tX * cosλ;
        // "cms" stands for "cos minus sin"
        cmsλ = tY * cosλ - tX * sinλ;
        cmsφ = (tZ + t * sinφ) * cosφ - spcλ * sinφ;
        scaleX = ANGULAR_SCALE / (ν * cosφ);
        scaleY = ANGULAR_SCALE / ρ;
        λt = λ + (cmsλ * scaleX);
        φt = φ + (cmsφ * scaleY);
        if (offset == null)
            break;
        /*
             * Following is executed only in InterpolatedMolodenskyTransform case.
             */
        grid.interpolateInCell(grid.normalizedToGridX(λt), grid.normalizedToGridY(φt), offset);
        tX = -offset[0];
        tY = -offset[1];
        tZ = -offset[2];
        offset = null;
    } while (true);
    if (dstPts != null) {
        dstPts[dstOff++] = λt;
        dstPts[dstOff++] = φt;
        if (isTarget3D) {
            // A term in the calculation of Δh
            double t1 = Δfmod * sin2φ;
            double t2 = Δa;
            if (!isAbridged) {
                // = Δf⋅(b/a)⋅ν⋅sin²φ
                t1 /= νden;
                // = Δa⋅(a/ν)
                t2 *= νden;
            }
            dstPts[dstOff++] = h + spcλ * cosφ + tZ * sinφ + t1 - t2;
        }
    }
    if (!derivate) {
        return null;
    }
    /*
         * At this point the (Abridged) Molodensky transformation is finished.
         * Code below this point is only for computing the derivative, if requested.
         * Note: variable names do not necessarily tell all the terms that they contain.
         */
    final Matrix matrix = Matrices.createDiagonal(getTargetDimensions(), getSourceDimensions());
    final double sinφcosφ = sinφ * cosφ;
    final double dν = eccentricitySquared * sinφcosφ / ν2den;
    final double dν3ρ = 3 * dν * (1 - eccentricitySquared) / ν2den;
    // double dXdλ     = spcλ;
    final double dYdλ = cmsλ * sinφ;
    final double dZdλ = cmsλ * cosφ;
    double dXdφ = dYdλ / cosφ;
    double dYdφ = -tZ * sinφ - cosφ * spcλ + t * (1 - 2 * sin2φ);
    double dZdφ = tZ * cosφ - sinφ * spcλ;
    if (isAbridged) {
        /*
             *   Δfmod  =  (a⋅Δf) + (f⋅Δa)
             *   t      =  2⋅Δfmod
             *   dXdh   =  0  so no need to set the matrix element.
             *   dYdh   =  0  so no need to set the matrix element.
             */
        dXdφ -= cmsλ * dν;
        dYdφ -= cmsφ * dν3ρ;
        dZdφ += t * cosφ * sinφ;
    } else {
        /*
             *   Δfmod  =  b⋅Δf
             *   t      =  Δf⋅[ν⋅(b/a) + ρ⋅(a/b)]    (real ν and ρ, without + h)
             *   ν         is actually ν + h
             *   ρ         is actually ρ + h
             */
        // Reminder: that ρ contains a h term.
        final double dρ = dν3ρ * νden * (semiMajor / ρ);
        // Reminder: that ν contains a h term.
        dXdφ -= dν * cmsλ * semiMajor / (νden * ν);
        dYdφ -= dρ * dZdφ - (Δfmod * (dν * 2 / (1 - eccentricitySquared) + (1 + 1 / ν2den) * (dν - dρ)) + Δa * (dν + 1) * eccentricitySquared) * sinφcosφ / νden;
        if (isSource3D) {
            final double dXdh = cmsλ / ν;
            final double dYdh = -cmsφ / ρ;
            matrix.setElement(0, 2, -dXdh * scaleX);
            matrix.setElement(1, 2, +dYdh * scaleY);
        }
        final double t1 = Δfmod * (dν * sin2φ + 2 * sinφcosφ);
        final double t2 = Δa * dν;
        dZdφ += t1 / νden + t2 * νden;
    }
    matrix.setElement(0, 0, 1 - spcλ * scaleX);
    matrix.setElement(1, 1, 1 + dYdφ * scaleY);
    matrix.setElement(0, 1, +dXdφ * scaleX);
    matrix.setElement(1, 0, -dYdλ * scaleY);
    if (isTarget3D) {
        matrix.setElement(2, 0, dZdλ);
        matrix.setElement(2, 1, dZdφ);
    }
    return matrix;
}