Search in sources :

Example 1 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class ScaleTransformTest method testExtendedPrecision.

/**
 * Verifies that {@link ScaleTransform} stores the error terms when they exist.
 */
@Test
@DependsOnMethod("testConstantDimension")
public void testExtendedPrecision() {
    final Number O = 0;
    final Number l = 1;
    final DoubleDouble r = DoubleDouble.createDegreesToRadians();
    final MatrixSIS matrix = Matrices.create(4, 4, new Number[] { r, O, O, O, O, r, O, O, O, O, l, O, O, O, O, l });
    final double[] elements = ((ExtendedPrecisionMatrix) matrix).getExtendedElements();
    assertTrue(r.value > r.error);
    // Paranoiac checks for making sure that next assertion will test something.
    assertFalse(r.error == 0);
    assertArrayEquals(new double[] { // Paranoiac check for making sure that getExtendedElements() is not broken.
    r.value, 0, 0, 0, 0, r.value, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, r.error, 0, 0, 0, 0, r.error, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, elements, 0);
    final ScaleTransform tr = new ScaleTransform(4, 4, elements);
    assertEquals("sourceDimensions", 3, tr.getSourceDimensions());
    assertEquals("targetDimensions", 3, tr.getTargetDimensions());
    Assert.assertMatrixEquals("matrix", matrix, tr.getMatrix(), 0.0);
    assertArrayEquals("elements", elements, tr.getExtendedElements(), 0.0);
    transform = tr;
    validate();
}
Also used : ExtendedPrecisionMatrix(org.apache.sis.internal.referencing.ExtendedPrecisionMatrix) MatrixSIS(org.apache.sis.referencing.operation.matrix.MatrixSIS) DoubleDouble(org.apache.sis.internal.util.DoubleDouble) Test(org.junit.Test) DependsOnMethod(org.apache.sis.test.DependsOnMethod)

Example 2 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class Equirectangular method createMathTransform.

/**
 * Creates an Equirectangular projection from the specified group of parameter values. This method is an
 * adaptation of {@link org.apache.sis.referencing.operation.projection.NormalizedProjection} constructor,
 * reproduced in this method because we will create an affine transform instead than the usual projection
 * classes.
 *
 * @param  factory     the factory to use if this constructor needs to create other math transforms.
 * @param  parameters  the parameter values that define the transform to create.
 * @return the map projection created from the given parameter values.
 * @throws FactoryException if an error occurred while creating the math transform.
 */
@Override
public MathTransform createMathTransform(final MathTransformFactory factory, final ParameterValueGroup parameters) throws FactoryException {
    final Parameters p = Parameters.castOrWrap(parameters);
    final ContextualParameters context = new ContextualParameters(this);
    double a = getAndStore(p, context, MapProjection.SEMI_MAJOR);
    double b = getAndStore(p, context, MapProjection.SEMI_MINOR);
    double λ0 = getAndStore(p, context, LONGITUDE_OF_ORIGIN);
    double φ0 = getAndStore(p, context, LATITUDE_OF_ORIGIN);
    double φ1 = getAndStore(p, context, STANDARD_PARALLEL);
    double fe = getAndStore(p, context, FALSE_EASTING);
    double fn = getAndStore(p, context, FALSE_NORTHING);
    /*
         * Perform following transformation, in that order. Note that following
         * AffineTransform convention, the Java code appears in reverse order:
         *
         *   1) Subtract φ0 to the latitude.
         *   2) Subtract λ0 to the longitude.
         *   3) Convert degrees to radians.
         *   4) Scale longitude by cos(φ1).
         */
    φ1 = toRadians(φ1);
    final MatrixSIS normalize = context.getMatrix(ContextualParameters.MatrixRole.NORMALIZATION);
    normalize.convertBefore(0, cos(φ1), null);
    context.normalizeGeographicInputs(λ0).convertBefore(1, null, -φ0);
    /*
         * At this point, we usually invoke 'denormalize.convertAfter(…, a, …)' where 'a' (the semi-major axis length)
         * is taken as the Earth radius (R). However quoting EPSG: "If the figure of the earth used is an ellipsoid
         * rather than a sphere then R should be calculated as the radius of the conformal sphere at the projection
         * origin at latitude φ1 using the formula for RC given in section 1.2, table 3".
         */
    if (a != b) {
        final double rs = b / a;
        final double sinφ1 = sin(φ1);
        a = b / (1 - (1 - rs * rs) * (sinφ1 * sinφ1));
    }
    final DoubleDouble k = new DoubleDouble(a);
    final MatrixSIS denormalize = context.getMatrix(ContextualParameters.MatrixRole.DENORMALIZATION);
    denormalize.convertAfter(0, k, new DoubleDouble(fe));
    denormalize.convertAfter(1, k, new DoubleDouble(fn));
    /*
         * Creates the ConcatenatedTransform, letting the factory returns the cached instance
         * if the caller already invoked this method previously (which usually do not happen).
         */
    MathTransform mt = context.completeTransform(factory, MathTransforms.identity(2));
    if (mt instanceof AffineTransform) {
        // Always true in Apache SIS implementation.
        mt = new ParameterizedAffine((AffineTransform) mt, context, true);
    }
    return mt;
}
Also used : Parameters(org.apache.sis.parameter.Parameters) ContextualParameters(org.apache.sis.referencing.operation.transform.ContextualParameters) ContextualParameters(org.apache.sis.referencing.operation.transform.ContextualParameters) MathTransform(org.opengis.referencing.operation.MathTransform) AffineTransform(java.awt.geom.AffineTransform) ParameterizedAffine(org.apache.sis.internal.referencing.j2d.ParameterizedAffine) MatrixSIS(org.apache.sis.referencing.operation.matrix.MatrixSIS) DoubleDouble(org.apache.sis.internal.util.DoubleDouble)

Example 3 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class ContextualParameters method normalizeGeographicInputs.

/**
 * Prepends a normalization step converting input ordinates in the two first dimensions from degrees to radians.
 * The normalization can optionally subtract the given λ₀ value (in degrees) from the longitude.
 *
 * <p>Invoking this method is equivalent to {@linkplain java.awt.geom.AffineTransform#concatenate concatenating}
 * the normalization matrix with the following matrix. This will have the effect of applying the conversion
 * described above before any other operation:</p>
 *
 * <center>{@include formulas.html#NormalizeGeographic}</center>
 *
 * @param  λ0  longitude of the central meridian, in degrees.
 * @return the normalization affine transform as a matrix.
 *         Callers can change that matrix directly if they want to apply additional normalization operations.
 * @throws IllegalStateException if this {@code ContextualParameter} has been made unmodifiable.
 */
public synchronized MatrixSIS normalizeGeographicInputs(final double λ0) {
    ensureModifiable();
    /*
         * In theory the check for (λ0 != 0) is useless. However Java has a notion of negative zero, and we want
         * to avoid negative zeros because we do not want them to appear in WKT formatting of matrix elements.
         */
    final DoubleDouble toRadians = DoubleDouble.createDegreesToRadians();
    DoubleDouble offset = null;
    if (λ0 != 0) {
        offset = new DoubleDouble(-λ0);
        offset.multiply(toRadians);
    }
    // Must be the same instance, not a copy.
    final MatrixSIS normalize = (MatrixSIS) this.normalize;
    normalize.convertBefore(0, toRadians, offset);
    normalize.convertBefore(1, toRadians, null);
    return normalize;
}
Also used : MatrixSIS(org.apache.sis.referencing.operation.matrix.MatrixSIS) DoubleDouble(org.apache.sis.internal.util.DoubleDouble)

Example 4 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class GeneralMatrix method setElements.

/**
 * Sets all matrix elements like {@link #setElements(double[])}, but from an array of {@code Number} instead
 * of {@code double}. The main purpose of this method is to fetch the {@link DoubleDouble#error} terms when
 * such instances are found.
 *
 * <p><b>Restrictions:</b></p>
 * <ul>
 *   <li>This matrix must use extended-precision elements as by {@link #createExtendedPrecision(int, int, boolean)}.</li>
 *   <li>If this method returns {@code false}, then error terms are <strong>not</strong> initialized - they
 *       may have any values.</li>
 * </ul>
 *
 * @param  newValues  the new matrix elements in a row-major array.
 * @return {@code true} if at leat one {@link DoubleDouble} instance has been found, in which case all
 *         errors terms have been initialized, or {@code false} otherwise, in which case no error term
 *         has been initialized (this is a <cite>all or nothing</cite> operation).
 * @throws IllegalArgumentException if the given array does not have the expected length.
 *
 * @see Matrices#create(int, int, Number[])
 */
final boolean setElements(final Number[] newValues) {
    // Protection against accidental changes.
    final int numRow = this.numRow;
    final int numCol = this.numCol;
    final int length = numRow * numCol;
    if (newValues.length != length) {
        throw new IllegalArgumentException(Errors.format(Errors.Keys.UnexpectedArrayLength_2, length, newValues.length));
    }
    boolean isExtended = false;
    for (int i = 0; i < length; i++) {
        Number value = newValues[i];
        if (DoubleDouble.shouldConvert(value)) {
            value = new DoubleDouble(value);
        }
        final double element = value.doubleValue();
        elements[i] = element;
        final double error;
        if (value instanceof DoubleDouble) {
            error = ((DoubleDouble) value).error;
            /*
                 * If this is the first time that we found an explicit error term, then we need to
                 * initialize all elements before the current one because they were left unitialized
                 * (i.e. we perform lazy initialization).
                 */
            if (!isExtended) {
                isExtended = true;
                for (int j = 0; j < i; j++) {
                    elements[j + length] = DoubleDouble.errorForWellKnownValue(elements[j]);
                }
            }
        } else {
            /*
                 * For any kind of numbers other than DoubleDoube, calculate the error term only if we know
                 * that the final matrix will use extended precision (i.e. we previously found at least one
                 * DoubleDouble instance). Otherwise skip the error calculation since maybe it will be discarded.
                 */
            if (!isExtended) {
                continue;
            }
            error = DoubleDouble.errorForWellKnownValue(element);
        }
        elements[i + length] = error;
    }
    return isExtended;
}
Also used : DoubleDouble(org.apache.sis.internal.util.DoubleDouble)

Example 5 with DoubleDouble

use of org.apache.sis.internal.util.DoubleDouble in project sis by apache.

the class GeneralMatrix method setToProduct.

/**
 * Sets this matrix to the product of the given matrices: {@code this = A × B}.
 * The matrix sizes much match - this is not verified unless assertions are enabled.
 */
final void setToProduct(final Matrix A, final Matrix B) {
    // Protection against accidental changes.
    final int numRow = this.numRow;
    final int numCol = this.numCol;
    final int nc = A.getNumCol();
    assert B.getNumRow() == nc;
    assert numRow == A.getNumRow() && numCol == B.getNumCol();
    /*
         * Get the matrix element values, together with the error terms if the matrix
         * use extended precision (double-double arithmetic).
         */
    final double[] eltA = getExtendedElements(A, numRow, nc, false);
    final double[] eltB = getExtendedElements(B, nc, numCol, false);
    // Where error terms start.
    final int errorOffset = numRow * numCol;
    final int errA = numRow * nc;
    final int errB = nc * numCol;
    /*
         * Compute the product, to be stored directly in 'this'.
         */
    final DoubleDouble dot = new DoubleDouble();
    final DoubleDouble sum = new DoubleDouble();
    for (int k = 0, j = 0; j < numRow; j++) {
        for (int i = 0; i < numCol; i++) {
            sum.clear();
            double max = 0;
            // Index of values in a single column of B.
            int iB = i;
            // Index of values in a single row of A.
            int iA = j * nc;
            final int nextRow = iA + nc;
            while (iA < nextRow) {
                dot.setFrom(eltA, iA, errA);
                dot.multiply(eltB, iB, errB);
                sum.add(dot);
                // Move to next row of B.
                iB += numCol;
                // Move to next column of A.
                iA++;
                final double value = Math.abs(dot.value);
                if (value > max)
                    max = value;
            }
            if (Math.abs(sum.value) < Math.ulp(max) * ZERO_THRESHOLD) {
                // Sum is not significant according double arithmetic.
                sum.clear();
            }
            sum.storeTo(elements, k++, errorOffset);
        }
    }
}
Also used : DoubleDouble(org.apache.sis.internal.util.DoubleDouble)

Aggregations

DoubleDouble (org.apache.sis.internal.util.DoubleDouble)39 MatrixSIS (org.apache.sis.referencing.operation.matrix.MatrixSIS)5 Test (org.junit.Test)3 ExtendedPrecisionMatrix (org.apache.sis.internal.referencing.ExtendedPrecisionMatrix)2 DirectPosition (org.opengis.geometry.DirectPosition)2 AxisDirection (org.opengis.referencing.cs.AxisDirection)2 AffineTransform (java.awt.geom.AffineTransform)1 BigDecimal (java.math.BigDecimal)1 IncommensurableException (javax.measure.IncommensurableException)1 ParameterizedAffine (org.apache.sis.internal.referencing.j2d.ParameterizedAffine)1 Parameters (org.apache.sis.parameter.Parameters)1 Matrix4 (org.apache.sis.referencing.operation.matrix.Matrix4)1 ContextualParameters (org.apache.sis.referencing.operation.transform.ContextualParameters)1 DependsOnMethod (org.apache.sis.test.DependsOnMethod)1 MismatchedDimensionException (org.opengis.geometry.MismatchedDimensionException)1 CoordinateSystem (org.opengis.referencing.cs.CoordinateSystem)1 MathTransform (org.opengis.referencing.operation.MathTransform)1