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

Example 16 with Matrix

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

the class BursaWolfParametersTest method testSetPositionVectorTransformation.

/**
 * Tests the {@link BursaWolfParameters#setPositionVectorTransformation(Matrix, double)} method.
 * This is an internal consistency test.
 */
@Test
@DependsOnMethod("testGetPositionVectorTransformation")
public void testSetPositionVectorTransformation() {
    final BursaWolfParameters bursaWolf = createED87_to_WGS84();
    final Matrix matrix = bursaWolf.getPositionVectorTransformation(null);
    final BursaWolfParameters actual = new BursaWolfParameters(bursaWolf.getTargetDatum(), bursaWolf.getDomainOfValidity());
    actual.setPositionVectorTransformation(matrix, 1E-10);
    assertEquals(bursaWolf, actual);
}

Example 17 with Matrix

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

the class BursaWolfParametersTest method testInvert.

/**
 * Tests {@link BursaWolfParameters#invert()}.
 *
 * @throws NoninvertibleMatrixException Should never happen.
 */
@Test
@DependsOnMethod("testProductOfInverse")
public void testInvert() throws NoninvertibleMatrixException {
    final BursaWolfParameters bursaWolf = createED87_to_WGS84();
    final Matrix original = getPositionVectorTransformation(bursaWolf).inverse();
    bursaWolf.invert();
    final Matrix inverse = getPositionVectorTransformation(bursaWolf);
    assertMatrixEquals("invert", original, inverse, 0.001);
}

Example 18 with Matrix

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

the class DefaultGeodeticDatumTest method testGetPositionVectorTransformation.

/**
 * Tests {@link DefaultGeodeticDatum#getPositionVectorTransformation(GeodeticDatum, Extent)}.
 */
@Test
@DependsOnMethod("testCreateAndSerialize")
public void testGetPositionVectorTransformation() {
    final Map<String, Object> properties = new HashMap<>();
    assertNull(properties.put(DefaultGeodeticDatum.NAME_KEY, "Invalid dummy datum"));
    /*
         * Associate two BursaWolfParameters, one valid only in a local area and the other one
         * valid globaly.  Note that we are building an invalid set of parameters, because the
         * source datum are not the same in both case. But for this test we are not interrested
         * in datum consistency - we only want any Bursa-Wolf parameters having different area
         * of validity.
         */
    // Local area (North Sea)
    final BursaWolfParameters local = BursaWolfParametersTest.createED87_to_WGS84();
    // Global area (World)
    final BursaWolfParameters global = BursaWolfParametersTest.createWGS72_to_WGS84();
    assertNull(properties.put(DefaultGeodeticDatum.BURSA_WOLF_KEY, new BursaWolfParameters[] { local, global }));
    /*
         * Build the datum using WGS 72 ellipsoid (so at least one of the BursaWolfParameters is real).
         */
    final DefaultGeodeticDatum datum = new DefaultGeodeticDatum(properties, GeodeticDatumMock.WGS72.getEllipsoid(), GeodeticDatumMock.WGS72.getPrimeMeridian());
    /*
         * Search for BursaWolfParameters around the North Sea area.
         */
    final DefaultGeographicBoundingBox areaOfInterest = new DefaultGeographicBoundingBox(-2, 8, 55, 60);
    final DefaultExtent extent = new DefaultExtent("Around the North Sea", areaOfInterest, null, null);
    Matrix matrix = datum.getPositionVectorTransformation(GeodeticDatumMock.NAD83, extent);
    assertNull("No BursaWolfParameters for NAD83", matrix);
    matrix = datum.getPositionVectorTransformation(GeodeticDatumMock.WGS84, extent);
    assertNotNull("BursaWolfParameters for WGS84", matrix);
    checkTransformationSignature(local, matrix, 0);
    /*
         * Expand the area of interest to something greater than North Sea, and test again.
         */
    areaOfInterest.setWestBoundLongitude(-8);
    matrix = datum.getPositionVectorTransformation(GeodeticDatumMock.WGS84, extent);
    assertNotNull("BursaWolfParameters for WGS84", matrix);
    checkTransformationSignature(global, matrix, 0);
    /*
         * Search in the reverse direction.
         */
    final DefaultGeodeticDatum targetDatum = new DefaultGeodeticDatum(GeodeticDatumMock.WGS84);
    matrix = targetDatum.getPositionVectorTransformation(datum, extent);
    // Expected result is the inverse.
    global.invert();
    checkTransformationSignature(global, matrix, 1E-6);
}

Example 19 with Matrix

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

the class MatricesTest method testCopy.

/**
 * Tests {@link Matrices#copy(Matrix)}
 */
@Test
public void testCopy() {
    final Matrix matrix = new Matrix3(10, 20, 30, 40, 50, 60, 70, 80, 90);
    final Matrix copy = Matrices.copy(matrix);
    assertNotSame("copy", matrix, copy);
    assertEquals("copy", matrix, copy);
}

Example 20 with Matrix

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

the class Envelopes method transform.

/**
 * Implementation of {@link #transform(MathTransform, Envelope)} with the opportunity to
 * save the projected center coordinate.
 *
 * @param  targetPt  after this method call, the center of the source envelope projected to the target CRS.
 *                   The length of this array must be the number of target dimensions.
 *                   May be {@code null} if this information is not needed.
 */
@SuppressWarnings("null")
private static GeneralEnvelope transform(final MathTransform transform, final Envelope envelope, final double[] targetPt) throws TransformException {
    if (transform.isIdentity()) {
        /*
             * Slight optimization: Just copy the envelope. Note that we need to set the CRS
             * to null because we don't know what the target CRS was supposed to be. Even if
             * an identity transform often imply that the target CRS is the same one than the
             * source CRS, it is not always the case. The metadata may be differents, or the
             * transform may be a datum shift without Bursa-Wolf parameters, etc.
             */
        final GeneralEnvelope transformed = new GeneralEnvelope(envelope);
        transformed.setCoordinateReferenceSystem(null);
        if (targetPt != null) {
            for (int i = envelope.getDimension(); --i >= 0; ) {
                targetPt[i] = transformed.getMedian(i);
            }
        }
        return transformed;
    }
    /*
         * Checks argument validity: envelope and math transform dimensions must be consistent.
         */
    final int sourceDim = transform.getSourceDimensions();
    final int targetDim = transform.getTargetDimensions();
    if (envelope.getDimension() != sourceDim) {
        throw new MismatchedDimensionException(Errors.format(Errors.Keys.MismatchedDimension_2, sourceDim, envelope.getDimension()));
    }
    /*
         * Allocates all needed objects. The value '3' below is because the following 'while'
         * loop uses a 'pointIndex' to be interpreted as a number in base 3 (see the comment
         * inside the loop).  The coordinate to transform must be initialized to the minimal
         * ordinate values. This coordinate will be updated in the 'switch' statement inside
         * the 'while' loop.
         */
    if (sourceDim >= 20) {
        // Maximal value supported by Formulas.pow3(int) is 19.
        throw new IllegalArgumentException(Errors.format(Errors.Keys.ExcessiveNumberOfDimensions_1));
    }
    int pointIndex = 0;
    boolean isDerivativeSupported = true;
    GeneralEnvelope transformed = null;
    final Matrix[] derivatives = new Matrix[Formulas.pow3(sourceDim)];
    final double[] ordinates = new double[derivatives.length * targetDim];
    final double[] sourcePt = new double[sourceDim];
    for (int i = sourceDim; --i >= 0; ) {
        sourcePt[i] = envelope.getMinimum(i);
    }
    // A window over a single coordinate in the 'ordinates' array.
    final DirectPositionView ordinatesView = new DirectPositionView.Double(ordinates, 0, targetDim);
    /*
         * Iterates over every minimal, maximal and median ordinate values (3 points) along each
         * dimension. The total number of iterations is 3 ^ (number of source dimensions).
         */
    transformPoint: while (true) {
        /*
             * Compute the derivative (optional operation). If this operation fails, we will
             * set a flag to 'false' so we don't try again for all remaining points. We try
             * to compute the derivative and the transformed point in a single operation if
             * we can. If we can not, we will compute those two information separately.
             *
             * Note that the very last point to be projected must be the envelope center.
             * There is usually no need to calculate the derivative for that last point,
             * but we let it does anyway for safety.
             */
        final int offset = pointIndex * targetDim;
        try {
            derivatives[pointIndex] = derivativeAndTransform(transform, sourcePt, ordinates, offset, isDerivativeSupported);
        } catch (TransformException e) {
            if (!isDerivativeSupported) {
                // Derivative were already disabled, so something went wrong.
                throw e;
            }
            isDerivativeSupported = false;
            transform.transform(sourcePt, 0, ordinates, offset, 1);
            // Log only if the above call was successful.
            recoverableException(Envelopes.class, e);
        }
        /*
             * The transformed point has been saved for future reuse after the enclosing
             * 'while' loop. Now add the transformed point to the destination envelope.
             */
        if (transformed == null) {
            transformed = new GeneralEnvelope(targetDim);
            for (int i = 0; i < targetDim; i++) {
                final double value = ordinates[offset + i];
                transformed.setRange(i, value, value);
            }
        } else {
            ordinatesView.offset = offset;
            transformed.add(ordinatesView);
        }
        /*
             * Get the next point coordinate. The 'coordinateIndex' variable is an index in base 3
             * having a number of digits equals to the number of source dimensions.  For example a
             * 4-D space have indexes ranging from "0000" to "2222" (numbers in base 3). The digits
             * are then mapped to minimal (0), maximal (1) or central (2) ordinates. The outer loop
             * stops when the counter roll back to "0000". Note that 'targetPt' must keep the value
             * of the last projected point, which must be the envelope center identified by "2222"
             * in the 4-D case.
             */
        int indexBase3 = ++pointIndex;
        for (int dim = sourceDim; --dim >= 0; indexBase3 /= 3) {
            switch(indexBase3 % 3) {
                // Continue the loop.
                case 0:
                    sourcePt[dim] = envelope.getMinimum(dim);
                    break;
                case 1:
                    sourcePt[dim] = envelope.getMaximum(dim);
                    continue transformPoint;
                case 2:
                    sourcePt[dim] = envelope.getMedian(dim);
                    continue transformPoint;
                // Should never happen
                default:
                    throw new AssertionError(indexBase3);
            }
        }
        break;
    }
    assert pointIndex == derivatives.length : pointIndex;
    /*
         * At this point we finished to build an envelope from all sampled positions. Now iterate
         * over all points. For each point, iterate over all line segments from that point to a
         * neighbor median point.  Use the derivate information for approximating the transform
         * behavior in that area by a cubic curve. We can then find analytically the curve extremum.
         *
         * The same technic is applied in transform(MathTransform, Rectangle2D), except that in
         * the Rectangle2D case the calculation was bundled right inside the main loop in order
         * to avoid the need for storage.
         */
    DirectPosition temporary = null;
    final DirectPositionView sourceView = new DirectPositionView.Double(sourcePt, 0, sourceDim);
    final CurveExtremum extremum = new CurveExtremum();
    for (pointIndex = 0; pointIndex < derivatives.length; pointIndex++) {
        final Matrix D1 = derivatives[pointIndex];
        if (D1 != null) {
            int indexBase3 = pointIndex, power3 = 1;
            for (int i = sourceDim; --i >= 0; indexBase3 /= 3, power3 *= 3) {
                final int digitBase3 = indexBase3 % 3;
                if (digitBase3 != 2) {
                    // Process only if we are not already located on the median along the dimension i.
                    final int medianIndex = pointIndex + power3 * (2 - digitBase3);
                    final Matrix D2 = derivatives[medianIndex];
                    if (D2 != null) {
                        final double xmin = envelope.getMinimum(i);
                        final double xmax = envelope.getMaximum(i);
                        final double x2 = envelope.getMedian(i);
                        final double x1 = (digitBase3 == 0) ? xmin : xmax;
                        final int offset1 = targetDim * pointIndex;
                        final int offset2 = targetDim * medianIndex;
                        for (int j = 0; j < targetDim; j++) {
                            extremum.resolve(x1, ordinates[offset1 + j], D1.getElement(j, i), x2, ordinates[offset2 + j], D2.getElement(j, i));
                            boolean isP2 = false;
                            do {
                                // Executed exactly twice, one for each extremum point.
                                final double x = isP2 ? extremum.ex2 : extremum.ex1;
                                if (x > xmin && x < xmax) {
                                    final double y = isP2 ? extremum.ey2 : extremum.ey1;
                                    if (y < transformed.getMinimum(j) || y > transformed.getMaximum(j)) {
                                        /*
                                             * At this point, we have determined that adding the extremum point
                                             * would expand the envelope. However we will not add that point
                                             * directly because its position may not be quite right (since we
                                             * used a cubic curve approximation). Instead, we project the point
                                             * on the envelope border which is located vis-à-vis the extremum.
                                             */
                                        for (int ib3 = pointIndex, dim = sourceDim; --dim >= 0; ib3 /= 3) {
                                            final double ordinate;
                                            if (dim == i) {
                                                // Position of the extremum.
                                                ordinate = x;
                                            } else
                                                switch(ib3 % 3) {
                                                    case 0:
                                                        ordinate = envelope.getMinimum(dim);
                                                        break;
                                                    case 1:
                                                        ordinate = envelope.getMaximum(dim);
                                                        break;
                                                    case 2:
                                                        ordinate = envelope.getMedian(dim);
                                                        break;
                                                    // Should never happen
                                                    default:
                                                        throw new AssertionError(ib3);
                                                }
                                            sourcePt[dim] = ordinate;
                                        }
                                        temporary = transform.transform(sourceView, temporary);
                                        transformed.add(temporary);
                                    }
                                }
                            } while ((isP2 = !isP2) == true);
                        }
                    }
                }
            }
            // Let GC do its job earlier.
            derivatives[pointIndex] = null;
        }
    }
    if (targetPt != null) {
        // Copy the coordinate of the center point.
        System.arraycopy(ordinates, ordinates.length - targetDim, targetPt, 0, targetDim);
    }
    return transformed;
}