Java Code Examples for org.opengis.referencing.operation.Matrix#getNumRow()

/**
 * If a sequence of 3 transforms are (inverse projection) → (affine) → (projection) where
 * the (projection) and (inverse projection) steps are the inverse of each other, returns
 * the matrix of the affine transform step. Otherwise returns {@code null}. This method
 * accepts also (projection) → (affine) → (inverse projection) sequence, but such sequences
 * should be much more unusual.
 *
 * @param  projection       either {@link NormalizedProjection} or {@link Inverse}.
 * @param  other            the arbitrary transforms to be concatenated with the given projection.
 * @param  applyOtherFirst  whether {@code other} is concatenated before {@code projection} or the converse.
 * @return the 3×3 matrix of the affine transform step, or {@code null} if none.
 */
private static Matrix getMiddleMatrix(final AbstractMathTransform projection, final MathTransform other,
        final boolean applyOtherFirst)
{
    final List<MathTransform> steps = MathTransforms.getSteps(other);
    if (steps.size() == 2) try {
        final int oi = applyOtherFirst ? 0 : 1;
        if (projection.equals(steps.get(oi).inverse(), ComparisonMode.IGNORE_METADATA)) {
            final Matrix m = MathTransforms.getMatrix(steps.get(oi ^ 1));
            if (Matrices.isAffine(m) && m.getNumRow() == DIMENSION+1 && m.getNumCol() == DIMENSION+1) {
                return m;
            }
        }
    } catch (NoninvertibleTransformException e) {
        Logging.recoverableException(Logging.getLogger(Loggers.COORDINATE_OPERATION),
                (projection instanceof NormalizedProjection) ? NormalizedProjection.class : projection.getClass(),
                "tryConcatenate", e);
    }
    return null;
}
 

/**
 * Returns {@code true} if the specified transform is likely to exists only for axis swapping
 * and/or unit conversions. The heuristic rule checks if the transform is backed by a square
 * matrix with exactly one non-null value in each row and each column.
 */
private static boolean isIgnorable(final MathTransform transform) {
    final Matrix matrix = MathTransforms.getMatrix(transform);
    if (matrix != null) {
        final int size = matrix.getNumRow();
        if (matrix.getNumCol() == size) {
            for (int j=0; j<size; j++) {
                int n1=0, n2=0;
                for (int i=0; i<size; i++) {
                    if (matrix.getElement(j,i) != 0) n1++;
                    if (matrix.getElement(i,j) != 0) n2++;
                }
                if (n1 != 1 || n2 != 1) {
                    return false;
                }
            }
            return true;
        }
    }
    return false;
}
 

/**
 * Asserts that the given matrix is equals to the expected one, up to the given tolerance value.
 *
 * @param message   Header of the exception message in case of failure, or {@code null} if none.
 * @param expected  The expected matrix, which may be {@code null}.
 * @param actual    The matrix to compare, or {@code null}.
 * @param tolerance The tolerance threshold.
 */
public static void assertMatrixEquals(final String message, final Matrix expected, final Matrix actual, final double tolerance) {
    if (isNull(message, expected, actual)) {
        return;
    }
    final int numRow = actual.getNumRow();
    final int numCol = actual.getNumCol();
    assertEquals("numRow", expected.getNumRow(), numRow);
    assertEquals("numCol", expected.getNumCol(), numCol);
    for (int j=0; j<numRow; j++) {
        for (int i=0; i<numCol; i++) {
            final double e = expected.getElement(j,i);
            final double a = actual.getElement(j,i);
            if (!(StrictMath.abs(e - a) <= tolerance) && Double.doubleToLongBits(a) != Double.doubleToLongBits(e)) {
                fail("Matrix.getElement(" + j + ", " + i + "): expected " + e + " but got " + a);
            }
        }
    }
}
 

/**
 * Asserts that the given matrix is diagonal, and that all elements on the diagonal are equal
 * to the given values. The matrix doesn't need to be square. The last row is handled especially
 * if the {@code affine} argument is {@code true}.
 *
 * @param expected   the values which are expected on the diagonal. If the length of this array
 *                   is less than the matrix size, then the last element in the array is repeated
 *                   for all remaining diagonal elements.
 * @param affine     if {@code true}, then the last row is expected to contains the value 1
 *                   in the last column, and all other columns set to 0.
 * @param matrix     the matrix to test.
 * @param tolerance  the tolerance threshold while comparing floating point values.
 */
public static void assertDiagonalEquals(final double[] expected, final boolean affine,
        final Matrix matrix, final double tolerance)
{
    final int numRows = matrix.getNumRow();
    final int numCols = matrix.getNumCol();
    final StringBuilder buffer = new StringBuilder("matrix(");
    final int bufferBase = buffer.length();
    for (int j=0; j<numRows; j++) {
        for (int i=0; i<numCols; i++) {
            final double e;
            if (affine && j == numRows-1) {
                e = (i == numCols-1) ? 1 : 0;
            } else if (i == j) {
                e = expected[min(expected.length-1, i)];
            } else {
                e = 0;
            }
            buffer.setLength(bufferBase);
            assertEquals(buffer.append(j).append(',').append(i).append(')').toString(),
                    e, matrix.getElement(j, i), tolerance);
        }
    }
}
 

/**
 * Creates a new matrix which is a copy of the given matrix.
 *
 * <div class="note"><b>Implementation note:</b>
 * For square matrix with a size between {@value org.apache.sis.referencing.operation.matrix.Matrix1#SIZE}
 * and {@value org.apache.sis.referencing.operation.matrix.Matrix4#SIZE} inclusive, the returned matrix is
 * usually an instance of one of {@link Matrix1} … {@link Matrix4} subtypes.</div>
 *
 * @param  matrix  the matrix to copy, or {@code null}.
 * @return a copy of the given matrix, or {@code null} if the given matrix was null.
 *
 * @see MatrixSIS#clone()
 * @see MatrixSIS#castOrCopy(Matrix)
 */
public static MatrixSIS copy(final Matrix matrix) {
    if (matrix == null) {
        return null;
    }
    final int size = matrix.getNumRow();
    if (size != matrix.getNumCol()) {
        return new NonSquareMatrix(matrix);
    }
    if (!isExtendedPrecision(matrix)) {
        switch (size) {
            case 1: return new Matrix1(matrix);
            case 2: return new Matrix2(matrix);
            case 3: return new Matrix3(matrix);
            case 4: return new Matrix4(matrix);
        }
    }
    return new GeneralMatrix(matrix);
}
 

/**
 * Solves {@code X} × <var>U</var> = {@code Y}.
 * This method is an adaptation of the {@code LUDecomposition} class of the JAMA matrix package.
 *
 * @param  X  the matrix to invert.
 * @param  Y  the desired result of {@code X} × <var>U</var>.
 * @throws NoninvertibleMatrixException if the {@code X} matrix is not square or singular.
 */
static MatrixSIS solve(final Matrix X, final Matrix Y) throws NoninvertibleMatrixException {
    final int size   = X.getNumRow();
    final int numCol = X.getNumCol();
    if (numCol != size) {
        throw new NoninvertibleMatrixException(Resources.format(Resources.Keys.NonInvertibleMatrix_2, size, numCol));
    }
    final int innerSize = Y.getNumCol();
    GeneralMatrix.ensureNumRowMatch(size, Y.getNumRow(), innerSize);
    double[] eltY = null;
    if (Y instanceof GeneralMatrix) {
        eltY = ((GeneralMatrix) Y).elements;
        if (eltY.length == size * innerSize) {
            eltY = null;                            // Matrix does not contains error terms.
        }
    }
    return solve(X, Y, eltY, size, innerSize, true);
}
 

/**
 * Creates a new test suite.
 */
public ProjectiveTransformTest() {
    this(new MathTransformFactoryBase() {
        @Override
        public MathTransform createAffineTransform(final Matrix matrix) {
            final ProjectiveTransform pt;
            if (matrix.getNumRow() == 3 && matrix.getNumCol() == 3) {
                pt = new ProjectiveTransform2D(matrix);
            } else {
                pt = new ProjectiveTransform(matrix);
            }
            MathTransform tr = pt.optimize();
            if (tr != pt) {
                /*
                 * Opportunistically tests ScaledTransform together with ProjectiveTransform.
                 * We takes ScaledTransform as a reference implementation since it is simpler.
                 */
                tr = new TransformResultComparator(tr, pt, 0);
            }
            return tr;
        }
    });
}
 

/**
 * Computes the conversion factors needed for calls to {@link DatumShiftGrid#interpolateInCell(double, double, double[])}.
 * This method takes only the {@value DatumShiftGrid#INTERPOLATED_DIMENSIONS} first dimensions. If a conversion factor can
 * not be computed, then it is set to NaN.
 */
@SuppressWarnings("fallthrough")
private void computeConversionFactors() {
    scaleX = Double.NaN;
    scaleY = Double.NaN;
    x0     = Double.NaN;
    y0     = Double.NaN;
    if (grid != null) {
        final LinearTransform coordinateToGrid = grid.getCoordinateToGrid();
        final double toStandardUnit = Units.toStandardUnit(grid.getCoordinateUnit());
        if (!Double.isNaN(toStandardUnit)) {
            final Matrix m = coordinateToGrid.getMatrix();
            if (Matrices.isAffine(m)) {
                final int n = m.getNumCol() - 1;
                switch (m.getNumRow()) {
                    default: y0 = m.getElement(1,n); scaleY = diagonal(m, 1, n) / toStandardUnit;   // Fall through
                    case 1:  x0 = m.getElement(0,n); scaleX = diagonal(m, 0, n) / toStandardUnit;
                    case 0:  break;
                }
            }
        }
    }
}
 

/**
 * Ensures that the given matrix has the given dimension.
 * This is a convenience method for subclasses.
 */
static void ensureSizeMatch(final int numRow, final int numCol, final Matrix matrix)
        throws MismatchedMatrixSizeException
{
    final int othRow = matrix.getNumRow();
    final int othCol = matrix.getNumCol();
    if (numRow != othRow || numCol != othCol) {
        throw new MismatchedMatrixSizeException(Errors.format(
                Errors.Keys.MismatchedMatrixSize_4, numRow, numCol, othRow, othCol));
    }
}
 

/**
 * Creates an arbitrary linear transform from the specified matrix. Usually the matrix
 * {@linkplain org.apache.sis.referencing.operation.matrix.MatrixSIS#isAffine() is affine},
 * but this is not mandatory. Non-affine matrix will define a projective transform.
 *
 * <p>If the transform input dimension is {@code M}, and output dimension is {@code N},
 * then the given matrix shall have size {@code [N+1][M+1]}.
 * The +1 in the matrix dimensions allows the matrix to do a shift, as well as a rotation.
 * The {@code [M][j]} element of the matrix will be the <var>j</var>'th coordinate of the moved origin.</p>
 *
 * @param  matrix  the matrix used to define the linear transform.
 * @return the linear (usually affine) transform.
 *
 * @see #getMatrix(MathTransform)
 * @see DefaultMathTransformFactory#createAffineTransform(Matrix)
 */
public static LinearTransform linear(final Matrix matrix) {
    ArgumentChecks.ensureNonNull("matrix", matrix);
    final int sourceDimension = matrix.getNumCol() - 1;
    final int targetDimension = matrix.getNumRow() - 1;
    if (sourceDimension == targetDimension) {
        if (matrix.isIdentity()) {
            return identity(sourceDimension);
        }
        if (Matrices.isAffine(matrix)) {
            switch (sourceDimension) {
                case 1: {
                    return linear(matrix.getElement(0,0), matrix.getElement(0,1));
                }
                case 2: {
                    if (matrix instanceof ExtendedPrecisionMatrix) {
                        return new AffineTransform2D(((ExtendedPrecisionMatrix) matrix).getExtendedElements());
                    } else {
                        return new AffineTransform2D(
                                matrix.getElement(0,0), matrix.getElement(1,0),
                                matrix.getElement(0,1), matrix.getElement(1,1),
                                matrix.getElement(0,2), matrix.getElement(1,2));
                    }
                }
            }
        } else if (sourceDimension == 2) {
            return new ProjectiveTransform2D(matrix);
        }
    }
    final LinearTransform candidate = CopyTransform.create(matrix);
    if (candidate != null) {
        return candidate;
    }
    return new ProjectiveTransform(matrix).optimize();
}
 

/**
 * Sets the {@link #scale} and {@link #offset} terms from the given function.
 *
 * @param  function  the transform to set.
 * @throws IllegalArgumentException if this method does not recognize the given transform.
 */
private void setLinearTerms(final LinearTransform function) throws IllegalArgumentException {
    final Matrix m = function.getMatrix();
    final int numRow = m.getNumRow();
    final int numCol = m.getNumCol();
    if (numRow != 2 || numCol != 2) {
        final Integer two = 2;
        throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedMatrixSize_4, two, two, numRow, numCol));
    }
    scale  = m.getElement(0, 0);
    offset = m.getElement(0, 1);
}
 

/**
 * Returns the inverse of the given matrix.
 *
 * @param  matrix  the matrix to inverse, or {@code null}.
 * @return the inverse of this matrix, or {@code null} if the given matrix was null.
 * @throws NoninvertibleMatrixException if the given matrix is not invertible.
 *
 * @see MatrixSIS#inverse()
 *
 * @since 0.6
 */
public static MatrixSIS inverse(final Matrix matrix) throws NoninvertibleMatrixException {
    if (matrix == null) {
        return null;
    } else if (matrix instanceof MatrixSIS) {
        return ((MatrixSIS) matrix).inverse();  // Maybe the subclass override that method.
    } else if (matrix.getNumRow() != matrix.getNumCol()) {
        return new NonSquareMatrix(matrix).inverse();
    } else {
        return Solver.inverse(matrix, true);
    }
}
 

/**
 * If a transform can be created from the given matrix, returns it.
 * Otherwise returns {@code null}.
 */
static CopyTransform create(final Matrix matrix) {
    final int srcDim = matrix.getNumCol() - 1;
    final int dstDim = matrix.getNumRow() - 1;
    for (int i=0; i <= srcDim; i++) {
        if (matrix.getElement(dstDim, i) != (i == srcDim ? 1 : 0)) {
            // Not an affine transform (ignoring if square or not).
            return null;
        }
    }
    final int[] indices = new int[dstDim];
    for (int j=0; j<dstDim; j++) {
        if (matrix.getElement(j, srcDim) != 0) {
            // The matrix has translation terms.
            return null;
        }
        boolean found = false;
        for (int i=0; i<srcDim; i++) {
            final double elt = matrix.getElement(j, i);
            if (elt != 0) {
                if (elt != 1 || found) {
                    // Not a simple copy operation.
                    return null;
                }
                indices[j] = i;
                found = true;
            }
        }
        if (!found) {
            // Target coordinate unconditionally set to 0 (not a copy).
            return null;
        }
    }
    return new CopyTransform(srcDim, indices);
}
 

/**
 * Compares the given matrices for equality, using the given relative or absolute tolerance threshold.
 * The matrix elements are compared as below:
 *
 * <ul>
 *   <li>{@link Double#NaN} values are considered equals to all other NaN values</li>
 *   <li>Infinite values are considered equal to other infinite values of the same sign</li>
 *   <li>All other values are considered equal if the absolute value of their difference is
 *       smaller than or equals to the threshold described below.</li>
 * </ul>
 *
 * If {@code relative} is {@code true}, then for any pair of values <var>v1</var><sub>j,i</sub>
 * and <var>v2</var><sub>j,i</sub> to compare, the tolerance threshold is scaled by
 * {@code max(abs(v1), abs(v2))}. Otherwise the threshold is used as-is.
 *
 * @param  m1        the first matrix to compare, or {@code null}.
 * @param  m2        the second matrix to compare, or {@code null}.
 * @param  epsilon   the tolerance value.
 * @param  relative  if {@code true}, then the tolerance value is relative to the magnitude
 *                   of the matrix elements being compared.
 * @return {@code true} if the values of the two matrix do not differ by a quantity greater
 *         than the given tolerance threshold.
 *
 * @see MatrixSIS#equals(Matrix, double)
 */
public static boolean equals(final Matrix m1, final Matrix m2, final double epsilon, final boolean relative) {
    if (m1 != m2) {
        if (m1 == null || m2 == null) {
            return false;
        }
        final int numRow = m1.getNumRow();
        if (numRow != m2.getNumRow()) {
            return false;
        }
        final int numCol = m1.getNumCol();
        if (numCol != m2.getNumCol()) {
            return false;
        }
        for (int j=0; j<numRow; j++) {
            for (int i=0; i<numCol; i++) {
                final double v1 = m1.getElement(j, i);
                final double v2 = m2.getElement(j, i);
                double tolerance = epsilon;
                if (relative) {
                    tolerance *= Math.max(Math.abs(v1), Math.abs(v2));
                }
                if (!(Math.abs(v1 - v2) <= tolerance)) {
                    if (Numerics.equals(v1, v2)) {
                        // Special case for NaN and infinite values.
                        continue;
                    }
                    return false;
                }
            }
        }
    }
    return true;
}
 

/**
 * Sets all parameter values to the element value in the specified matrix.
 * After this method call, {@link #values} will returns only the elements
 * different from the default value.
 *
 * @param  matrix  the matrix to copy in this group of parameters.
 */
final void setMatrix(final Matrix matrix) {
    final int numRow = matrix.getNumRow();
    final int numCol = matrix.getNumCol();
    dimensions[0].setValue(numRow);
    dimensions[1].setValue(numCol);
    values = null;
    final int[] indices = new int[2];
    for (int j=0; j<numRow; j++) {
        indices[0] = j;
        ParameterValue<?>[] row = null;
        for (int i=0; i<numCol; i++) {
            indices[1] = i;
            ParameterDescriptor<E> descriptor = descriptors.getElementDescriptor(indices);
            final E def = descriptor.getDefaultValue();
            final double element = matrix.getElement(j,i);
            if (!(def instanceof Number) || !Numerics.equalsIgnoreZeroSign(element, ((Number) def).doubleValue())) {
                final ParameterValue<?> value = descriptor.createValue();
                value.setValue(element);
                if (row == null) {
                    row = new ParameterValue<?>[numCol];
                    if (values == null) {
                        values = new ParameterValue<?>[numRow][];
                    }
                    values[j] = row;
                }
                row[i] = value;
            }
        }
    }
}
 

/**
 * Estimates the derivative at the given grid indices. Derivatives must be consistent with values given by
 * {@link #interpolateInCell(double, double, double[])} at adjacent positions. For a two-dimensional grid,
 * {@code tₐ(x,y)} an abbreviation for {@code interpolateInCell(gridX, gridY, …)[a]} and for <var>x</var>
 * and <var>y</var> integers, the derivative is:
 *
 * {@preformat math
 *   ┌                   ┐         ┌                                                        ┐
 *   │  ∂t₀/∂x   ∂t₀/∂y  │    =    │  t₀(x+1,y) - t₀(x,y) + 1      t₀(x,y+1) - t₀(x,y)      │
 *   │  ∂t₁/∂x   ∂t₁/∂y  │         │  t₁(x+1,y) - t₁(x,y)          t₁(x,y+1) - t₁(x,y) + 1  │
 *   └                   ┘         └                                                        ┘
 * }
 *
 * <h4>Extrapolations</h4>
 * Derivatives must be consistent with {@link #interpolateInCell(double, double, double[])} even when the
 * given coordinates are outside the grid. The {@code interpolateInCell(…)} contract in such cases is to
 * compute the translation vector at the closest position in the grid. A consequence of this contract is
 * that translation vectors stay constant when moving along at least one direction outside the grid.
 * Consequences on the derivative matrix are as below:
 *
 * <ul>
 *   <li>If both {@code gridX} and {@code gridY} are outside the grid, then the derivative is the identity matrix.</li>
 *   <li>If only {@code gridX} is outside the grid, then only the first column is set to [1, 0, …].
 *       The second column is set to the derivative of the closest cell at {@code gridY} position.</li>
 *   <li>If only {@code gridY} is outside the grid, then only the second column is set to [0, 1, …].
 *       The first column is set to the derivative of the closest cell at {@code gridX} position.</li>
 * </ul>
 *
 * @param  gridX  first grid coordinate of the point for which to get the translation.
 * @param  gridY  second grid coordinate of the point for which to get the translation.
 * @return the derivative at the given location.
 *
 * @see #isCellInGrid(double, double)
 * @see #interpolateInCell(double, double, double[])
 */
public Matrix derivativeInCell(double gridX, double gridY) {
    final int xmax = gridSize[0] - 2;
    final int ymax = gridSize[1] - 2;
    int ix = (int) gridX;
    int iy = (int) gridY;
    if (ix < 0 || ix > xmax || iy < 0 || iy > ymax) {
        final double[] gridCoordinates = {gridX, gridY};
        replaceOutsideGridCoordinates(gridCoordinates);
        gridX = gridCoordinates[0];
        gridY = gridCoordinates[1];
        ix = Math.max(0, Math.min(xmax, (int) gridX));
        iy = Math.max(0, Math.min(ymax, (int) gridY));
    }
    gridX -= ix;
    gridY -= iy;
    final boolean skipX = (gridX < 0 || gridX > 1);
    final boolean skipY = (gridY < 0 || gridY > 1);
    final Matrix derivative = Matrices.createDiagonal(getTranslationDimensions(), gridSize.length);
    for (int j=derivative.getNumRow(); --j>=0;) {
        final double r00 = getCellValue(j, ix,   iy  );
        final double r01 = getCellValue(j, ix+1, iy  );       // Naming convention: ryx (row index first, like matrix).
        final double r10 = getCellValue(j, ix,   iy+1);
        final double r11 = getCellValue(j, ix+1, iy+1);
        if (!skipX) {
            double dx = r01 - r00;
            dx += (r11 - r10 - dx) * gridX;
            derivative.setElement(j, 0, derivative.getElement(j, 0) + dx);
        }
        if (!skipY) {
            double dy = r10 - r00;
            dy += (r11 - r01 - dy) * gridY;
            derivative.setElement(j, 1, derivative.getElement(j, 1) + dy);
        }
    }
    return derivative;
}
 

/**
 * Removes the sources dimensions that are not required for computing the target dimensions.
 * This method is invoked only if {@link #sourceDimensions} is non-null at {@link #separate()} invocation time.
 * This method can operate only on the first transform of a transformation chain.
 * If this method succeed, then {@link #sourceDimensions} will be updated.
 *
 * <p>This method can process only linear transforms (potentially indirectly through a concatenated transform).
 * Actually it would be possible to also process pass-through transform followed by a linear transform, but this
 * case should have been optimized during transform concatenation. If it is not the case, consider improving the
 * {@link PassThroughTransform#tryConcatenate(boolean, MathTransform, MathTransformFactory)} method instead then
 * this one.</p>
 *
 * @param  head  the first transform of a transformation chain.
 * @return the reduced transform, or {@code head} if this method did not reduced the transform.
 */
private MathTransform removeUnusedSourceDimensions(final MathTransform head) {
    Matrix m = MathTransforms.getMatrix(head);
    if (m != null) {
        int[] retainedDimensions = ArraysExt.EMPTY_INT;
        final int dimension = m.getNumCol() - 1;            // Number of source dimensions (ignore translations column).
        final int numRows   = m.getNumRow();                // Number of target dimensions + 1.
        for (int i=0; i<dimension; i++) {
            for (int j=0; j<numRows; j++) {
                if (m.getElement(j,i) != 0) {
                    // Found a source dimension which is required by target dimension.
                    final int length = retainedDimensions.length;
                    retainedDimensions = Arrays.copyOf(retainedDimensions, length+1);
                    retainedDimensions[length] = i;
                    break;
                }
            }
        }
        if (retainedDimensions.length != dimension) {
            /*
             * If we do not retain all dimensions, remove the matrix columns corresponding to the excluded
             * source dimensions and create a new transform. We remove consecutive columns in single calls
             * to 'removeColumns', from 'lower' inclusive to 'upper' exclusive.
             */
            int upper = dimension;
            for (int i = retainedDimensions.length; --i >= -1;) {
                final int keep = (i >= 0) ? retainedDimensions[i] : -1;
                final int lower = keep + 1;                                     // First column to exclude.
                if (lower != upper) {
                    // Remove source dimensions that are not retained.
                    m = MatrixSIS.castOrCopy(m).removeColumns(lower, upper);
                }
                upper = keep;
            }
            /*
             * If the user specified source dimensions, the indices need to be adjusted.
             * This loop has no effect if all source dimensions were kept before this method call.
             */
            for (int i=0; i<retainedDimensions.length; i++) {
                retainedDimensions[i] = sourceDimensions[retainedDimensions[i]];
            }
            sourceDimensions = retainedDimensions;
            return MathTransforms.linear(m);
        }
    } else if (head instanceof ConcatenatedTransform) {
        final MathTransform transform1 = ((ConcatenatedTransform) head).transform1;
        final MathTransform reduced = removeUnusedSourceDimensions(transform1);
        if (reduced != transform1) {
            return MathTransforms.concatenate(reduced, ((ConcatenatedTransform) head).transform2);
        }
    }
    return head;
}
 

/**
 * Sets all Bursa-Wolf parameters from the given <cite>Position Vector transformation</cite> matrix.
 * The matrix shall comply to the following constraints:
 *
 * <ul>
 *   <li>The matrix shall be {@linkplain org.apache.sis.referencing.operation.matrix.MatrixSIS#isAffine() affine}.</li>
 *   <li>The sub-matrix defined by {@code matrix} without the last row and last column shall be
 *       <a href="http://en.wikipedia.org/wiki/Skew-symmetric_matrix">skew-symmetric</a> (a.k.a. antisymmetric).</li>
 * </ul>
 *
 * @param  matrix     the matrix from which to get Bursa-Wolf parameters.
 * @param  tolerance  the tolerance error for the skew-symmetric matrix test, in units of PPM or arc-seconds (e.g. 1E-8).
 * @throws IllegalArgumentException if the specified matrix does not met the conditions.
 *
 * @see #getPositionVectorTransformation(Date)
 */
public void setPositionVectorTransformation(final Matrix matrix, final double tolerance) throws IllegalArgumentException {
    final int numRow = matrix.getNumRow();
    final int numCol = matrix.getNumCol();
    if (numRow != SIZE || numCol != SIZE) {
        final Integer n = SIZE;
        throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedMatrixSize_4, n, n, numRow, numCol));
    }
    if (!Matrices.isAffine(matrix)) {
        throw new IllegalArgumentException(Resources.format(Resources.Keys.NotAnAffineTransform));
    }
    /*
     * Translation terms, taken "as-is".
     * If the matrix contains only translation terms (which is often the case), we are done.
     */
    tX = matrix.getElement(0,3);
    tY = matrix.getElement(1,3);
    tZ = matrix.getElement(2,3);
    if (Matrices.isTranslation(matrix)) {                   // Optimization for a common case.
        return;
    }
    /*
     * Scale factor: take the average of elements on the diagonal. All those
     * elements should have the same value, but we tolerate slight deviation
     * (this will be verified later).
     */
    final DoubleDouble S = DoubleDouble.castOrCopy(getNumber(matrix, 0,0));
    S.addGuessError(getNumber(matrix, 1,1));
    S.addGuessError(getNumber(matrix, 2,2));
    S.divide(3);
    /*
     * Computes: RS = S * toRadians(1″)
     *           dS = (S-1) * PPM
     */
    final DoubleDouble RS = DoubleDouble.createSecondsToRadians();
    RS.multiply(S);
    S.add(-1);
    S.multiply(PPM);
    dS = S.doubleValue();
    /*
     * Rotation terms. Each rotation terms appear twice, with one value being the negative of the other value.
     * We verify this skew symmetric aspect in the loop. We also opportunistically verify that the scale terms
     * are uniform.
     */
    for (int j=0; j < SIZE-1; j++) {
        if (!(abs((matrix.getElement(j,j) - 1)*PPM - dS) <= tolerance)) {
            throw new IllegalArgumentException(Resources.format(Resources.Keys.NonUniformScale));
        }
        for (int i = j+1; i < SIZE-1; i++) {
            S.setFrom(RS);
            S.inverseDivideGuessError(getNumber(matrix, j,i));          // Negative rotation term.
            double value = S.value;
            double error = S.error;
            S.setFrom(RS);
            S.inverseDivideGuessError(getNumber(matrix, i,j));          // Positive rotation term.
            if (!(abs(value + S.value) <= tolerance)) {                 // We expect r1 ≈ -r2
                throw new IllegalArgumentException(Resources.format(Resources.Keys.NotASkewSymmetricMatrix));
            }
            S.subtract(value, error);
            S.multiply(0.5);
            value = S.doubleValue();                                    // Average of the two rotation terms.
            switch (j*SIZE + i) {
                case 1: rZ =  value; break;
                case 2: rY = -value; break;
                case 6: rX =  value; break;
            }
        }
    }
}
 

/**
 * Returns a matrix with the same content than the given matrix but a different size, assuming an affine transform.
 * This method can be invoked for adding or removing the <strong>last</strong> dimensions of an affine transform.
 * More specifically:
 *
 * <ul class="verbose">
 *   <li>If the given {@code numCol} is <var>n</var> less than the number of columns in the given matrix,
 *       then the <var>n</var> columns <em>before the last column</em> are removed.
 *       The last column is left unchanged because it is assumed to contain the translation terms.</li>
 *   <li>If the given {@code numCol} is <var>n</var> more than the number of columns in the given matrix,
 *       then <var>n</var> columns are inserted <em>before the last column</em>.
 *       All values in the new columns will be zero.</li>
 *   <li>If the given {@code numRow} is <var>n</var> less than the number of rows in the given matrix,
 *       then the <var>n</var> rows <em>before the last row</em> are removed.
 *       The last row is left unchanged because it is assumed to contain the usual [0 0 0 … 1] terms.</li>
 *   <li>If the given {@code numRow} is <var>n</var> more than the number of rows in the given matrix,
 *       then <var>n</var> rows are inserted <em>before the last row</em>.
 *       The corresponding offset and scale factors will be 0 and 1 respectively.
 *       In other words, new dimensions are propagated unchanged.</li>
 * </ul>
 *
 * @param  matrix  the matrix to resize. This matrix will never be changed.
 * @param  numRow  the new number of rows. This is equal to the desired number of target dimensions plus 1.
 * @param  numCol  the new number of columns. This is equal to the desired number of source dimensions plus 1.
 * @return a new matrix of the given size, or the given {@code matrix} if no resizing was needed.
 */
public static Matrix resizeAffine(final Matrix matrix, int numRow, int numCol) {
    ArgumentChecks.ensureNonNull         ("matrix", matrix);
    ArgumentChecks.ensureStrictlyPositive("numRow", numRow);
    ArgumentChecks.ensureStrictlyPositive("numCol", numCol);
    int srcRow = matrix.getNumRow();
    int srcCol = matrix.getNumCol();
    if (numRow == srcRow && numCol == srcCol) {
        return matrix;
    }
    final int      stride  = srcCol;
    final int      length  = srcCol * srcRow;
    final double[] sources = getExtendedElements(matrix);
    final DoubleDouble transfer = (sources.length > length) ? new DoubleDouble() : null;
    final MatrixSIS resized = createZero(numRow, numCol, transfer != null);
    final int copyRow = Math.min(--numRow, --srcRow);
    final int copyCol = Math.min(--numCol, --srcCol);
    for (int j=copyRow; j<numRow; j++) {
        resized.setElement(j, j, 1);
    }
    resized.setElements(sources, length, stride, transfer, 0,      0,       0,      0,       copyRow, copyCol);    // Shear and scale terms.
    resized.setElements(sources, length, stride, transfer, 0,      srcCol,  0,      numCol,  copyRow, 1);          // Translation column.
    resized.setElements(sources, length, stride, transfer, srcRow, 0,       numRow, 0,       1,       copyCol);    // Last row.
    resized.setElements(sources, length, stride, transfer, srcRow, srcCol,  numRow, numCol,  1,       1);          // Last row.
    return resized;
}
 

/**
 * If the given matrix to be concatenated to this transform, can be concatenated to the
 * sub-transform instead, returns the matrix to be concatenated to the sub-transform.
 * Otherwise returns {@code null}.
 *
 * <p>This method does not verify if the matrix size is compatible with this transform dimension.</p>
 *
 * @param  applyOtherFirst  {@code true} if the transformation order is {@code matrix} followed by {@code this}, or
 *                          {@code false} if the transformation order is {@code this} followed by {@code matrix}.
 */
private Matrix toSubMatrix(final boolean applyOtherFirst, final Matrix matrix) {
    final int numRow = matrix.getNumRow();
    final int numCol = matrix.getNumCol();
    if (numRow != numCol) {
        // Current implementation requires a square matrix.
        return null;
    }
    final int subDim = applyOtherFirst ? subTransform.getSourceDimensions()
                                       : subTransform.getTargetDimensions();
    final MatrixSIS sub = Matrices.createIdentity(subDim + 1);
    /*
     * Ensure that every dimensions which are scaled by the affine transform are one
     * of the dimensions modified by the sub-transform, and not any other dimension.
     */
    for (int j=numRow; --j>=0;) {
        final int sj = j - firstAffectedCoordinate;
        for (int i=numCol; --i>=0;) {
            final double element = matrix.getElement(j, i);
            if (sj >= 0 && sj < subDim) {
                final int si;
                final boolean copy;
                if (i == numCol-1) {                    // Translation term (last column)
                    si = subDim;
                    copy = true;
                } else {                                // Any term other than translation.
                    si = i - firstAffectedCoordinate;
                    copy = (si >= 0 && si < subDim);
                }
                if (copy) {
                    sub.setElement(sj, si, element);
                    continue;
                }
            }
            if (element != (i == j ? 1 : 0)) {
                // Found a dimension which perform some scaling or translation.
                return null;
            }
        }
    }
    return sub;
}