You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@sis.apache.org by de...@apache.org on 2023/01/17 14:30:03 UTC
[sis] 01/02: Opportunistically use the division by `w` in ProjectiveTransform for reducing rounding errors with fractional matrix element values.
This is an automated email from the ASF dual-hosted git repository.
desruisseaux pushed a commit to branch geoapi-3.1
in repository https://gitbox.apache.org/repos/asf/sis.git
commit 3014fe6003c66cc4c54a3fc144c3cc339a67e191
Author: Martin Desruisseaux <ma...@geomatys.com>
AuthorDate: Tue Jan 17 15:19:31 2023 +0100
Opportunistically use the division by `w` in ProjectiveTransform
for reducing rounding errors with fractional matrix element values.
---
.../transform/AbstractLinearTransform.java | 4 +-
.../operation/transform/ProjectiveTransform.java | 221 ++++++++++++++++-----
.../transform/ProjectiveTransformTest.java | 12 ++
3 files changed, 182 insertions(+), 55 deletions(-)
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
index f37950fc84..f8cc1d777e 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/AbstractLinearTransform.java
@@ -80,8 +80,8 @@ abstract class AbstractLinearTransform extends AbstractMathTransform implements
for (int i=0; i<elements.length; i++) {
final double element = elements[i];
if (element != 0) {
- final int vi = (int) element; // Check if we can store as integer.
- numbers[i] = (vi == element) ? Integer.valueOf(vi) : Double.valueOf(element);
+ final int ie = (int) element; // Check if we can store as integer.
+ numbers[i] = (ie == element) ? Integer.valueOf(ie) : Double.valueOf(element);
}
}
return numbers;
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
index 5188c12e66..0e6f023b0c 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
@@ -23,9 +23,11 @@ import org.apache.sis.internal.referencing.DirectPositionView;
import org.apache.sis.referencing.operation.matrix.Matrices;
import org.apache.sis.referencing.operation.matrix.MatrixSIS;
import org.apache.sis.internal.referencing.ExtendedPrecisionMatrix;
+import org.apache.sis.internal.referencing.Arithmetic;
import org.apache.sis.internal.referencing.Formulas;
-import org.apache.sis.util.resources.Errors;
+import org.apache.sis.internal.util.Numerics;
import org.apache.sis.util.ArgumentChecks;
+import org.apache.sis.math.Fraction;
/**
@@ -48,7 +50,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
/**
* Serial number for inter-operability with different versions.
*/
- private static final long serialVersionUID = -5616239625370371256L;
+ private static final long serialVersionUID = -4813507361303377148L;
/**
* The number of rows.
@@ -61,12 +63,26 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
private final int numCol;
/**
- * Elements of the matrix. Column indices vary fastest.
+ * 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 than {@link #elt} with potentially extended precision.
+ * Same numbers than {@link #elt} excluding the denominators and with potentially extended precision.
* Zero values <em>shall</em> be represented by null elements.
*/
private final Number[] numbers;
@@ -74,47 +90,110 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
/**
* 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.
+ * @param matrix the matrix containing the coefficients of this projective transform.
*/
protected ProjectiveTransform(final Matrix matrix) {
numRow = matrix.getNumRow();
numCol = matrix.getNumCol();
- Number[] numbers = null;
+ /*
+ * 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);
- }
- if (matrix instanceof MatrixSIS) {
+ hasNumbers = true;
+ } else if (matrix instanceof MatrixSIS) {
final MatrixSIS m = (MatrixSIS) matrix;
- elt = m.getElements();
- if (elt.length != numRow * numCol) {
- throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedArrayLengths));
- }
- if (numbers == null) {
- numbers = new Number[elt.length];
- for (int i=0; i<numbers.length; i++) {
- final Number element = m.getNumber(i / numCol, i % numCol);
- if (!ExtendedPrecisionMatrix.isZero(element)) {
- numbers[i] = element;
- }
+ 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 {
- elt = new double[numRow * numCol];
- for (int i=0; i<elt.length; i++) {
- elt[i] = matrix.getElement(i / numCol, i % numCol);
+ 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 (numbers == null) {
- numbers = wrap(elt);
+ if (i != dstDim) {
+ elt[k++] = 1;
+ } else {
+ assert k == elt.length;
}
}
- this.numbers = numbers;
- if (elt.length != numbers.length) {
- throw new IllegalArgumentException(Errors.format(Errors.Keys.MismatchedArrayLengths));
+ /*
+ * 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.
@@ -125,7 +204,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
}
final int n = (numRow - 1) * numCol;
for (int i = 0; i < numCol;) {
- if (elt[n + i] != (++i == numCol ? 1 : 0)) {
+ if (!isIdentity(numbers[n + i], ++i == numCol)) {
return this; // Transform is not affine (ignoring if square or not).
}
}
@@ -142,11 +221,11 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
boolean isTranslation = (numRow == numCol); // TranslationTransform is restricted to square matrix.
final int lastColumn = numCol - 1;
for (int i=0; i<n; i++) {
- final double v = elt[i];
+ final Number element = numbers[i];
final int row = (i / numCol);
final int col = (i % numCol);
- isScale &= (col == row) || (v == 0);
- isTranslation &= (col == lastColumn) || (v == (col == row ? 1 : 0));
+ if (col != row) isScale &= (element == null);
+ if (col != lastColumn) isTranslation &= isIdentity(element, col == row);
if (!(isScale | isTranslation)) {
return this;
}
@@ -158,6 +237,19 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
}
}
+ /**
+ * 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.
*/
@@ -217,28 +309,35 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
public final double getElement(final int row, final int column) {
ArgumentChecks.ensureBetween("row", 0, numRow - 1, row);
ArgumentChecks.ensureBetween("column", 0, numCol - 1, column);
- return elt[row * numCol + 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.
*
- * <div class="note"><b>Note:</b> this method should always returns {@code false}, since
+ * <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.</div>
+ * instead of using a factory method.
*/
@Override
public final boolean isIdentity() {
if (numRow != numCol) {
return false;
}
- int mix = 0;
- for (int j=0; j<numRow; j++) {
- for (int i=0; i<numCol; i++) {
- if (elt[mix++] != (i == j ? 1 : 0)) {
- return false;
- }
+ for (int i=0; i < numbers.length; i++) {
+ if (!isIdentity(numbers[i], (i / numCol) == (i % numCol))) {
+ return false;
}
}
return true;
@@ -328,12 +427,15 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
}
}
buffer[j] = sum;
- mix++;
+ 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;
+ dstPts[dstOff + j] = buffer[j] / (w * elt[k]);
+ k += rowStride;
}
srcOff += srcInc;
dstOff += dstInc;
@@ -396,11 +498,14 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
}
}
buffer[j] = sum;
- mix++;
+ 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);
+ dstPts[dstOff + j] = (float) (buffer[j] / (w * elt[k]));
+ k += rowStride;
}
srcOff += srcInc;
dstOff += dstInc;
@@ -436,11 +541,14 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
}
}
buffer[j] = sum;
- mix++;
+ 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);
+ dstPts[dstOff++] = (float) (buffer[j] / (w * elt[k]));
+ k += rowStride;
}
srcOff += srcDim;
}
@@ -475,11 +583,14 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
}
}
buffer[j] = sum;
- mix++;
+ 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;
+ dstPts[dstOff++] = buffer[j] / (w * elt[k]);
+ k += rowStride;
}
srcOff += srcDim;
}
@@ -494,14 +605,15 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
*/
@Override
public final Matrix derivative(final DirectPosition point) {
- final int srcDim = numCol - 1;
- final int dstDim = numRow - 1;
+ 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 * numCol;
+ 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++];
@@ -514,12 +626,15 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
* 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++] / w);
+ matrix.setElement(j, i, elt[mix++] / r);
}
- mix++; // Skip translation column.
+ mix += NON_SCALE_COLUMNS; // Skip the translation column and the denominator column.
+ k += rowStride;
}
return matrix;
}
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java
index 7ed92c62a3..db8f374a07 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/ProjectiveTransformTest.java
@@ -177,6 +177,13 @@ public class ProjectiveTransformTest extends AffineTransformTest {
* Tests the concatenation of transforms that would result in rounding errors
* in extended-precision matrix operations were not used.
*
+ * Actually there is two sources of rounding errors tested by this method.
+ * The first source is rounding errors caused by matrix multiplications.
+ * The other source is rounding errors inside the {@code transform(…)} methods,
+ * which is reduced by a denominator column in {@link ProjectiveTransform#elt}.
+ * For demonstrating the latter rounding errors, it may be necessary to set the
+ * {@link org.apache.sis.internal.referencing.Formulas#USE_FMA} flag to {@code false}.
+ *
* @throws FactoryException if the transform cannot be created.
* @throws TransformException if a coordinate conversion failed.
*
@@ -204,6 +211,11 @@ public class ProjectiveTransformTest extends AffineTransformTest {
assertEquals(new Fraction(2, 37), m.getElementOrNull(0,0));
assertEquals(DoubleDouble.of(325).divide(100000), m.getElementOrNull(1,1));
assertEquals(new Fraction(-17, 127), m.getElementOrNull(2,2));
+ assertNull ( m.getElementOrNull(0,1));
+ assertEquals( 0, m.getElement(0,1), STRICT);
+ assertEquals( 2d / 37, m.getElement(0,0), 1E-15);
+ assertEquals( 0.00325, m.getElement(1,1), 1E-15);
+ assertEquals(-17d / 127, m.getElement(2,2), 1E-15);
/*
* Test a transformation, expecting exact result.
*/