sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch geoapi-4.0 updated: Calculation of derivative outside the grid must be consistent with the way values are extrapolated by `interpolateInCell`.
Date Mon, 15 Apr 2019 10:55:10 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


The following commit(s) were added to refs/heads/geoapi-4.0 by this push:
     new 1f24dff  Calculation of derivative outside the grid must be consistent with the way
values are extrapolated by `interpolateInCell`.
1f24dff is described below

commit 1f24dff960024706b4c7041e970de3a982548c4c
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Apr 15 12:53:55 2019 +0200

    Calculation of derivative outside the grid must be consistent with the way values are
extrapolated by `interpolateInCell`.
---
 .../provider/DatumShiftGridCompressed.java         | 43 +++++++----
 .../sis/referencing/datum/DatumShiftGrid.java      | 90 ++++++++++++++++------
 .../provider/DatumShiftGridFileTest.java           | 52 ++++++++++++-
 .../transform/InterpolatedTransformTest.java       | 33 +++++++-
 .../operation/transform/SinusoidalShiftGrid.java   |  5 ++
 5 files changed, 178 insertions(+), 45 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 4eee2d8..3527244 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
@@ -184,35 +184,50 @@ final class DatumShiftGridCompressed<C extends Quantity<C>,
T extends Quantity<T
      */
     @Override
     public void interpolateInCell(double gridX, double gridY, double[] vector) {
-        if (gridX < 0) gridX = 0;
-        if (gridY < 0) gridY = 0;
+        boolean skipX = false;
+        boolean skipY = false;                          // Whether to skip derivative calculation
for X or Y.
+        if (gridX < 0) {gridX = 0; skipX = true;}
+        if (gridY < 0) {gridY = 0; skipY = true;}
         int ix = (int) gridX;  gridX -= ix;
         int iy = (int) gridY;  gridY -= iy;
         if (ix > nx - 2) {
-            ix = nx - 2;
-            gridX = 1;
+            skipX |= (ix != nx-1 || gridX != 0);        // Keep value 'false' if gridX ==
gridSize[0] - 1.
+            ix     = nx - 2;
+            gridX  = 1;
         }
-        if (iy > ymax) {             // Subtraction of 2 already done by the constructor.
-            iy = ymax;
-            gridY = 1;
+        if (iy > ymax) {                                // Subtraction of 2 already done
by the constructor.
+            skipY |= (iy != ymax+1 || gridY != 0);      // Keep value 'false' if gridY ==
gridSize[1] - 1.
+            iy     = ymax;
+            gridY  = 1;
         }
         final int p00 = nx*iy + ix;
         final int p10 = nx + p00;
         final int n   = data.length;
-        boolean derivative = (vector.length >= n + INTERPOLATED_DIMENSIONS*INTERPOLATED_DIMENSIONS);
+        boolean derivative = (vector.length >= n + INTERPOLATED_DIMENSIONS * INTERPOLATED_DIMENSIONS);
         for (int dim = 0; dim < n; dim++) {
-            double dx;
+            double dx, dy;
             final short[] values = data[dim];
             final double r00 = values[p00    ];
-            final double r01 = values[p00 + 1];         // Naming convention: ryx (row index
first, like matrix).
+            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);
+            final double r0x = r00 + gridX * (dx = r01 - r00);      // TODO: use Math.fma
on JDK9.
+            final double r1x = r10 + gridX * (dy = r11 - r10);      // Not really "dy" measurement
yet, will become dy later.
             vector[dim] = (gridY * (r1x - r0x) + r0x) * scale + averages[dim];
             if (derivative) {
-                double dy = (r10 - r00) * scale;
-                dx *= scale;
+                if (skipX) {
+                    dx = 0;
+                } else {
+                    dx += (dy - dx) * gridX;
+                    dx *= scale;
+                }
+                if (skipY) {
+                    dy = 0;
+                } else {
+                    dy  =  r10 - r00;
+                    dy += (r11 - r01 - dy) * gridY;
+                    dy *= scale;
+                }
                 int i = n;
                 if (dim == 0) {
                     dx++;
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 a34d548..e9efe3c 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
@@ -445,36 +445,50 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
      * @see #isCellInGrid(double, double)
      */
     public void interpolateInCell(double gridX, double gridY, final double[] vector) {
-        if (gridX < 0) gridX = 0;
-        if (gridY < 0) gridY = 0;
+        boolean skipX = false;
+        boolean skipY = false;                          // Whether to skip derivative calculation
for X or Y.
+        if (gridX < 0) {gridX = 0; skipX = true;}
+        if (gridY < 0) {gridY = 0; skipY = true;}
         int ix = (int) gridX;  gridX -= ix;
         int iy = (int) gridY;  gridY -= iy;
         int n;
         if (ix > (n = gridSize[0] - 2)) {
-            ix = n;
-            gridX = 1;
+            skipX |= (ix != n+1 || gridX != 0);         // Keep value 'false' if gridX ==
gridSize[0] - 1.
+            ix     = n;
+            gridX  = 1;
         }
         if (iy > (n = gridSize[1] - 2)) {
-            iy = n;
-            gridY = 1;
+            skipY |= (iy != n+1 || gridY != 0);         // Keep value 'false' if gridY ==
gridSize[1] - 1.
+            iy     = n;
+            gridY  = 1;
         }
         n = getTranslationDimensions();
         boolean derivative = (vector.length >= n + INTERPOLATED_DIMENSIONS * INTERPOLATED_DIMENSIONS);
         for (int dim = 0; dim < n; dim++) {
-            double dx;
+            double dx, dy;
             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);
+            final double r0x = r00 + gridX * (dx = r01 - r00);      // TODO: use Math.fma
on JDK9.
+            final double r1x = r10 + gridX * (dy = r11 - r10);      // Not really "dy" measurement
yet, will become dy later.
             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;
+                if (skipX) {
+                    dx = 0;
+                } else {
+                    dx += (dy - dx) * gridX;
+                }
+                if (skipY) {
+                    dy = 0;
+                } else {
+                    dy  =  r10 - r00;
+                    dy += (r11 - r01 - dy) * gridY;
+                }
                 int i = n;
                 if (dim == 0) {
                     dx++;
@@ -490,14 +504,18 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
     }
 
     /**
-     * Returns the derivative at the given grid indices.
-     * If the given coordinates is outside the grid, then this method returns the derivative
at the closest cell.
+     * Estimates the derivative at the given grid indices.
      *
-     * <div class="section">Default implementation</div>
-     * The current implementation assumes that the derivative is constant everywhere in the
cell
-     * at the given indices. It does not yet take in account the fractional part of {@code
gridX}
-     * and {@code gridY}, because empirical tests suggest that the accuracy of such interpolation
-     * is uncertain.
+     * <div class="section">Extrapolations</div>
+     * If the given coordinates is outside the grid, then the derivative will have some columns
set to identity.
+     *
+     * <ul>
+     *   <li>If both {@code gridX} and {@code gridY} are outside the grid, then the
derivative is the identity matrix.</li>
+     *   <li>If only {@code gridX} is outside the grid, then only the first column
is set to [1, 0, …].
+     *       The second column is set to the derivative of the closest cell at {@code gridY}
position.</li>
+     *   <li>If only {@code gridY} is outside the grid, then only the second column
is set to [0, 1, …].
+     *       The first column is set to the derivative of the closest cell at {@code gridX}
position.</li>
+     * </ul>
      *
      * @param  gridX  first grid coordinate of the point for which to get the translation.
      * @param  gridY  second grid coordinate of the point for which to get the translation.
@@ -506,16 +524,29 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
      * @see #isCellInGrid(double, double)
      * @see #interpolateInCell(double, double, double[])
      */
-    public Matrix derivativeInCell(final double gridX, final double gridY) {
+    public Matrix derivativeInCell(double gridX, double gridY) {
         final int ix = Math.max(0, Math.min(gridSize[0] - 2, (int) gridX));
         final int iy = Math.max(0, Math.min(gridSize[1] - 2, (int) gridY));
+        gridX -= ix;
+        gridY -= iy;
+        final boolean skipX = (gridX < 0 || gridX > 1);
+        final boolean skipY = (gridY < 0 || gridY > 1);
         final Matrix derivative = Matrices.createDiagonal(getTranslationDimensions(), gridSize.length);
         for (int j=derivative.getNumRow(); --j>=0;) {
             final double r00 = getCellValue(j, ix,   iy  );
-            final double 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);
+            final double r01 = getCellValue(j, ix+1, iy  );       // Naming convention: ryx
(row index first, like matrix).
+            final double r10 = getCellValue(j, ix,   iy+1);
+            final double r11 = getCellValue(j, ix+1, iy+1);
+            if (!skipX) {
+                double dx = r01 - r00;
+                dx += (r11 - r10 - dx) * gridX;
+                derivative.setElement(j, 0, derivative.getElement(j, 0) + dx);
+            }
+            if (!skipY) {
+                double dy = r10 - r00;
+                dy += (r11 - r01 - dy) * gridY;
+                derivative.setElement(j, 1, derivative.getElement(j, 1) + dy);
+            }
         }
         return derivative;
     }
@@ -533,11 +564,16 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
      *       {@linkplain #DatumShiftGrid(Unit, LinearTransform, int[], boolean, Unit) constructor
javadoc}).</li>
      * </ul>
      *
+     * Caller must ensure that all arguments given to this method are in their expected ranges.
+     * The behavior of this method is undefined if any argument value is out-of-range.
+     * (this method is not required to validate arguments, for performance reasons).
+     *
      * @param  dim    the dimension of the translation vector component to get,
      *                from 0 inclusive to {@link #getTranslationDimensions()} exclusive.
      * @param  gridX  the grid index on the <var>x</var> axis, from 0 inclusive
to {@code gridSize[0]} exclusive.
      * @param  gridY  the grid index on the <var>y</var> axis, from 0 inclusive
to {@code gridSize[1]} exclusive.
      * @return the translation for the given dimension in the grid cell at the given index.
+     * @throws IndexOutOfBoundsException may be thrown (but is not guaranteed to be throw)
if an argument is out of range.
      */
     public abstract double getCellValue(int dim, int gridX, int gridY);
 
@@ -685,9 +721,13 @@ public abstract class DatumShiftGrid<C extends Quantity<C>,
T extends Quantity<T
      */
     @Override
     public String toString() {
-        final Parameters p = Parameters.castOrWrap(getParameterDescriptors().createValue());
-        getParameterValues(p);
-        return p.toString();
+        final ParameterDescriptorGroup d = getParameterDescriptors();
+        if (d != null) {
+            final Parameters p = Parameters.castOrWrap(d.createValue());
+            getParameterValues(p);
+            return p.toString();
+        }
+        return super.toString();
     }
 
     /**
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
index cad1710..2c30927 100644
--- 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
@@ -20,6 +20,7 @@ import java.util.Random;
 import java.awt.geom.Point2D;
 import java.awt.geom.AffineTransform;
 import javax.measure.quantity.Dimensionless;
+import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.TransformException;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.apache.sis.measure.Units;
@@ -172,10 +173,11 @@ public strictfp class DatumShiftGridFileTest extends TestCase {
         init(30);
         for (int i=0; i<NUM_TESTS; i++) {
             /*
-             * Compute the reference point.
+             * Compute the reference point. We need to avoid points outside the grid because
derivates at those
+             * locations are partially fixed to identity, which is different than affine
transform coefficients.
              */
-            final int x = random.nextInt(WIDTH);
-            final int y = random.nextInt(HEIGHT);
+            final int x = random.nextInt(WIDTH  - 1);
+            final int y = random.nextInt(HEIGHT - 1);
             point.x = x;
             point.y = y;
             assertSame(point, reference.transform(point, point));
@@ -189,6 +191,50 @@ public strictfp class DatumShiftGridFileTest extends TestCase {
             assertEquals("m01", reference.getShearX(), vector[3], TOLERANCE);
             assertEquals("m10", reference.getShearY(), vector[4], TOLERANCE);
             assertEquals("m11", reference.getScaleY(), vector[5], TOLERANCE);
+            assertSameDerivative(x, y, vector);
+        }
+    }
+
+    /**
+     * Tests {@link DatumShiftGridFile#interpolateAt(double...)} with some values outside
the grid.
+     * Derivatives outside the grid have different coefficients than derivatives inside the
grid.
+     * This test verifies that methods computing derivatives are self-consistent.
+     *
+     * @throws TransformException if an error occurred while transforming a coordinates.
+     */
+    @Test
+    @DependsOnMethod("testInterpolateAndDerivative")
+    public void testExtrapolation() throws TransformException {
+        final Random random = TestUtilities.createRandomNumberGenerator();
+        final Point2D.Float point = new Point2D.Float();
+        final double[] vector = new double[6];
+        init(50);
+        for (int i=0; i<NUM_TESTS; i++) {
+            final int x = random.nextInt(WIDTH  * 2) - WIDTH  / 4;
+            final int y = random.nextInt(HEIGHT * 2) - HEIGHT / 4;
+            point.x = x;
+            point.y = y;
+            grid.interpolateInCell(x, y, vector);
+            if (x >= 0 && x < WIDTH && y >= 0 && y <
HEIGHT) {
+                assertSame(point, reference.transform(point, point));
+                assertEquals("x", point.x, (float) (vector[0] + x), TOLERANCE);
+                assertEquals("y", point.y, (float) (vector[1] + y), TOLERANCE);
+            }
+            assertSameDerivative(x, y, vector);
         }
     }
+
+    /**
+     * Verifies that the matrix returned by {@link DatumShiftGridFile#derivativeInCell(double,
double)}
+     * contains coefficients identical to the ones in the given vector.
+     */
+    private void assertSameDerivative(final int x, final int y, final double[] vector) {
+        Matrix m = grid.derivativeInCell(x, y);
+        assertEquals("numRow", 2, m.getNumRow());
+        assertEquals("numCol", 2, m.getNumCol());
+        assertEquals("m00", m.getElement(0,0), vector[2], STRICT);
+        assertEquals("m01", m.getElement(0,1), vector[3], STRICT);
+        assertEquals("m10", m.getElement(1,0), vector[4], STRICT);
+        assertEquals("m11", m.getElement(1,1), vector[5], STRICT);
+    }
 }
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 fd1e2c7..e673b04 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
@@ -150,7 +150,7 @@ public final strictfp class InterpolatedTransformTest extends MathTransformTestC
     }
 
     /**
-     * Tests the derivatives at the sample point. This method compares the derivatives computed
by
+     * Tests the derivatives at the sample points. This method compares the derivatives computed
by
      * the transform with an estimation of derivatives computed by the finite differences
method.
      *
      * @throws TransformException if an error occurred while transforming the coordinate.
@@ -168,7 +168,7 @@ public final strictfp class InterpolatedTransformTest extends MathTransformTestC
              * for the same reason than in `testForwardTransform()`. The matrix values are
close to ±1.
              */
             System.arraycopy(samplePoints[0], i, point, 0, SinusoidalShiftGrid.DIMENSION);
-            tolerance = (i < SinusoidalShiftGrid.FIRST_FRACTIONAL_COORDINATE) ? 1E-10
: 0.02;
+            tolerance = (i < SinusoidalShiftGrid.FIRST_FRACTIONAL_COORDINATE) ? 1E-10
: 0.01;
             verifyDerivative(point);
             /*
              * Verify derivative at the same point but using inverse transform,
@@ -184,6 +184,34 @@ public final strictfp class InterpolatedTransformTest extends MathTransformTestC
     }
 
     /**
+     * Tests the derivatives at sample points outside the grid. Those derivatives must be
consistent
+     * in order to allow inverse transformation to work when the initial point is outside
the grid.
+     *
+     * @throws TransformException if an error occurred while transforming the coordinate.
+     */
+    @Test
+    @DependsOnMethod("testDerivative")
+    public void testExtrapolations() throws TransformException {
+        createSinusoidal(-50);
+        final double[] point = new double[SinusoidalShiftGrid.DIMENSION];
+        derivativeDeltas = new double[] {0.002, 0.002};
+        isInverseTransformSupported = false;                                            
   // For focusing on a single aspect.
+        tolerance = 1E-10;
+        for (int i=0; i<=4; i++) {
+            switch (i) {
+                default: throw new AssertionError(i);
+                case 0: point[0] = -50; point[1] =  40; break;       // Point outside grid
on the left.
+                case 1: point[0] = 200; point[1] =  60; break;       // Point outside grid
on the right.
+                case 2: point[0] =  20; point[1] = -50; break;       // Point outside grid
on the top.
+                case 3: point[0] = -80; point[1] = 230; break;       // Point outside grid
two sides.
+                case 4: point[0] =  80; point[1] = 185;              // Point outside grid
on the bottom.
+                        tolerance = 0.3; break;
+            }
+            verifyDerivative(point);
+        }
+    }
+
+    /**
      * 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.
      *
@@ -209,7 +237,6 @@ public final strictfp class InterpolatedTransformTest extends MathTransformTestC
         // Forward derivative
         transform        = transform.inverse();
         derivativeDeltas = new double[] {0.2, 0.2};
-        tolerance        = 5E-6;                        // Empirical value.
         verifyDerivative(FranceGeocentricInterpolationTest.samplePoint(1));
 
         // Inverse derivative
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/SinusoidalShiftGrid.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/SinusoidalShiftGrid.java
index c88ad6e..1a8e93e 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/SinusoidalShiftGrid.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/SinusoidalShiftGrid.java
@@ -25,6 +25,8 @@ import org.apache.sis.measure.Units;
 import org.apache.sis.parameter.DefaultParameterDescriptorGroup;
 import org.apache.sis.parameter.Parameters;
 
+import static org.junit.Assert.*;
+
 
 /**
  * Dummy implementation of {@link DatumShiftGrid} containing translation vectors that are
computed by a sinusoidal function.
@@ -136,6 +138,9 @@ final strictfp class SinusoidalShiftGrid extends DatumShiftGrid<Dimensionless,Di
      */
     @Override
     public double getCellValue(int dim, int gridX, int gridY) {
+        assertTrue(dim   >= 0 && dim < DIMENSION);
+        assertTrue(gridX >= 0 && gridX < WIDTH);
+        assertTrue(gridY >= 0 && gridY < HEIGHT);
         buffer[0] = gridX;
         buffer[1] = gridY;
         transform(buffer, buffer);


Mime
View raw message