This is an automated email from the ASF dual-hosted git repository. desruisseaux pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit 6092c0acb78a807021d15d12836cce12105fb6dc Author: Martin Desruisseaux <[email protected]> AuthorDate: Fri Apr 12 20:21:39 2019 +0200 First draft of an InterpolatedTransform algorithm using derivative (Jacobian matrix) for deciding in which direction to compensate for errors. Not yet enabled. --- .../provider/DatumShiftGridCompressed.java | 35 +++- .../sis/referencing/datum/DatumShiftGrid.java | 72 ++++++-- .../operation/transform/InterpolatedTransform.java | 140 +++++++++------ .../provider/DatumShiftGridCompressedTest.java | 46 +++++ .../provider/DatumShiftGridFileTest.java | 194 +++++++++++++++++++++ .../transform/InterpolatedTransformTest.java | 167 ++++++++++++------ .../operation/transform/MathTransformTestCase.java | 27 ++- .../operation/transform/QuadraticShiftGrid.java | 150 ++++++++++++++++ .../sis/test/suite/ReferencingTestSuite.java | 2 + 9 files changed, 690 insertions(+), 143 deletions(-) diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java index 55ef3b2..97cef2e 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressed.java @@ -196,15 +196,34 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T iy = ymax; gridY = 1; } - final int p0 = nx*iy + ix; - final int p1 = nx + p0; - for (int dim = 0; dim < data.length; dim++) { + final int p00 = nx*iy + ix; + final int p10 = nx + p00; + final int n = data.length; + boolean derivative = (vector.length >= n + 4); + for (int dim = 0; dim < n; dim++) { + double dx; final short[] values = data[dim]; - double r0 = values[p0]; - double r1 = values[p1]; - r0 += gridX * (values[p0+1] - r0); - r1 += gridX * (values[p1+1] - r1); - vector[dim] = (gridY * (r1 - r0) + r0) * scale + averages[dim]; + final double r00 = values[p00 ]; + final double r01 = values[p00 + 1]; // Naming convention: ryx (row index first, like matrix). + final double r10 = values[p10 ]; + final double r11 = values[p10 + 1]; + final double r0x = r00 + gridX * (dx = r01 - r00); // TODO: use Math.fma on JDK9. + final double r1x = r10 + gridX * ( r11 - r10); + vector[dim] = (gridY * (r1x - r0x) + r0x) * scale + averages[dim]; + if (derivative) { + double dy = (r10 - r00) * scale; + dx *= scale; + int i = n; + if (dim == 0) { + dx++; + } else { + dy++; + i += 2; + derivative = false; + } + vector[i ] = dx; + vector[i+1] = dy; + } } } diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java index 1f77e28..1c3cdd3 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/datum/DatumShiftGrid.java @@ -141,6 +141,12 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T private static final long serialVersionUID = 8405276545243175808L; /** + * Number of dimensions in which interpolations are applied. The grid may have more dimensions, + * but only this number of dimensions will be used in interpolations. + */ + private static final int INTERPOLATED_DIMENSIONS = 2; + + /** * The unit of measurements of input values, before conversion to grid indices by {@link #coordinateToGrid}. * The coordinate unit is typically {@link org.apache.sis.measure.Units#DEGREE}. * @@ -174,8 +180,8 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T private final boolean isCellValueRatio; /** - * Number of grid cells in along each dimension. This is usually an array of length 2 containing - * the number of grid cells along the <var>x</var> and <var>y</var> axes. + * Number of grid cells along each dimension. This is usually an array of length {@value #INTERPOLATED_DIMENSIONS} + * containing the number of grid cells along the <var>x</var> and <var>y</var> axes. */ private final int[] gridSize; @@ -232,8 +238,8 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T /** * Computes the conversion factors needed by {@link #interpolateInCell(double, double, double[])}. - * This method takes only the 2 first dimensions. If a conversion factor can not be computed, - * then it is set to NaN. + * This method takes only the {@value #INTERPOLATED_DIMENSIONS} first dimensions. If a conversion + * factor can not be computed, then it is set to NaN. */ @SuppressWarnings("fallthrough") private void computeConversionFactors() { @@ -466,10 +472,24 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T * the given {@code vector} array, which shall have a length of at least {@link #getTranslationDimensions()}. * The output unit of measurement is the same than the one documented in {@link #getCellValue(int, int, int)}. * - * <p>If the given coordinates are outside this grid, then this method computes the translation vector at the + * <div class="section">Extrapolations</div> + * If the given coordinates are outside this grid, then this method computes the translation vector at the * closest position in the grid. Applying translations on points outside the grid is a kind of extrapolation, * but some amount of extrapolations are necessary for operations like transforming an envelope before to compute - * its intersection with another envelope.</p> + * its intersection with another envelope. + * + * <div class="section">Derivative (Jacobian matrix)</div> + * If the length of the given array is at least <var>n</var> + 4 where <var>n</var> = {@link #getTranslationDimensions()}, + * then this method appends the derivative (approximated) at the given grid indices. This is the same derivative than the + * one computed by {@link #derivativeInCell(double, double)}, opportunistically computed here for performance reasons. + * The matrix layout is as below, where <var>u</var> and <var>v</var> are the coordinates after translation. + * + * {@preformat math + * ┌ ┐ ┌ ┐ + * │ ∂u/∂x ∂u/∂y │ = │ vector[n+0] vector[n+1] │ + * │ ∂v/∂x ∂v/∂y │ │ vector[n+2] vector[n+3] │ + * └ ┘ └ ┘ + * } * * <div class="section">Default implementation</div> * The default implementation performs the following steps for each dimension <var>dim</var>, @@ -480,6 +500,7 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T * <li>Clamp the {@code gridY} index into the [0 … {@code gridSize[1]} - 1] range, inclusive.</li> * <li>Using {@link #getCellValue}, get the cell values around the given indices.</li> * <li>Apply a bilinear interpolation and store the result in {@code vector[dim]}.</li> + * <li>Appends Jacobian matrix coefficients if the {@code vector} length is sufficient.</li> * </ol> * * @param gridX first grid coordinate of the point for which to get the translation. @@ -503,12 +524,33 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T gridY = 1; } n = getTranslationDimensions(); + boolean derivative = (vector.length >= n + INTERPOLATED_DIMENSIONS * INTERPOLATED_DIMENSIONS); for (int dim = 0; dim < n; dim++) { - double r0 = getCellValue(dim, ix, iy ); - double r1 = getCellValue(dim, ix, iy+1); - r0 += gridX * (getCellValue(dim, ix+1, iy ) - r0); - r1 += gridX * (getCellValue(dim, ix+1, iy+1) - r1); - vector[dim] = gridY * (r1 - r0) + r0; + double dx; + final double r00 = getCellValue(dim, ix, iy ); + final double r01 = getCellValue(dim, ix+1, iy ); // Naming convention: ryx (row index first, like matrix). + final double r10 = getCellValue(dim, ix, iy+1); + final double r11 = getCellValue(dim, ix+1, iy+1); + final double r0x = r00 + gridX * (dx = r01 - r00); // TODO: use Math.fma on JDK9. + final double r1x = r10 + gridX * ( r11 - r10); + vector[dim] = gridY * (r1x - r0x) + r0x; + if (derivative) { + /* + * Following code appends the same values than the ones computed by derivativeInCell(gridX, gridY), + * but reusing some of the values that we already fetched for computing the interpolation. + */ + double dy = r10 - r00; + int i = n; + if (dim == 0) { + dx++; + } else { + dy++; + i += INTERPOLATED_DIMENSIONS; + derivative = false; + } + vector[i ] = dx; + vector[i+1] = dy; + } } } @@ -534,9 +576,11 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T final int iy = Math.max(0, Math.min(gridSize[1] - 2, (int) gridY)); final Matrix derivative = Matrices.createDiagonal(getTranslationDimensions(), gridSize.length); for (int j=derivative.getNumRow(); --j>=0;) { - final double orig = getCellValue(j, ix, iy); - derivative.setElement(j, 0, derivative.getElement(j, 0) + (getCellValue(j, ix+1, iy) - orig)); - derivative.setElement(j, 1, derivative.getElement(j, 1) + (getCellValue(j, ix, iy+1) - orig)); + final double r00 = getCellValue(j, ix, iy ); + final double dx = getCellValue(j, ix+1, iy ) - r00; + final double dy = getCellValue(j, ix, iy+1) - r00; + derivative.setElement(j, 0, derivative.getElement(j, 0) + dx); + derivative.setElement(j, 1, derivative.getElement(j, 1) + dy); } return derivative; } diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedTransform.java index 3738512..8fc5bcc 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedTransform.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/InterpolatedTransform.java @@ -70,7 +70,7 @@ import org.apache.sis.internal.referencing.provider.DatumShiftGridFile; * @author Martin Desruisseaux (IRD, Geomatys) * @author Simon Reynard (Geomatys) * @author Rueben Schulz (UBC) - * @version 0.8 + * @version 1.0 * * @see DatumShiftGrid * @see org.apache.sis.referencing.operation.builder.LocalizationGridBuilder @@ -432,6 +432,20 @@ public class InterpolatedTransform extends DatumShiftTransform { private static final long serialVersionUID = 4335801994727826360L; /** + * Set to {@code true} for debugging. + */ + private static final boolean DEBUG = false; + + /** + * Whether to use a simple version of this inverse transform (where the Jacobian matrix is close to identity + * in every points) or a more complex version using derivatives (Jacobian matrices) for adjusting the errors. + * The simple mode is okay for transformations based on e.g. NADCON or NTv2 grids, but is not sufficient for + * more curved grids like the ones read from netCDF files. We provide this switch for allowing performance + * comparisons. + */ + private static final boolean SIMPLE = true; + + /** * Maximum number of iterations. This is set to a higher value than {@link Formulas#MAXIMUM_ITERATIONS} because * the data used in {@link InterpolatedTransform} grid may come from anywhere, in particular localization grids * in netCDF files. Deformations in those grids may be much higher than e.g. {@link DatumShiftTransform} grids. @@ -490,51 +504,12 @@ public class InterpolatedTransform extends DatumShiftTransform { dstPts = new double[dimension]; dstOff = 0; } - double xi, yi; - final double x = xi = srcPts[srcOff ]; - final double y = yi = srcPts[srcOff+1]; - final double[] vector = new double[dimension]; - int it = MAXIMUM_ITERATIONS; - double tol = tolerance; - do { - /* - * We want (xi, yi) such as the following conditions hold: - * - * xi + vector[0] ≈ x ⟶ xi ≈ x - vector[0] - * yi + vector[1] ≈ y ⟶ yi ≈ y - vector[1] - */ - forward.grid.interpolateInCell(xi, yi, vector); - final double ox = xi; - final double oy = yi; - xi = x - vector[0]; - yi = y - vector[1]; - if (!(Math.abs(xi - ox) > tol || // Use '!' for catching NaN. - Math.abs(yi - oy) > tol)) - { - if (dimension > GRID_DIMENSION) { - System.arraycopy(srcPts, srcOff + GRID_DIMENSION, - dstPts, dstOff + GRID_DIMENSION, - dimension - GRID_DIMENSION); - /* - * We can not use srcPts[srcOff + i] = dstPts[dstOff + i] + offset[i] - * because the arrays may overlap. The contract said that this method - * must behave as if all input coordinate values have been read before - * we write outputs, which is the reason for System.arraycopy(…) call. - */ - int i = dimension; - do dstPts[dstOff + --i] -= vector[i]; - while (i > GRID_DIMENSION); - } - dstPts[dstOff ] = xi; // Shall not be done before above loop. - dstPts[dstOff+1] = yi; - if (derivate) { - return Matrices.inverse(forward.derivative( - new DirectPositionView.Double(dstPts, dstOff, dimension))); - } - return null; - } - } while (--it >= 0 || (tol = tryAgain(it, xi, yi)) > 0); - throw new TransformException(Resources.format(Resources.Keys.NoConvergence)); + transform(srcPts, srcOff, dstPts, dstOff, 1); + if (derivate) { + return Matrices.inverse(forward.derivative(new DirectPositionView.Double(dstPts, dstOff, dimension))); + } else { + return null; + } } /** @@ -543,6 +518,7 @@ public class InterpolatedTransform extends DatumShiftTransform { * @throws TransformException if a point can not be transformed. */ @Override + @SuppressWarnings("UseOfSystemOutOrSystemErr") // Used if DEBUG = true only. public final void transform(double[] srcPts, int srcOff, final double[] dstPts, int dstOff, int numPts) throws TransformException { @@ -566,22 +542,76 @@ public class InterpolatedTransform extends DatumShiftTransform { } } } - final double[] vector = new double[dimension]; + final double[] vector = new double[SIMPLE ? dimension : dimension + 4]; // +4 is for derivate coefficients. nextPoint: while (--numPts >= 0) { - double xi, yi; + double xi, yi; // (x,y) values at iteration i. final double x = xi = srcPts[srcOff ]; final double y = yi = srcPts[srcOff+1]; + if (DEBUG) { + System.out.format("Searching source coordinates for target = (%g %g)%n", x, y); + } int it = MAXIMUM_ITERATIONS; double tol = tolerance; do { forward.grid.interpolateInCell(xi, yi, vector); - final double ox = xi; - final double oy = yi; - xi = x - vector[0]; - yi = y - vector[1]; - if (!(Math.abs(xi - ox) > tol || // Use '!' for catching NaN. - Math.abs(yi - oy) > tol)) - { + final boolean more; + if (SIMPLE) { + /* + * We want (xi, yi) such as the following conditions hold + * (see next commnt for the simplification applied here): + * + * xi + vector[0] ≈ x ⟶ xi ≈ x - vector[0] + * yi + vector[1] ≈ y ⟶ yi ≈ y - vector[1] + */ + final double ox = xi; + final double oy = yi; + xi = x - vector[0]; + yi = y - vector[1]; + more = Math.abs(xi - ox) > tol || Math.abs(yi - oy) > tol; + } else { + /* + * The error between the new position (xi + tx) and the desired position x is measured + * in target units. In order to apply a correction in opposite direction, we need to + * convert the translation vector from target units to source units. The above simple + * case was assuming that this conversion is close to identity transform, allowing us + * to write xi = xi - ex with the error ex = (xi + tx) - x, which simplified as + * xi = x - tx. But if the grid is more curved, we can not assume anymore that error + * ex in target units is approximately equal to error dx in source units. A conversion + * needs to be applied, by multiplying the translation error by the derivative of the + * "target to source" transform. Since we do not know that derivative, we use inverse + * of the "source to target" transform derivative. + */ + final double tx = vector[0]; + final double ty = vector[1]; + final double m00 = vector[dimension ]; + final double m01 = vector[dimension + 1]; + final double m10 = vector[dimension + 2]; + final double m11 = vector[dimension + 3]; + final double det = m00 * m11 - m01 * m10; + final double ex = (xi + tx) - x; // Errors in target units. + final double ey = (yi + ty) - y; + final double dx = (ex * m11 - ey * m01) / det; // Errors in source units. + final double dy = (ey * m00 - ex * m10) / det; + if (DEBUG) { + Matrix m = forward.grid.derivativeInCell(xi, yi); + assert m.getElement(0,0) == m00; + assert m.getElement(0,1) == m01; + assert m.getElement(1,0) == m10; + assert m.getElement(1,1) == m11; + + m = Matrices.inverse(m); + double[] c = ((MatrixSIS) m).multiply(new double[] {ex, ey}); + assert Math.abs(dx - c[0]) < 1E-5 : Arrays.toString(c) + " : " + dx; + assert Math.abs(dy - c[1]) < 1E-5 : Arrays.toString(c) + " : " + dy; + + System.out.printf(" source=(%8.2f %8.2f) target=(%8.2f %8.2f) error=(%8.2f %8.2f) → (%7.2f %7.2f)%n", + xi, yi, (xi+tx), (yi+ty), ex, ey, dx, dy); + } + xi -= dx; + yi -= dy; + more = Math.abs(ex) > tol || Math.abs(ey) > tol; + } + if (!more) { // Use '!' for catching NaN. if (dimension > GRID_DIMENSION) { System.arraycopy(srcPts, srcOff + GRID_DIMENSION, dstPts, dstOff + GRID_DIMENSION, diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressedTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressedTest.java new file mode 100644 index 0000000..2b434b8 --- /dev/null +++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/DatumShiftGridCompressedTest.java @@ -0,0 +1,46 @@ +/* + * 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.internal.referencing.provider; + +import javax.measure.quantity.Dimensionless; +import org.opengis.referencing.operation.NoninvertibleTransformException; + +import static org.opengis.test.Assert.*; + + +/** + * Tests {@link DatumShiftGridCompressed}. This class creates a grid using values computed by an affine transform, + * and compare values computed by the grid using the affine transform as a reference. + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.0 + * @since 1.0 + * @module + */ +public final strictfp class DatumShiftGridCompressedTest extends DatumShiftGridFileTest { + /** + * Creates a new grid using an affine transform as a reference. + * + * @param rotation ignored. + */ + @Override + void init(final double rotation) throws NoninvertibleTransformException { + super.init(0); // No rotation in order to have integer values. + grid = DatumShiftGridCompressed.compress((DatumShiftGridFile.Float<Dimensionless,Dimensionless>) grid, null, 0.5); + assertInstanceOf("grid", DatumShiftGridCompressed.class, grid); + } +} diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFileTest.java b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFileTest.java new file mode 100644 index 0000000..cad1710 --- /dev/null +++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFileTest.java @@ -0,0 +1,194 @@ +/* + * 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.internal.referencing.provider; + +import java.util.Random; +import java.awt.geom.Point2D; +import java.awt.geom.AffineTransform; +import javax.measure.quantity.Dimensionless; +import org.opengis.referencing.operation.TransformException; +import org.opengis.referencing.operation.NoninvertibleTransformException; +import org.apache.sis.measure.Units; +import org.apache.sis.test.TestCase; +import org.apache.sis.test.TestUtilities; +import org.apache.sis.test.DependsOnMethod; +import org.junit.Test; + +import static org.junit.Assert.*; + + +/** + * Tests {@link DatumShiftGridFile}. This class creates a grid using values computed by an affine transform, + * and compare values computed by the grid using the affine transform as a reference. + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.0 + * @since 1.0 + * @module + */ +public strictfp class DatumShiftGridFileTest extends TestCase { + /** + * Size of the grid created for testing purpose. + */ + private static final int WIDTH = 30, HEIGHT = 45; + + /** + * Tolerance threshold used in this test. We use a relatively high threshold because translations are stored + * in {@code float[]} arrays, so the tolerance needs to be slightly higher than {@code float} precision. + */ + private static final float TOLERANCE = 2E-5f; + + /** + * Number of random points to test in each method. + */ + private static final int NUM_TESTS = 20; + + /** + * The transform to use as a reference. + */ + private AffineTransform reference; + + /** + * The grid wrapping computed from affine transform. + */ + DatumShiftGridFile<Dimensionless,Dimensionless> grid; + + /** + * Creates a new grid using an affine transform as a reference. + * An arbitrary non-uniform scale is applied on axes. + * + * @param rotation rotation angle in degrees to apply on affine transform. + */ + void init(final double rotation) throws NoninvertibleTransformException { + reference = AffineTransform.getRotateInstance(StrictMath.toRadians(rotation), WIDTH/2, HEIGHT/2); + reference.scale(2, 5); + final DatumShiftGridFile.Float<Dimensionless,Dimensionless> grid = new DatumShiftGridFile.Float<>( + 2, Units.UNITY, Units.UNITY, true, 0, 0, 1, 1, WIDTH, HEIGHT, null); + assertEquals(2, grid.offsets.length); + final Point2D.Float point = new Point2D.Float(); + int i = 0; + for (int y=0; y<HEIGHT; y++) { + for (int x=0; x<WIDTH; x++) { + point.x = x; + point.y = y; + assertSame(point, reference.transform(point, point)); + grid.offsets[0][i] = point.x - x; + grid.offsets[1][i] = point.y - y; + i++; + } + } + assertEquals(grid.offsets[0].length, i); + assertEquals(grid.offsets[1].length, i); + this.grid = grid; + } + + /** + * Tests {@link DatumShiftGridFile#interpolateInCell(double, double, double[])}. + * This test does not perform interpolations and does not compute derivatives. + * + * @throws TransformException if an error occurred while transforming a coordinates. + */ + @Test + public void testInterpolateAtIntegers() throws TransformException { + final Random random = TestUtilities.createRandomNumberGenerator(); + final Point2D.Float point = new Point2D.Float(); + final double[] vector = new double[2]; + init(15); + for (int i=0; i<NUM_TESTS; i++) { + /* + * Compute the reference point. + */ + final int x = random.nextInt(WIDTH); + final int y = random.nextInt(HEIGHT); + point.x = x; + point.y = y; + assertSame(point, reference.transform(point, point)); + /* + * Compute the actual point and compare. + */ + grid.interpolateInCell(x, y, vector); + assertEquals("x", point.x, (float) (vector[0] + x), TOLERANCE); + assertEquals("y", point.y, (float) (vector[1] + y), TOLERANCE); + } + } + + /** + * Tests {@link DatumShiftGridFile#interpolateAt(double...)}. + * This tests include interpolations. + * + * @throws TransformException if an error occurred while transforming a coordinates. + */ + @Test + @DependsOnMethod("testInterpolateAtIntegers") + public void testInterpolateAtReals() throws TransformException { + final Random random = TestUtilities.createRandomNumberGenerator(); + final Point2D.Float point = new Point2D.Float(); + final double[] vector = new double[2]; + init(0); // No rotation for having same interpolations. + for (int i=0; i<NUM_TESTS; i++) { + /* + * Compute the reference point. + */ + final float x = random.nextFloat() * (WIDTH - 1); + final float y = random.nextFloat() * (HEIGHT - 1); + point.x = x; + point.y = y; + assertSame(point, reference.transform(point, point)); + /* + * Compute the actual point and compare. + */ + grid.interpolateInCell(x, y, vector); + assertEquals("x", point.x, (float) (vector[0] + x), TOLERANCE); + assertEquals("y", point.y, (float) (vector[1] + y), TOLERANCE); + } + } + + /** + * Tests {@link DatumShiftGridFile#interpolateAt(double...)} with opportunistic derivative calculations. + * Since the grid is computed from an affine transform, the derivative should be constant everywhere. + * + * @throws TransformException if an error occurred while transforming a coordinates. + */ + @Test + @DependsOnMethod("testInterpolateAtIntegers") + public void testInterpolateAndDerivative() throws TransformException { + final Random random = TestUtilities.createRandomNumberGenerator(); + final Point2D.Float point = new Point2D.Float(); + final double[] vector = new double[6]; + init(30); + for (int i=0; i<NUM_TESTS; i++) { + /* + * Compute the reference point. + */ + final int x = random.nextInt(WIDTH); + final int y = random.nextInt(HEIGHT); + point.x = x; + point.y = y; + assertSame(point, reference.transform(point, point)); + /* + * Compute the actual point, compare, then check derivative. + */ + grid.interpolateInCell(x, y, vector); + assertEquals("x", point.x, (float) (vector[0] + x), TOLERANCE); + assertEquals("y", point.y, (float) (vector[1] + y), TOLERANCE); + assertEquals("m00", reference.getScaleX(), vector[2], TOLERANCE); + assertEquals("m01", reference.getShearX(), vector[3], TOLERANCE); + assertEquals("m10", reference.getShearY(), vector[4], TOLERANCE); + assertEquals("m11", reference.getScaleY(), vector[5], TOLERANCE); + } + } +} diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/InterpolatedTransformTest.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/InterpolatedTransformTest.java index 1e273b4..a983095 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/InterpolatedTransformTest.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/InterpolatedTransformTest.java @@ -17,6 +17,7 @@ package org.apache.sis.referencing.operation.transform; import java.net.URL; +import java.util.Arrays; import org.opengis.util.FactoryException; import org.opengis.parameter.ParameterValueGroup; import org.opengis.referencing.operation.MathTransformFactory; @@ -35,12 +36,18 @@ import org.apache.sis.test.DependsOn; import org.junit.Test; - /** * Tests {@link InterpolatedTransform}. + * Tested transformations are: + * + * <ul> + * <li>Simple case based on linear calculations (easier to debug).</li> + * <li>From NTF to RGF93 using a NTv2 grid.</li> + * <li>From NAD27 to NAD83 using a NADCON grid.</li> + * </ul> * * @author Martin Desruisseaux (Geomatys) - * @version 0.7 + * @version 1.0 * @since 0.7 * @module */ @@ -50,6 +57,21 @@ import org.junit.Test; }) public final strictfp class InterpolatedTransformTest extends MathTransformTestCase { /** + * Creates an {@link InterpolatedTransform} derived from a quadratic formula. + * We do not really need {@code InterpolatedTransform} for quadratic formulas, + * but we use them for testing purpose since they are easier to debug. + * + * @param rotation rotation angle, in degrees. Use 0 for debugging a simple case. + * @return suggested points to use for testing purposes as an array of length 2, + * with source coordinates in the first array and target coordinates in the second array. + */ + private double[][] createQuadratic(final double rotation) throws TransformException { + final QuadraticShiftGrid grid = new QuadraticShiftGrid(rotation); + transform = new InterpolatedTransform(grid); + return grid.samplePoints(); + } + + /** * Creates the same transformation than <cite>"France geocentric interpolation"</cite> transform * (approximately), but using shifts in geographic domain instead than in geocentric domain. * @@ -66,7 +88,7 @@ public final strictfp class InterpolatedTransformTest extends MathTransformTestC } /** - * Creates a transformation from NAD27 to NAD93 + * Creates a transformation from NAD27 to NAD93. * * @throws FactoryException if an error occurred while loading the grid. */ @@ -83,90 +105,131 @@ public final strictfp class InterpolatedTransformTest extends MathTransformTestC } /** - * Tests forward transformation of sample points. Tested transformations are: - * <ul> - * <li>From NTF to RGF93 using a NTv2 grid.</li> - * <li>From NAD27 to NAD83 using a NADCON grid.</li> - * </ul> + * Tests forward transformation of sample points. * - * @throws FactoryException if an error occurred while loading a grid. * @throws TransformException if an error occurred while transforming a coordinate. - * - * @see InterpolatedGeocentricTransformTest#testInverseTransform() */ @Test - public void testForwardTransform() throws FactoryException, TransformException { - isInverseTransformSupported = false; - createRGF93(); - verifyTransform(FranceGeocentricInterpolationTest.samplePoint(1), - FranceGeocentricInterpolationTest.samplePoint(3)); - createNADCON(); - verifyTransform(NADCONTest.samplePoint(1), - NADCONTest.samplePoint(3)); + public void testForwardTransform() throws TransformException { + isInverseTransformSupported = false; // For focusing on a single aspect. + final double[][] samplePoints = createQuadratic(-15); + tolerance = 1E-10; + verifyTransform(Arrays.copyOf(samplePoints[0], QuadraticShiftGrid.FIRST_FRACTIONAL_COORDINATE), + Arrays.copyOf(samplePoints[1], QuadraticShiftGrid.FIRST_FRACTIONAL_COORDINATE)); + + tolerance = 0.003; // Because of interpolations in fractional coordinates. + verifyTransform(samplePoints[0], samplePoints[1]); } /** - * Tests inverse transformation of sample points. Tested transformations are: - * <ul> - * <li>From RGF93 to NTF using a NTv2 grid.</li> - * <li>From NAD83 to NAD27 using a NADCON grid.</li> - * </ul> + * Tests inverse transformation of sample points. + * Inverse transform requires derivative. * - * @throws FactoryException if an error occurred while loading a grid. * @throws TransformException if an error occurred while transforming a coordinate. */ @Test - @DependsOnMethod("testForwardTransform") - public void testInverseTransform() throws FactoryException, TransformException { - isInverseTransformSupported = false; - createRGF93(); - transform = transform.inverse(); - verifyTransform(FranceGeocentricInterpolationTest.samplePoint(3), - FranceGeocentricInterpolationTest.samplePoint(1)); - createNADCON(); + @org.junit.Ignore("Debugging still in progress") + @DependsOnMethod({"testForwardTransform", "testDerivative"}) + public void testInverseTransform() throws TransformException { + isInverseTransformSupported = false; // For focusing on a single aspect. + final double[][] samplePoints = createQuadratic(-20); transform = transform.inverse(); - verifyTransform(NADCONTest.samplePoint(3), - NADCONTest.samplePoint(1)); + tolerance = QuadraticShiftGrid.PRECISION; + verifyTransform(samplePoints[1], samplePoints[0]); } /** * Tests the derivatives at the sample point. This method compares the derivatives computed by * the transform with an estimation of derivatives computed by the finite differences method. * - * @throws FactoryException if an error occurred while loading the grid. * @throws TransformException if an error occurred while transforming the coordinate. */ @Test @DependsOnMethod("testForwardTransform") - public void testForwardDerivative() throws FactoryException, TransformException { - createRGF93(); - final double delta = 0.2; - derivativeDeltas = new double[] {delta, delta}; - tolerance = 5E-6; // Empirical value. - verifyDerivative(FranceGeocentricInterpolationTest.samplePoint(1)); + public void testDerivative() throws TransformException { + final double[][] samplePoints = createQuadratic(-40); + final double[] point = new double[QuadraticShiftGrid.DIMENSION]; // A single point from 'samplePoints' + derivativeDeltas = new double[] {0.002, 0.002}; + isInverseTransformSupported = false; // For focusing on a single aspect. + for (int i=0; i < samplePoints[0].length; i += QuadraticShiftGrid.DIMENSION) { + System.arraycopy(samplePoints[0], i, point, 0, QuadraticShiftGrid.DIMENSION); + if (i < QuadraticShiftGrid.FIRST_FRACTIONAL_COORDINATE) { + tolerance = 1E-10; // Empirical value. + } else { + tolerance = 0.003; // Because current implementation does not yet interpolate derivatives. + } + verifyDerivative(point); + /* + * Verify derivative at the same point but using inverse transform, + * done in same loop for easier to comparisons during debugging. + */ + if (isInverseTransformSupported) { + transform = transform.inverse(); + System.arraycopy(samplePoints[1], i, point, 0, QuadraticShiftGrid.DIMENSION); + verifyDerivative(point); + transform = transform.inverse(); // Back to forward transform. + } + } } /** - * Tests the derivatives at the sample point. This method compares the derivatives computed by - * the transform with an estimation of derivatives computed by the finite differences method. + * Performs the tests using the same transformation than <cite>"France geocentric interpolation"</cite> + * transform (approximately), but using shifts in geographic domain instead than in geocentric domain. * - * @throws FactoryException if an error occurred while loading the grid. - * @throws TransformException if an error occurred while transforming the coordinate. + * @throws FactoryException if an error occurred while creating a transform. + * @throws TransformException if an error occurred while transforming a coordinate. + * + * @see InterpolatedGeocentricTransformTest#testInverseTransform() */ @Test - @DependsOnMethod("testInverseTransform") - public void testInverseDerivative() throws FactoryException, TransformException { + public void testRGF93() throws FactoryException, TransformException { createRGF93(); + + // Forward transform + isInverseTransformSupported = true; // Set to 'false' for testing one direction at time. + verifyTransform(FranceGeocentricInterpolationTest.samplePoint(1), + FranceGeocentricInterpolationTest.samplePoint(3)); + + // Inverse transform + transform = transform.inverse(); + verifyTransform(FranceGeocentricInterpolationTest.samplePoint(3), + FranceGeocentricInterpolationTest.samplePoint(1)); + + // Forward derivative + transform = transform.inverse(); + derivativeDeltas = new double[] {0.2, 0.2}; + tolerance = 5E-6; // Empirical value. + verifyDerivative(FranceGeocentricInterpolationTest.samplePoint(1)); + + // Inverse derivative transform = transform.inverse(); - final double delta = 0.2; - derivativeDeltas = new double[] {delta, delta}; - tolerance = 5E-6; // Empirical value. verifyDerivative(FranceGeocentricInterpolationTest.samplePoint(3)); } /** + * Performs the tests using the transformation from NAD27 to NAD93. + * + * @throws FactoryException if an error occurred while creating a transform. + * @throws TransformException if an error occurred while transforming a coordinate. + */ + @Test + public void testNADCON() throws FactoryException, TransformException { + createNADCON(); + + // Forward transform + isInverseTransformSupported = true; // Set to 'false' for testing one direction at time. + verifyTransform(NADCONTest.samplePoint(1), + NADCONTest.samplePoint(3)); + + // Inverse transform + transform = transform.inverse(); + verifyTransform(NADCONTest.samplePoint(3), + NADCONTest.samplePoint(1)); + } + + /** * Tests the Well Known Text (version 1) formatting. - * The result is what we show to users, but may quite different than what SIS has in memory. + * The result is what we show to users, but may be quite different than what SIS has in memory. * * @throws FactoryException if an error occurred while creating a transform. * @throws TransformException should never happen. diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformTestCase.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformTestCase.java index dc2be7b..00a7c59 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformTestCase.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/MathTransformTestCase.java @@ -77,7 +77,7 @@ import org.opengis.test.CalculationType; */ public abstract strictfp class MathTransformTestCase extends TransformTestCase { /** - * The number of ordinates to use for stressing the math transform. We use a number that + * The number of coordinates to use for stressing the math transform. We use a number that * encompass at least 2 time the default buffer size in order to test the code that use * the buffer. We add an arbitrary number just for making the transform job harder. */ @@ -85,7 +85,7 @@ public abstract strictfp class MathTransformTestCase extends TransformTestCase { /** * The dimension of longitude, or {@code null} if none. If non-null, then the comparison of - * ordinate values along that dimension will ignore 360° offsets. + * coordinate values along that dimension will ignore 360° offsets. * * <p>The first array element is the dimension during forward transforms, and the second * array element is the dimension during inverse transforms (can be omitted if the later @@ -105,7 +105,7 @@ public abstract strictfp class MathTransformTestCase extends TransformTestCase { /** * The tolerance level for height above the ellipsoid. This tolerance is usually higher - * than the {@linkplain #tolerance tolerance} level for horizontal ordinate values. + * than the {@linkplain #tolerance tolerance} level for horizontal coordinate values. */ protected double zTolerance; @@ -158,8 +158,8 @@ public abstract strictfp class MathTransformTestCase extends TransformTestCase { * The SIS implementation ensures that longitude values are contained in the ±180° range, * applying 360° shifts if needed. * - * @param expected the expected ordinate value provided by the test case. - * @param actual the ordinate value computed by the {@linkplain #transform transform} being tested. + * @param expected the expected coordinate values provided by the test case. + * @param actual the coordinate values computed by the {@linkplain #transform transform} being tested. * @param mode indicates if the coordinates being compared are the result of a direct * or inverse transform, or if strict equality is requested. */ @@ -219,18 +219,17 @@ public abstract strictfp class MathTransformTestCase extends TransformTestCase { /** * Transforms the given coordinates and verifies that the result is equals (within a positive delta) - * to the expected ones. If the difference between an expected and actual ordinate value is greater + * to the expected ones. If the difference between an expected and actual coordinate value is greater * than the {@linkplain #tolerance tolerance} threshold, then the assertion fails. * * <p>If {@link #isInverseTransformSupported} is {@code true}, then this method will also transform - * the expected coordinate points using the {@linkplain MathTransform#inverse() inverse transform} and - * compare with the source coordinates.</p> + * the expected points using the {@linkplain MathTransform#inverse() inverse transform} and compare + * with the source coordinates.</p> * * <p>This method verifies also the consistency of {@code MathTransform.transform(…)} method variants.</p> * - * @param coordinates the coordinate points to transform. - * @param expected the expect result of the transformation, or - * {@code null} if {@code coordinates} is expected to be null. + * @param coordinates the points to transform. + * @param expected the expected transformation results, or {@code null} if {@code coordinates} is expected to be null. * @throws TransformException if the transformation failed. */ @Override @@ -250,10 +249,10 @@ public abstract strictfp class MathTransformTestCase extends TransformTestCase { } /* * The comparison below needs a higher tolerance threshold, because we converted the source - * ordinates to floating points which induce a lost of precision. The multiplication factor + * coordinates to floating points which induce a lost of precision. The multiplication factor * used here has been determined empirically. The value is quite high, but this is only an * oportunist check anyway. The "real" test is the one performed by 'verifyConsistency'. - * We do not perform this check for non-linear transforms, because the difference in input + * We do not perform this check for non-linear transforms, because the differences in input * have too unpredictable consequences on the output. */ if (transform instanceof LinearTransform) { @@ -269,7 +268,7 @@ public abstract strictfp class MathTransformTestCase extends TransformTestCase { } /** - * Stress the current {@linkplain #transform transform} using random ordinates in the given domain. + * Stress the current {@linkplain #transform transform} using random coordinates in the given domain. * First, this method creates a grid of regularly spaced points along all dimensions in the given domain. * Next, this method adds small random displacements to every points and shuffle the coordinates in random order. * Finally this method delegates the resulting array of coordinates to the following methods: diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/QuadraticShiftGrid.java b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/QuadraticShiftGrid.java new file mode 100644 index 0000000..b7ff5d4 --- /dev/null +++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/QuadraticShiftGrid.java @@ -0,0 +1,150 @@ +/* + * 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.awt.geom.AffineTransform; +import javax.measure.quantity.Dimensionless; +import org.apache.sis.referencing.datum.DatumShiftGrid; +import org.apache.sis.measure.Units; + + +/** + * Dummy implementation of {@link DatumShiftGrid} containing translation vectors that are computed by a quadratic function. + * This class has no computational interest compared to a direct implementation of quadratic formulas, but is convenient + * for debugging {@link InterpolatedTransform} because of its predictable behavior (easier to see with rotation of 0°). + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.0 + * @since 1.0 + * @module + */ +@SuppressWarnings("serial") // Not intended to be serialized. +final strictfp class QuadraticShiftGrid extends DatumShiftGrid<Dimensionless,Dimensionless> { + /** + * Number of source and target dimensions of the grid. + */ + static final int DIMENSION = 2; + + /** + * The precision used for inverse transforms, in units of cells. + * A smaller numbers will cause more iterations to occur. + */ + static final double PRECISION = 1E-7; + + /** + * Index of the first point in {@link #samplePoints()} which contains a fractional part. + * Those points need a different tolerance threshold because the interpolations performed + * by {@link DatumShiftGrid} is not the same than the formula used in this test class. + * + * @see #samplePoints() + */ + static final int FIRST_FRACTIONAL_COORDINATE = 3 * DIMENSION; + + /** + * Grid size in number of pixels. The translation vectors in the middle of this grid will be (0,0). + */ + private static final int WIDTH = 200, HEIGHT = 2000; + + /** + * An internal step performed during computation of translation vectors at a given point. + */ + private final AffineTransform offsets; + + /** + * A temporary buffer for calculations with {@link #offsets}. + */ + private final double[] buffer; + + /** + * Creates a new grid mock of size {@value #WIDTH} × {@value #HEIGHT} pixels. + * The translation vectors will be (0,0) in the middle of that grid and increase + * as we get further from the center. + * + * @param rotation the rotation angle, in degrees. + */ + QuadraticShiftGrid(final double rotation) { + super(Units.UNITY, MathTransforms.identity(DIMENSION), new int[] {WIDTH, HEIGHT}, true, Units.UNITY); + offsets = AffineTransform.getRotateInstance(StrictMath.toRadians(rotation)); + offsets.scale(-0.1, 0.025); + offsets.translate(-0.5*WIDTH, -0.5*HEIGHT); + buffer = new double[DIMENSION]; + } + + /** + * Returns the number of dimensions of the translation vectors interpolated by this shift grid. + */ + @Override + public int getTranslationDimensions() { + return DIMENSION; + } + + /** + * Returns suggested points to use for testing purposes. This method returns an array of length 2, + * with source coordinates in the first array and target coordinates in the second array. + * + * @see #FIRST_FRACTIONAL_COORDINATE + */ + final double[][] samplePoints() { + final double[] sources = { + /*[0]*/ WIDTH/2, HEIGHT/2, + /*[1]*/ 50, 1400, + /*[2]*/ 175, 200, + /*[3]*/ 10.356, 30.642 // FIRST_FRACTIONAL_COORDINATE must point here. + }; + final double[] targets = new double[sources.length]; + transform(sources, targets); + return new double[][] {sources, targets}; + } + + /** + * Applies an arbitrary non-linear transform on the given source coordinates + * and stores the results in the given target array. + */ + private void transform(final double[] sources, final double[] targets) { + offsets.transform(sources, 0, targets, 0, sources.length / DIMENSION); + for (int i=0; i<targets.length; i++) { + final double t = targets[i]; + targets[i] = StrictMath.copySign(t*t, t); + } + for (int i=0; i<targets.length;) { + targets[i++] += WIDTH / 2; + targets[i++] += HEIGHT / 2; + } + } + + /** + * Returns the cell value at the given dimension and grid index. + * Those values are components of <em>translation</em> vectors. + */ + @Override + public double getCellValue(int dim, int gridX, int gridY) { + buffer[0] = gridX; + buffer[1] = gridY; + transform(buffer, buffer); + buffer[0] -= gridX; + buffer[1] -= gridY; + return buffer[dim]; + } + + /** + * Returns an arbitrary cell precision. This determines when the iteration algorithm stops. + */ + @Override + public double getCellPrecision() { + return PRECISION; + } +} diff --git a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java index 5494363..040e08d 100644 --- a/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java +++ b/core/sis-referencing/src/test/java/org/apache/sis/test/suite/ReferencingTestSuite.java @@ -148,6 +148,8 @@ import org.junit.BeforeClass; org.apache.sis.internal.referencing.provider.PositionVector7ParamTest.class, org.apache.sis.internal.referencing.provider.CoordinateFrameRotationTest.class, org.apache.sis.internal.referencing.provider.MolodenskyTest.class, + org.apache.sis.internal.referencing.provider.DatumShiftGridFileTest.class, + org.apache.sis.internal.referencing.provider.DatumShiftGridCompressedTest.class, org.apache.sis.internal.referencing.provider.FranceGeocentricInterpolationTest.class, org.apache.sis.internal.referencing.provider.NTv2Test.class, org.apache.sis.internal.referencing.provider.NADCONTest.class,
