sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 03/03: First draft of an InterpolatedTransform algorithm using derivative (Jacobian matrix) for deciding in which direction to compensate for errors. Not yet enabled.
Date Fri, 12 Apr 2019 18:22:46 GMT
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 <martin.desruisseaux@geomatys.com>
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,


Mime
View raw message