| /* |
| * Licensed to the Apache Software Foundation (ASF) under one or more |
| * contributor license agreements. See the NOTICE file distributed with |
| * this work for additional information regarding copyright ownership. |
| * The ASF licenses this file to You under the Apache License, Version 2.0 |
| * (the "License"); you may not use this file except in compliance with |
| * the License. You may obtain a copy of the License at |
| * |
| * http://www.apache.org/licenses/LICENSE-2.0 |
| * |
| * Unless required by applicable law or agreed to in writing, software |
| * distributed under the License is distributed on an "AS IS" BASIS, |
| * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. |
| * See the License for the specific language governing permissions and |
| * limitations under the License. |
| */ |
| package org.apache.sis.referencing.operation.transform; |
| |
| import java.util.Arrays; |
| import org.opengis.geometry.DirectPosition; |
| import org.opengis.referencing.operation.Matrix; |
| import org.apache.sis.referencing.privy.DirectPositionView; |
| import org.apache.sis.referencing.privy.ExtendedPrecisionMatrix; |
| import org.apache.sis.referencing.privy.Formulas; |
| import org.apache.sis.referencing.operation.matrix.Matrices; |
| import org.apache.sis.referencing.operation.matrix.MatrixSIS; |
| import org.apache.sis.referencing.internal.Arithmetic; |
| import org.apache.sis.util.ArgumentChecks; |
| import org.apache.sis.util.privy.Numerics; |
| import org.apache.sis.math.Fraction; |
| |
| |
| /** |
| * A usually affine, or otherwise a projective transform for the generic cases. |
| * This implementation is used for cases other than identity, translation only, |
| * scale only, 1D transform, 2D transform or axis swapping. |
| * |
| * <p>A projective transform is capable of mapping an arbitrary quadrilateral into another arbitrary quadrilateral, |
| * while preserving the straightness of lines. In the special case where the transform is affine, the parallelism of |
| * lines in the source is preserved in the output.</p> |
| * |
| * @author Martin Desruisseaux (IRD, Geomatys) |
| * |
| * @see java.awt.geom.AffineTransform |
| */ |
| class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPrecisionMatrix { |
| /** |
| * Serial number for inter-operability with different versions. |
| */ |
| private static final long serialVersionUID = -4813507361303377148L; |
| |
| /** |
| * The number of rows. |
| */ |
| private final int numRow; |
| |
| /** |
| * The number of columns. |
| */ |
| private final int numCol; |
| |
| /** |
| * Number of columns for coefficients other than scales and shears in the {@link #elt} array. |
| * There is one column for translation terms and one column for denominators. |
| */ |
| private static final int NON_SCALE_COLUMNS = 2; |
| |
| /** |
| * Elements of the matrix augmented with one column containing common denominators. |
| * Column indices vary fastest. |
| * |
| * <h4>Denominator column</h4> |
| * An additional column is appended after the translation column. |
| * That column contains a denominator inferred from fractional values found on the row. |
| * All elements in the matrix row shall be multiplied by that denominator. |
| * The intend is to increase the chances that matrix elements are integer values. |
| * If no fractional value is found, the default denominator value is 1. |
| */ |
| private final double[] elt; |
| |
| /** |
| * Same numbers as {@link #elt} excluding the denominators and with potentially extended precision. |
| * Zero values <em>shall</em> be represented by null elements. |
| */ |
| private final Number[] numbers; |
| |
| /** |
| * Constructs a transform from the specified matrix. |
| * The matrix is usually square and affine, but this is not enforced. |
| * Non-affine transforms (e.g. projective transforms) are accepted but may not be invertible. |
| * |
| * @param matrix the matrix containing the coefficients of this projective transform. |
| */ |
| protected ProjectiveTransform(final Matrix matrix) { |
| numRow = matrix.getNumRow(); |
| numCol = matrix.getNumCol(); |
| /* |
| * Get the matrix elements as `Number` instances if possible. |
| * Those instances allow better precision than `double` values. |
| * Those numbers are available only through SIS-specific API. |
| */ |
| final boolean hasNumbers; |
| if (matrix instanceof ExtendedPrecisionMatrix) { |
| // Use `writable = true` because we need a copy protected from changes. |
| numbers = ((ExtendedPrecisionMatrix) matrix).getElementAsNumbers(true); |
| hasNumbers = true; |
| } else if (matrix instanceof MatrixSIS) { |
| final MatrixSIS m = (MatrixSIS) matrix; |
| numbers = new Number[numRow * numCol]; |
| for (int i=0; i<numbers.length; i++) { |
| final Number e = m.getNumber(i / numCol, i % numCol); |
| if (!ExtendedPrecisionMatrix.isZero(e)) { |
| numbers[i] = e; |
| } |
| } |
| hasNumbers = true; |
| } else { |
| numbers = new Number[numRow * numCol]; |
| hasNumbers = false; |
| } |
| /* |
| * Get the matrix elements as `double` values through the standard matrix API. |
| * We do that as a way to preserve negative zero, which is lost in `numbers`. |
| * The `numbers` array is either completed or compared for consistency. |
| */ |
| final int dstDim = numRow - 1; // Last row is [0 0 … 1] in an affine transform. |
| final int rowStride = numCol + 1; // The `elt` array has an extra column for denominators. |
| elt = new double[numRow * rowStride - 1]; // We don't need the denominator of the [0 0 … 1] row. |
| for (int k=0,i=0; i < numRow; i++) { |
| for (int j=0; j < numCol; j++) { |
| final double e = matrix.getElement(i, j); |
| elt[k++] = e; // May be negative zero. |
| if (hasNumbers) { |
| assert epsilonEqual(e, numbers[i*numCol + j]); |
| } else if (e != 0) { |
| final int v = (int) e; // Check if we can store as integer. |
| numbers[i*numCol + j] = (v == e) ? Integer.valueOf(v) : Double.valueOf(e); |
| } |
| } |
| if (i != dstDim) { |
| elt[k++] = 1; |
| } else { |
| assert k == elt.length; |
| } |
| } |
| /* |
| * At this point, this `ProjectiveTransform` is initialized and valid. |
| * Optionally update the elements values for reducing rounding errors |
| * when a denominator can be identified. This is where the denominator |
| * column in the `elt` array may get values different than 1. |
| */ |
| if (hasNumbers) { |
| for (int row=0; row < dstDim; row++) { |
| final int lower = numCol * row; |
| final int upper = numCol + lower; |
| final Integer denominator; |
| try { |
| Fraction sum = null; |
| for (int i=lower; i<upper; i++) { |
| final Number element = numbers[i]; |
| if (element instanceof Fraction) { |
| final Fraction f = (Fraction) element; |
| sum = (sum != null) ? sum.add(f) : f; |
| } |
| } |
| if (sum == null) { |
| continue; |
| } |
| denominator = sum.denominator; |
| } catch (ArithmeticException e) { |
| continue; |
| } |
| int k = row * rowStride; |
| for (int i=lower; i<upper; i++) { |
| final Number element = Arithmetic.multiply(numbers[i], denominator); |
| if (element != null) { |
| elt[k] = element.doubleValue(); |
| } |
| k++; |
| } |
| elt[k] = denominator.doubleValue(); |
| } |
| } |
| } |
| |
| /** |
| * Returns whether the given number are equal, with a tolerance of 1 ULP. |
| * A null {@code Number} is interpreted as zero. |
| */ |
| private static boolean epsilonEqual(final double e, final Number v) { |
| return Numerics.epsilonEqual(e, (v != null) ? v.doubleValue() : 0, Math.ulp(e)); |
| } |
| |
| /** |
| * If a more efficient implementation of this math transform can be used, returns it. |
| * Otherwise returns {@code this} unchanged. |
| */ |
| final LinearTransform optimize() { |
| if (numCol < numRow) { |
| return this; |
| } |
| final int n = (numRow - 1) * numCol; |
| for (int i = 0; i < numCol;) { |
| if (!isIdentity(numbers[n + i], ++i == numCol)) { |
| return this; // Transform is not affine (ignoring if square or not). |
| } |
| } |
| /* |
| * Note: we could check for CopyTransform case here, but this check is rather done in |
| * MathTransforms.linear(Matrix) in order to avoid ProjectiveTransform instantiation. |
| * |
| * Check which elements are non-zero. For ScaleTransform, they must be on the diagonal. |
| * For TranslationTransform, they must be in the last column. Note that the transform |
| * should not be identity (except for testing purpose) since the identity case should |
| * have been handled by MathTransforms.linear(Matrix) before to reach this point. |
| */ |
| boolean isScale = true; // ScaleTransform accepts non-square matrix. |
| boolean isTranslation = (numRow == numCol); // TranslationTransform is restricted to square matrix. |
| final int lastColumn = numCol - 1; |
| for (int i=0; i<n; i++) { |
| final Number element = numbers[i]; |
| final int row = (i / numCol); |
| final int col = (i % numCol); |
| if (col != row) isScale &= (element == null); |
| if (col != lastColumn) isTranslation &= isIdentity(element, col == row); |
| if (!(isScale | isTranslation)) { |
| return this; |
| } |
| } |
| if (isTranslation) { |
| return new TranslationTransform(numRow, numbers); |
| } else { |
| return new ScaleTransform(numRow, numCol, numbers); |
| } |
| } |
| |
| /** |
| * Returns {@code true} if the given element has the expected value of an identity matrix. |
| * |
| * @param element the element to test. |
| * @param diagonal whether the element is on the diagonal. |
| * @return whether the given element is an element of an identity matrix. |
| * |
| * @see #isIdentity() |
| */ |
| private static boolean isIdentity(final Number element, final boolean diagonal) { |
| return diagonal ? (element != null && element.doubleValue() == 1) : (element == null); |
| } |
| |
| /** |
| * Gets the dimension of input points. |
| */ |
| @Override |
| public final int getSourceDimensions() { |
| return numCol - 1; |
| } |
| |
| /** |
| * Gets the dimension of output points. |
| */ |
| @Override |
| public final int getTargetDimensions() { |
| return numRow - 1; |
| } |
| |
| /** |
| * Gets the number of rows in the matrix. |
| */ |
| @Override |
| public final int getNumRow() { |
| return numRow; |
| } |
| |
| /** |
| * Gets the number of columns in the matrix. |
| */ |
| @Override |
| public final int getNumCol() { |
| return numCol; |
| } |
| |
| /** |
| * Returns a copy of all matrix elements in a flat, row-major array. Zero values <em>shall</em> be null. |
| * Callers can write in the returned array if and only if the {@code writable} argument is {@code true}. |
| */ |
| @Override |
| public final Number[] getElementAsNumbers(final boolean writable) { |
| return writable ? numbers.clone() : numbers; |
| } |
| |
| /** |
| * Retrieves the value at the specified row and column of the matrix. |
| * If the value is zero, then this method <em>shall</em> return {@code null}. |
| */ |
| @Override |
| public final Number getElementOrNull(final int row, final int column) { |
| ArgumentChecks.ensureBetween("row", 0, numRow - 1, row); |
| ArgumentChecks.ensureBetween("column", 0, numCol - 1, column); |
| return numbers[row * numCol + column]; |
| } |
| |
| /** |
| * Returns the matrix element at the given index. |
| */ |
| @Override |
| public final double getElement(final int row, final int column) { |
| ArgumentChecks.ensureBetween("row", 0, numRow - 1, row); |
| ArgumentChecks.ensureBetween("column", 0, numCol - 1, column); |
| final Number element = numbers[row * numCol + column]; |
| if (element != null) { |
| return element.doubleValue(); |
| } |
| /* |
| * Fallback on the `elt` array only for 0 values for avoiding the need to divide by the denominator. |
| * Do not return a hard-coded 0 value in order to preserve the sign of negative zero. |
| */ |
| final int rowStride = numCol + 1; |
| return elt[row * rowStride + column]; |
| } |
| |
| /** |
| * Tests whether this transform does not move any points. |
| * |
| * <h4>Note</h4> |
| * This method should always returns {@code false}, because |
| * {@code MathTransforms.linear(…)} should have created specialized implementations for identity cases. |
| * Nevertheless we perform the full check as a safety, in case someone instantiated this class directly |
| * instead of using a factory method. |
| */ |
| @Override |
| public final boolean isIdentity() { |
| if (numRow != numCol) { |
| return false; |
| } |
| for (int i=0; i < numbers.length; i++) { |
| if (!isIdentity(numbers[i], (i / numCol) == (i % numCol))) { |
| return false; |
| } |
| } |
| return true; |
| } |
| |
| /** |
| * Converts a single coordinate tuple in a list of coordinate values, |
| * and optionally computes the derivative at that location. |
| * |
| * @return {@inheritDoc} |
| */ |
| @Override |
| public final Matrix transform(final double[] srcPts, final int srcOff, |
| final double[] dstPts, final int dstOff, |
| final boolean derivate) |
| { |
| if (!derivate) { |
| transform(srcPts, srcOff, dstPts, dstOff, 1); |
| return null; |
| } |
| // A non-null position is required in case this transform is non-affine. |
| Matrix derivative = derivative(new DirectPositionView.Double(srcPts, srcOff, getSourceDimensions())); |
| transform(srcPts, srcOff, dstPts, dstOff, 1); |
| return derivative; |
| } |
| |
| /** |
| * Transforms an array of floating point coordinates by this matrix. Point coordinates must have a dimension |
| * equal to <code>{@link Matrix#getNumCol}-1</code>. For example, for square matrix of size 4×4, coordinate |
| * points are three-dimensional and stored in the arrays starting at the specified offset ({@code srcOff}) in |
| * the order |
| * <code>[x₀, y₀, z₀, |
| * x₁, y₁, z₁..., |
| * x<sub>n</sub>, y<sub>n</sub>, z<sub>n</sub>]</code>. |
| * |
| * @param srcPts the array containing the source point coordinates. |
| * @param srcOff the offset to the first point to be transformed in the source array. |
| * @param dstPts the array into which the transformed point coordinates are returned. |
| * @param dstOff the offset to the location of the first transformed point that is stored in the |
| * destination array. The source and destination array sections can overlap. |
| * @param numPts the number of points to be transformed. |
| */ |
| @Override |
| public final void transform(double[] srcPts, int srcOff, final double[] dstPts, int dstOff, int numPts) { |
| final int srcDim, dstDim; |
| int srcInc = srcDim = numCol - 1; // The last coordinate will be assumed equal to 1. |
| int dstInc = dstDim = numRow - 1; |
| if (srcPts == dstPts) { |
| switch (IterationStrategy.suggest(srcOff, srcDim, dstOff, dstDim, numPts)) { |
| case ASCENDING: { |
| break; |
| } |
| case DESCENDING: { |
| srcOff += (numPts - 1) * srcDim; |
| dstOff += (numPts - 1) * dstDim; |
| srcInc = -srcInc; |
| dstInc = -dstInc; |
| break; |
| } |
| default: { |
| srcPts = Arrays.copyOfRange(srcPts, srcOff, srcOff + numPts*srcDim); |
| srcOff = 0; |
| break; |
| } |
| } |
| } |
| final double[] buffer = new double[numRow]; |
| while (--numPts >= 0) { |
| int mix = 0; |
| for (int j=0; j<numRow; j++) { |
| double sum = elt[mix + srcDim]; // Initialize to translation term. |
| for (int i=0; i<srcDim; i++) { |
| final double e = elt[mix++]; |
| if (e != 0) { |
| /* |
| * The purpose of the test for non-zero value is not performance (it is actually more likely |
| * to slow down the calculation), but to get a valid sum even if some source coordinates are NaN. |
| * This occurs when the ProjectiveTransform is used for excluding some dimensions, for example |
| * getting 2D points from 3D points. In such case, the fact that the excluded dimensions had |
| * NaN values should not force the retained dimensions to get NaN values. |
| */ |
| if (Formulas.USE_FMA) { |
| sum = Math.fma(srcPts[srcOff + i], e, sum); |
| } else { |
| sum += srcPts[srcOff + i] * e; |
| } |
| } |
| } |
| buffer[j] = sum; |
| mix += NON_SCALE_COLUMNS; // Skip the translation column and the denominator column. |
| } |
| int k = numCol; |
| final int rowStride = numCol + 1; |
| final double w = buffer[dstDim]; |
| for (int j=0; j<dstDim; j++) { |
| // `w` is equal to 1 if the transform is affine. |
| dstPts[dstOff + j] = buffer[j] / (w * elt[k]); |
| k += rowStride; |
| } |
| srcOff += srcInc; |
| dstOff += dstInc; |
| } |
| } |
| |
| /** |
| * Transforms an array of floating point coordinates by this matrix. Point coordinates must have a dimension |
| * equal to <code>{@link Matrix#getNumCol()} - 1</code>. For example, for square matrix of size 4×4, coordinate |
| * points are three-dimensional and stored in the arrays starting at the specified offset ({@code srcOff}) |
| * in the order |
| * <code>[x₀, y₀, z₀, |
| * x₁, y₁, z₁..., |
| * x<sub>n</sub>, y<sub>n</sub>, z<sub>n</sub>]</code>. |
| * |
| * @param srcPts the array containing the source point coordinates. |
| * @param srcOff the offset to the first point to be transformed in the source array. |
| * @param dstPts the array into which the transformed point coordinates are returned. |
| * @param dstOff the offset to the location of the first transformed point that is stored in the |
| * destination array. The source and destination array sections can overlap. |
| * @param numPts the number of points to be transformed. |
| */ |
| @Override |
| public final void transform(float[] srcPts, int srcOff, final float[] dstPts, int dstOff, int numPts) { |
| final int srcDim, dstDim; |
| int srcInc = srcDim = numCol - 1; |
| int dstInc = dstDim = numRow - 1; |
| if (srcPts == dstPts) { |
| switch (IterationStrategy.suggest(srcOff, srcDim, dstOff, dstDim, numPts)) { |
| case ASCENDING: { |
| break; |
| } |
| case DESCENDING: { |
| srcOff += (numPts - 1) * srcDim; |
| dstOff += (numPts - 1) * dstDim; |
| srcInc = -srcInc; |
| dstInc = -dstInc; |
| break; |
| } |
| default: { |
| srcPts = Arrays.copyOfRange(srcPts, srcOff, srcOff + numPts*srcDim); |
| srcOff = 0; |
| break; |
| } |
| } |
| } |
| final double[] buffer = new double[numRow]; |
| while (--numPts >= 0) { |
| int mix = 0; |
| for (int j=0; j<numRow; j++) { |
| double sum = elt[mix + srcDim]; |
| for (int i=0; i<srcDim; i++) { |
| final double e = elt[mix++]; |
| if (e != 0) { // See comment in transform(double[], ...) |
| if (Formulas.USE_FMA) { |
| sum = Math.fma(srcPts[srcOff + i], e, sum); |
| } else { |
| sum += srcPts[srcOff + i] * e; |
| } |
| } |
| } |
| buffer[j] = sum; |
| mix += NON_SCALE_COLUMNS; |
| } |
| int k = numCol; |
| final int rowStride = numCol + 1; |
| final double w = buffer[dstDim]; |
| for (int j=0; j<dstDim; j++) { |
| dstPts[dstOff + j] = (float) (buffer[j] / (w * elt[k])); |
| k += rowStride; |
| } |
| srcOff += srcInc; |
| dstOff += dstInc; |
| } |
| } |
| |
| /** |
| * Transforms an array of floating point coordinates by this matrix. |
| * |
| * @param srcPts the array containing the source point coordinates. |
| * @param srcOff the offset to the first point to be transformed in the source array. |
| * @param dstPts the array into which the transformed point coordinates are returned. |
| * @param dstOff the offset to the location of the first transformed point that is stored in the destination array. |
| * @param numPts the number of points to be transformed. |
| */ |
| @Override |
| public final void transform(final double[] srcPts, int srcOff, final float[] dstPts, int dstOff, int numPts) { |
| final int srcDim = numCol - 1; |
| final int dstDim = numRow - 1; |
| final double[] buffer = new double[numRow]; |
| while (--numPts >= 0) { |
| int mix = 0; |
| for (int j=0; j<numRow; j++) { |
| double sum = elt[mix + srcDim]; |
| for (int i=0; i<srcDim; i++) { |
| final double e = elt[mix++]; |
| if (e != 0) { // See comment in transform(double[], ...) |
| if (Formulas.USE_FMA) { |
| sum = Math.fma(srcPts[srcOff + i], e, sum); |
| } else { |
| sum += srcPts[srcOff + i] * e; |
| } |
| } |
| } |
| buffer[j] = sum; |
| mix += NON_SCALE_COLUMNS; |
| } |
| int k = numCol; |
| final int rowStride = numCol + 1; |
| final double w = buffer[dstDim]; |
| for (int j=0; j<dstDim; j++) { |
| dstPts[dstOff++] = (float) (buffer[j] / (w * elt[k])); |
| k += rowStride; |
| } |
| srcOff += srcDim; |
| } |
| } |
| |
| /** |
| * Transforms an array of floating point coordinates by this matrix. |
| * |
| * @param srcPts the array containing the source point coordinates. |
| * @param srcOff the offset to the first point to be transformed in the source array. |
| * @param dstPts the array into which the transformed point coordinates are returned. |
| * @param dstOff the offset to the location of the first transformed point that is stored in the destination array. |
| * @param numPts the number of points to be transformed. |
| */ |
| @Override |
| public final void transform(final float[] srcPts, int srcOff, final double[] dstPts, int dstOff, int numPts) { |
| final int srcDim = numCol - 1; |
| final int dstDim = numRow - 1; |
| final double[] buffer = new double[numRow]; |
| while (--numPts >= 0) { |
| int mix = 0; |
| for (int j=0; j<numRow; j++) { |
| double sum = elt[mix + srcDim]; |
| for (int i=0; i<srcDim; i++) { |
| final double e = elt[mix++]; |
| if (e != 0) { // See comment in transform(double[], ...) |
| if (Formulas.USE_FMA) { |
| sum = Math.fma(srcPts[srcOff + i], e, sum); |
| } else { |
| sum += srcPts[srcOff + i] * e; |
| } |
| } |
| } |
| buffer[j] = sum; |
| mix += NON_SCALE_COLUMNS; |
| } |
| int k = numCol; |
| final int rowStride = numCol + 1; |
| final double w = buffer[dstDim]; |
| for (int j=0; j<dstDim; j++) { |
| dstPts[dstOff++] = buffer[j] / (w * elt[k]); |
| k += rowStride; |
| } |
| srcOff += srcDim; |
| } |
| } |
| |
| /** |
| * Gets the derivative of this transform at a point. For an affine transform, |
| * the derivative is the same everywhere and the given point is ignored. |
| * For a perspective transform, the given point is used. |
| * |
| * @param point the coordinate tuple where to evaluate the derivative. |
| */ |
| @Override |
| public final Matrix derivative(final DirectPosition point) { |
| final int srcDim = numCol - 1; |
| final int dstDim = numRow - 1; |
| final int rowStride = numCol + 1; |
| /* |
| * In the `transform(…)` method, all coordinate values are divided by a `w` coefficient |
| * which depends on the position. We need to reproduce that division here. Note that `w` |
| * coefficient is different than 1 only if the transform is non-affine. |
| */ |
| int mix = dstDim * rowStride; |
| double w = elt[mix + srcDim]; // `w` is equal to 1 if the transform is affine. |
| for (int i=0; i<srcDim; i++) { |
| final double e = elt[mix++]; |
| if (e != 0) { // For avoiding NullPointerException if affine. |
| w += point.getCoordinate(i) * e; |
| } |
| } |
| /* |
| * In the usual affine case (w=1), the derivative is a copy of this matrix |
| * with last row and last column omitted. |
| */ |
| mix = 0; |
| int k = numCol; |
| final MatrixSIS matrix = Matrices.createZero(dstDim, srcDim); |
| for (int j=0; j<dstDim; j++) { |
| final double r = w * elt[k]; |
| for (int i=0; i<srcDim; i++) { |
| matrix.setElement(j, i, elt[mix++] / r); |
| } |
| mix += NON_SCALE_COLUMNS; // Skip the translation column and the denominator column. |
| k += rowStride; |
| } |
| return matrix; |
| } |
| |
| /** |
| * {@inheritDoc} |
| * |
| * @return {@inheritDoc} |
| */ |
| @Override |
| protected int computeHashCode() { |
| return Arrays.hashCode(elt) + 31 * super.computeHashCode(); |
| } |
| |
| /** |
| * Compares this math transform with an object which is known to be an instance of the same class. |
| */ |
| @Override |
| protected boolean equalsSameClass(final Object object) { |
| final ProjectiveTransform that = (ProjectiveTransform) object; |
| return numRow == that.numRow && |
| numCol == that.numCol && |
| Arrays.equals(elt, that.elt) && |
| Arrays.equals(numbers, that.numbers); |
| } |
| } |
