Examples with DoubleDouble org.apache.sis.internal.util.DoubleDouble used on opensource projects

Example 6 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;
}

Example 7 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;
}

Example 8 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);
        }
    }
}

Example 9 with DoubleDouble

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

the class Matrices method createTransform.

/**
 * Implementation of {@code createTransform(…)} public methods expecting envelopes and/or axis directions.
 * Argument validity shall be verified by the caller.
 *
 * @param useEnvelopes {@code true} if source and destination envelopes shall be taken in account.
 *        If {@code false}, then source and destination envelopes will be ignored and can be null.
 */
@SuppressWarnings("null")
private static MatrixSIS createTransform(final Envelope srcEnvelope, final AxisDirection[] srcAxes, final Envelope dstEnvelope, final AxisDirection[] dstAxes, final boolean useEnvelopes) {
    final DirectPosition dstCorner, srcCorner, srcOppositeCorner;
    if (useEnvelopes) {
        dstCorner = dstEnvelope.getLowerCorner();
        srcCorner = srcEnvelope.getLowerCorner();
        srcOppositeCorner = srcEnvelope.getUpperCorner();
    } else {
        dstCorner = srcCorner = srcOppositeCorner = null;
    }
    /*
         * Unconditionally create extended precision matrix even if standard precision would be
         * enough because callers in other package may perform additional arithmetic operations
         * on it (for example org.apache.sis.referencing.cs.CoordinateSystems.swapAndScaleAxes).
         */
    final MatrixSIS matrix = new GeneralMatrix(dstAxes.length + 1, srcAxes.length + 1, false, 2);
    /*
         * Maps source axes to destination axes. If no axis is moved (for example if the user
         * want to transform (NORTH,EAST) to (SOUTH,EAST)), then source and destination index
         * will be equal.   If some axes are moved (for example if the user want to transform
         * (NORTH,EAST) to (EAST,NORTH)), then ordinates at index {@code srcIndex} will have
         * to be moved at index {@code dstIndex}.
         */
    for (int dstIndex = 0; dstIndex < dstAxes.length; dstIndex++) {
        boolean hasFound = false;
        final AxisDirection dstDir = dstAxes[dstIndex];
        final AxisDirection search = AxisDirections.absolute(dstDir);
        for (int srcIndex = 0; srcIndex < srcAxes.length; srcIndex++) {
            final AxisDirection srcDir = srcAxes[srcIndex];
            if (search.equals(AxisDirections.absolute(srcDir))) {
                if (hasFound) {
                    throw new IllegalArgumentException(Resources.format(Resources.Keys.ColinearAxisDirections_2, srcDir, dstDir));
                }
                hasFound = true;
                /*
                     * Set the matrix elements. Some matrix elements will never be set.
                     * They will be left to zero, which is their desired value.
                     */
                final boolean same = srcDir.equals(dstDir);
                if (useEnvelopes) {
                    /*
                         * See the comment in transform(Envelope, Envelope) for an explanation about why
                         * we use the lower/upper corners instead than getMinimum()/getMaximum() methods.
                         */
                    final DoubleDouble scale = new DoubleDouble(same ? +1 : -1, 0);
                    scale.multiply(dstEnvelope.getSpan(dstIndex));
                    scale.divide(srcEnvelope.getSpan(srcIndex));
                    final DoubleDouble translate = new DoubleDouble(scale);
                    translate.multiply((same ? srcCorner : srcOppositeCorner).getOrdinate(srcIndex));
                    translate.negate();
                    translate.add(dstCorner.getOrdinate(dstIndex));
                    matrix.setNumber(dstIndex, srcIndex, scale);
                    matrix.setNumber(dstIndex, srcAxes.length, translate);
                } else {
                    matrix.setElement(dstIndex, srcIndex, same ? +1 : -1);
                }
            }
        }
        if (!hasFound) {
            throw new IllegalArgumentException(Resources.format(Resources.Keys.CanNotMapAxisToDirection_1, dstAxes[dstIndex]));
        }
    }
    matrix.setElement(dstAxes.length, srcAxes.length, 1);
    return matrix;
}

Example 10 with DoubleDouble

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

the class MatrixSIS method multiply.

/**
 * Returns a new vector which is the result of multiplying this matrix with the specified vector.
 * In other words, returns {@code this} × {@code vector}. The length of the given vector must be
 * equal to the number of columns in this matrix, and the length of the returned vector will be
 * equal to the number of rows in this matrix.
 *
 * <div class="section">Relationship with coordinate operations</div>
 * In the context of coordinate operations, {@code Matrix.multiply(vector)} is related to
 * <code>{@linkplain AffineTransform#transform(double[], int, double[], int, int) AffineTransform.transform}(…)</code>
 * except that the last {@code vector} number is implicitly 1 in {@code AffineTransform} operations.
 * While this {@code multiply(double[])} method could be used for coordinate transformation, it is not its purpose.
 * This method is designed for occasional uses when accuracy is more important than performance.
 *
 * @param  vector  the vector to multiply to this matrix.
 * @return the result of {@code this} × {@code vector}.
 * @throws MismatchedMatrixSizeException if the length of the given vector is not equals to the
 *         number of columns in this matrix.
 *
 * @since 0.8
 */
public double[] multiply(final double[] vector) {
    final int numCol = getNumCol();
    if (vector.length != numCol) {
        throw new MismatchedMatrixSizeException(Errors.format(Errors.Keys.UnexpectedArrayLength_2, numCol, vector.length));
    }
    final double[] target = new double[getNumRow()];
    final DoubleDouble ele = new DoubleDouble();
    final DoubleDouble sum = new DoubleDouble();
    for (int j = 0; j < target.length; j++) {
        for (int i = 0; i < numCol; i++) {
            get(j, i, ele);
            ele.multiply(vector[i]);
            sum.add(ele);
        }
        target[j] = sum.value;
        sum.clear();
    }
    return target;
}