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 <[email protected]>
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);
