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

Example 11 with MatrixSIS

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

the class CoordinateSystems method swapAndScaleAxes.

/**
 * Returns an affine transform between two coordinate systems.
 * Only units and axes order (e.g. transforming from
 * ({@linkplain AxisDirection#NORTH North}, {@linkplain AxisDirection#WEST West}) to
 * ({@linkplain AxisDirection#EAST East}, {@linkplain AxisDirection#NORTH North})
 * are taken in account by this method.
 *
 * <div class="section">Conditions</div>
 * The two coordinate systems must implement the same GeoAPI coordinate system interface.
 * For example if {@code sourceCS} is a {@link org.opengis.referencing.cs.CartesianCS},
 * then {@code targetCS} must be a {@code CartesianCS} too.
 *
 * <div class="note"><b>Example:</b>
 * If coordinates in {@code sourceCS} are (<var>x</var>,<var>y</var>) tuples in metres
 * and coordinates in {@code targetCS} are (<var>-y</var>,<var>x</var>) tuples in centimetres,
 * then the transformation can be performed as below:
 *
 * {@preformat math
 *     ┌      ┐   ┌                ┐ ┌     ┐
 *     │-y(cm)│   │   0  -100    0 │ │ x(m)│
 *     │ x(cm)│ = │ 100     0    0 │ │ y(m)│
 *     │ 1    │   │   0     0    1 │ │ 1   │
 *     └      ┘   └                ┘ └     ┘
 * }
 * </div>
 *
 * @param  sourceCS  the source coordinate system.
 * @param  targetCS  the target coordinate system.
 * @return the conversion from {@code sourceCS} to {@code targetCS} as an affine transform.
 *         Only axis direction and units are taken in account.
 * @throws IllegalArgumentException if the CS are not of the same type, or axes do not match.
 * @throws IncommensurableException if the units are not compatible, or the conversion is non-linear.
 *
 * @see Matrices#createTransform(AxisDirection[], AxisDirection[])
 */
@SuppressWarnings("fallthrough")
public static Matrix swapAndScaleAxes(final CoordinateSystem sourceCS, final CoordinateSystem targetCS) throws IllegalArgumentException, IncommensurableException {
    ensureNonNull("sourceCS", sourceCS);
    ensureNonNull("targetCS", targetCS);
    if (!Classes.implementSameInterfaces(sourceCS.getClass(), targetCS.getClass(), CoordinateSystem.class)) {
        throw new IllegalArgumentException(Resources.format(Resources.Keys.IncompatibleCoordinateSystemTypes));
    }
    final AxisDirection[] srcAxes = getAxisDirections(sourceCS);
    final AxisDirection[] dstAxes = getAxisDirections(targetCS);
    final MatrixSIS matrix = Matrices.createTransform(srcAxes, dstAxes);
    assert Arrays.equals(srcAxes, dstAxes) == matrix.isIdentity() : matrix;
    /*
         * The previous code computed a matrix for swapping axes. Usually, this
         * matrix contains only 0 and 1 values with only one "1" value by row.
         * For example, the matrix operation for swapping x and y axes is:
         *          ┌ ┐   ┌         ┐ ┌ ┐
         *          │y│   │ 0  1  0 │ │x│
         *          │x│ = │ 1  0  0 │ │y│
         *          │1│   │ 0  0  1 │ │1│
         *          └ ┘   └         ┘ └ ┘
         * Now, take in account units conversions. Each matrix's element (j,i)
         * is multiplied by the conversion factor from sourceCS.getUnit(i) to
         * targetCS.getUnit(j). This is an element-by-element multiplication,
         * not a matrix multiplication. The last column is processed in a special
         * way, since it contains the offset values.
         */
    // == sourceCS.getDimension()
    final int sourceDim = matrix.getNumCol() - 1;
    // == targetCS.getDimension()
    final int targetDim = matrix.getNumRow() - 1;
    for (int j = 0; j < targetDim; j++) {
        final Unit<?> targetUnit = targetCS.getAxis(j).getUnit();
        for (int i = 0; i < sourceDim; i++) {
            if (matrix.getElement(j, i) == 0) {
                // (i.e. axes are orthogonal).
                continue;
            }
            final Unit<?> sourceUnit = sourceCS.getAxis(i).getUnit();
            if (Objects.equals(sourceUnit, targetUnit)) {
                // between source[i] and target[j].
                continue;
            }
            Number scale = 1;
            Number offset = 0;
            final Number[] coefficients = Units.coefficients(sourceUnit.getConverterToAny(targetUnit));
            switch(coefficients != null ? coefficients.length : -1) {
                // Fall through
                case 2:
                    scale = coefficients[1];
                // Fall through
                case 1:
                    offset = coefficients[0];
                case 0:
                    break;
                default:
                    throw new IncommensurableException(Resources.format(Resources.Keys.NonLinearUnitConversion_2, sourceUnit, targetUnit));
            }
            final DoubleDouble element = DoubleDouble.castOrCopy(matrix.getNumber(j, i));
            final DoubleDouble r = new DoubleDouble(element);
            r.multiply(scale);
            matrix.setNumber(j, i, r);
            r.setFrom(element);
            r.multiply(offset);
            r.add(matrix.getNumber(j, sourceDim));
            matrix.setNumber(j, sourceDim, r);
        }
    }
    return matrix;
}

Example 12 with MatrixSIS

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

the class LinearTransformBuilder method create.

/**
 * Creates a linear transform approximation from the source positions to the target positions.
 * This method assumes that source positions are precise and that all uncertainty is in the target positions.
 *
 * @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 fitted linear transform.
 * @throws FactoryException if the transform can not be created,
 *         for example because the source or target points have not be specified.
 *
 * @since 0.8
 */
@Override
@SuppressWarnings("serial")
public LinearTransform create(final MathTransformFactory factory) throws FactoryException {
    if (transform == null) {
        // Protect from changes.
        final double[][] sources = this.sources;
        final double[][] targets = this.targets;
        if (targets == null) {
            throw new InvalidGeodeticParameterException(noData());
        }
        final int sourceDim = (sources != null) ? sources.length : gridSize.length;
        final int targetDim = targets.length;
        correlation = new double[targetDim];
        final MatrixSIS matrix = Matrices.create(targetDim + 1, sourceDim + 1, ExtendedPrecisionMatrix.ZERO);
        matrix.setElement(targetDim, sourceDim, 1);
        for (int j = 0; j < targetDim; j++) {
            final double c;
            switch(sourceDim) {
                case 1:
                    {
                        final int row = j;
                        final Line line = new Line() {

                            @Override
                            public void setEquation(final Number slope, final Number y0) {
                                super.setEquation(slope, y0);
                                // Preserve the extended precision (double-double).
                                matrix.setNumber(row, 0, slope);
                                matrix.setNumber(row, 1, y0);
                            }
                        };
                        if (sources != null) {
                            c = line.fit(vector(sources[0]), vector(targets[j]));
                        } else {
                            c = line.fit(Vector.createSequence(0, 1, gridSize[0]), Vector.create(targets[j], false));
                        }
                        break;
                    }
                case 2:
                    {
                        final int row = j;
                        final Plane plan = new Plane() {

                            @Override
                            public void setEquation(final Number sx, final Number sy, final Number z0) {
                                super.setEquation(sx, sy, z0);
                                // Preserve the extended precision (double-double).
                                matrix.setNumber(row, 0, sx);
                                matrix.setNumber(row, 1, sy);
                                matrix.setNumber(row, 2, z0);
                            }
                        };
                        if (sources != null) {
                            c = plan.fit(vector(sources[0]), vector(sources[1]), vector(targets[j]));
                        } else
                            try {
                                c = plan.fit(gridSize[0], gridSize[1], Vector.create(targets[j], false));
                            } catch (IllegalArgumentException e) {
                                // This may happen if the z vector still contain some "NaN" values.
                                throw new InvalidGeodeticParameterException(noData(), e);
                            }
                        break;
                    }
                default:
                    {
                        throw new FactoryException(Errors.format(Errors.Keys.ExcessiveNumberOfDimensions_1, sourceDim));
                    }
            }
            correlation[j] = c;
        }
        transform = (LinearTransform) nonNull(factory).createAffineTransform(matrix);
    }
    return transform;
}

Example 13 with MatrixSIS

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

the class TimeDependentBWPTest method testSetPositionVectorTransformation.

/**
 * Tests the {@link TimeDependentBWP#setPositionVectorTransformation(Matrix, double)} method
 * using the example given in the EPSG database for operation method EPSG:1053.
 *
 * @throws NoninvertibleMatrixException Should not happen.
 */
@Test
@DependsOnMethod("testEpsgCalculation")
public void testSetPositionVectorTransformation() throws NoninvertibleMatrixException {
    /*
         * The transformation that we are going to test use as input
         * geocentric coordinates on ITRF2008 at epoch 2013.9.
         */
    final TimeDependentBWP p = create();
    final Date time = p.getTimeReference();
    time.setTime(time.getTime() + StrictMath.round((2013.9 - 1994) * JULIAN_YEAR_LENGTH));
    assertEquals(date("2013-11-25 11:24:00"), time);
    /*
         * Transform the point given in the EPSG example and compare with the coordinate
         * that we obtain if we do the calculation ourself using the intermediate values
         * given in EPSG example.
         */
    final MatrixSIS toGDA94 = MatrixSIS.castOrCopy(p.getPositionVectorTransformation(time));
    final MatrixSIS toITRF2008 = toGDA94.inverse();
    final MatrixSIS source = Matrices.create(4, 1, new double[] { -3789470.702, 4841770.411, -1690893.950, 1 });
    final MatrixSIS target = Matrices.create(4, 1, new double[] { -3789470.008, 4841770.685, -1690895.103, 1 });
    final MatrixSIS actual = toGDA94.multiply(source);
    compareWithExplicitCalculation(actual);
    /*
         * Now compare with the expected value given in the EPSG example. The 0.013 tolerance threshold
         * is for the X ordinate and has been determined empirically in testEpsgCalculation().
         */
    assertMatrixEquals("toGDA94", target, actual, 0.013);
    assertMatrixEquals("toITRF2008", source, toITRF2008.multiply(target), 0.013);
}

Example 14 with MatrixSIS

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

the class BursaWolfParametersTest method testProductOfInverse.

/**
 * Multiplies the <cite>ED87 to WGS 84</cite> parameters (EPSG:1146) transformation by its inverse and
 * verifies that the result is very close to the identity matrix, thanks to the double-double arithmetic.
 * This is an internal consistency test.
 *
 * @throws NoninvertibleMatrixException Should never happen.
 */
@Test
@DependsOnMethod("testGetPositionVectorTransformation")
public void testProductOfInverse() throws NoninvertibleMatrixException {
    final BursaWolfParameters bursaWolf = createED87_to_WGS84();
    final MatrixSIS toWGS84 = getPositionVectorTransformation(bursaWolf);
    final MatrixSIS toED87 = getPositionVectorTransformation(bursaWolf).inverse();
    final MatrixSIS product = toWGS84.multiply(toED87);
    assertTrue(product.isIdentity());
}

Example 15 with MatrixSIS

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

the class BursaWolfParametersTest method testGetPositionVectorTransformation.

/**
 * Tests {@link BursaWolfParameters#getPositionVectorTransformation(Date)}.
 * This test transform a point from WGS72 to WGS84, and conversely,
 * as documented in the example section of EPSG operation method 9606.
 *
 * @throws NoninvertibleMatrixException Should never happen.
 */
@Test
public void testGetPositionVectorTransformation() throws NoninvertibleMatrixException {
    final BursaWolfParameters bursaWolf = createWGS72_to_WGS84();
    final MatrixSIS toWGS84 = getPositionVectorTransformation(bursaWolf);
    final MatrixSIS toWGS72 = toWGS84.inverse();
    final MatrixSIS source = Matrices.create(4, 1, new double[] { 3657660.66, 255768.55, 5201382.11, 1 });
    final MatrixSIS target = Matrices.create(4, 1, new double[] { 3657660.78, 255778.43, 5201387.75, 1 });
    assertMatrixEquals("toWGS84", target, toWGS84.multiply(source), 0.01);
    assertMatrixEquals("toWGS72", source, toWGS72.multiply(target), 0.01);
    /*
         * Tests the optimized path for translation-only parameters.
         * Matrices having only translation terms are much easier to predict.
         */
    assertMatrixEquals("Translation only", new Matrix4(1, 0, 0, -168, 0, 1, 0, -60, 0, 0, 1, 320, 0, 0, 0, 1), getPositionVectorTransformation(createNTF_to_WGS84()), 0);
}