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 4685195 Enable InterpolatedTransform.Inverse method based on Jacobian
matrix (but still needs more test).
4685195 is described below
commit 468519583b54f5499b0c11c8799a272a629ee966
Author: Martin Desruisseaux <[email protected]>
AuthorDate: Sun Apr 14 17:15:29 2019 +0200
Enable InterpolatedTransform.Inverse method based on Jacobian matrix (but
still needs more test).
---
.../operation/transform/InterpolatedTransform.java | 108 ++++++++++-----------
.../transform/InterpolatedTransformTest.java | 55 ++++++-----
...aticShiftGrid.java => SinusoidalShiftGrid.java} | 36 +++----
3 files changed, 102 insertions(+), 97 deletions(-)
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 62ccf8b..0f3d4cd 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
@@ -427,7 +427,7 @@ public class InterpolatedTransform extends
DatumShiftTransform {
* more curved grids like the ones read from netCDF files. We provide
this switch for allowing performance
* comparisons.
*/
- private static final boolean SIMPLE = true;
+ private static final boolean SIMPLE = false;
/**
* Maximum number of iterations. This is set to a higher value than
{@link Formulas#MAXIMUM_ITERATIONS} because
@@ -438,8 +438,6 @@ public class InterpolatedTransform extends
DatumShiftTransform {
* values, for example close to 1000 while we usually expect values
smaller than 1. Behavior with such grids may
* be unpredictable, sometime with the {@code abs(xi - ox)} or {@code
abs(yi - oy)} errors staying high for a
* long time before to suddenly fall to zero.
- *
- * @see #tryAgain(int, double, double)
*/
private static final int MAXIMUM_ITERATIONS =
Formulas.MAXIMUM_ITERATIONS * 4;
@@ -527,18 +525,16 @@ public class InterpolatedTransform extends
DatumShiftTransform {
}
}
final double[] vector = new double[SIMPLE ? dimension : dimension
+ 4]; // +4 is for derivate coefficients.
-nextPoint: while (--numPts >= 0) {
+ while (--numPts >= 0) {
double xi, yi;
// (x,y) values at iteration i.
final double x = xi = srcPts[srcOff ];
final double y = yi = srcPts[srcOff+1];
+ double tol = tolerance;
if (DEBUG) {
System.out.format("Searching source coordinates for target
= (%g %g)%n", x, y);
}
- int it = MAXIMUM_ITERATIONS;
- double tol = tolerance;
- do {
+ for (int it = MAXIMUM_ITERATIONS;;) {
forward.grid.interpolateInCell(xi, yi, vector);
- final boolean more;
if (SIMPLE) {
/*
* We want (xi, yi) such as the following conditions
hold
@@ -551,7 +547,7 @@ nextPoint: while (--numPts >= 0) {
final double oy = yi;
xi = x - vector[0];
yi = y - vector[1];
- more = Math.abs(xi - ox) > tol || Math.abs(yi - oy) >
tol;
+ if (!(Math.abs(xi - ox) > tol || Math.abs(yi - oy) >
tol)) break; // Use '!' for catching NaN.
} else {
/*
* The error between the new position (xi + tx) and
the desired position x is measured
@@ -588,65 +584,63 @@ nextPoint: while (--numPts >= 0) {
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",
+ System.out.printf(" source=(%9.3f %9.3f)
target=(%9.3f %9.3f) error=(%11.6f %11.6f) → (%11.6f %11.6f)%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 (!(Math.abs(ex) > tol || Math.abs(ey) > tol))
break; // Use '!' for catching NaN.
}
- if (!more) {
// Use '!' for catching NaN.
- 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);
+ /*
+ * At this point we determined that we need to iterate
more. If iteration does not converge, we may relaxe
+ * threshold in last resort but nevertheless aim for an
accuracy of 0.5 of cell size in order to keep some
+ * consistency with forward transform. If the point was
inside the grid, we assume (for well-formed grid)
+ * that iteration should have converged. But during
extrapolations since there is no authoritative results,
+ * we consider that a more approximate result is okay. In
particular it does not make sense to require a
+ * 1E-7 accuracy (relative to cell size) if we don't
really know what the answer should be.
+ */
+ if (--it < 0) {
+ if (it < -1) {
+ xi = yi = Double.NaN;
+ break;
+ }
+ if (forward.grid.isCellInGrid(xi, yi)) {
+ throw new
TransformException(Resources.format(Resources.Keys.NoConvergence));
}
- dstPts[dstOff ] = xi; // Shall not be done
before above loop.
- dstPts[dstOff+1] = yi;
- dstOff += inc;
- srcOff += inc;
- continue nextPoint;
+ tol = 0.5;
}
- } while (--it >= 0 || (tol = tryAgain(it, xi, yi)) > 0);
- throw new
TransformException(Resources.format(Resources.Keys.NoConvergence));
- }
- }
-
- /**
- * If iteration did not converge, tells whether we should perform
another try with a more permissive threshold.
- * We start relaxing threshold only in last resort, and nevertheless
aim for an accuracy of 0.5 of cell size in
- * order to keep some consistency with forward transform. We may relax
more in case of extrapolations.
- *
- * @param it the iteration counter. Should be negative since we
exhausted the normal number of iterations.
- * @param xi best <var>x</var> estimation so far.
- * @param yi best <var>y</var> estimation so far.
- * @return the new tolerance threshold, or {@link Double#NaN} if no
more try should be allowed.
- *
- * @see #MAXIMUM_ITERATIONS
- */
- private double tryAgain(final int it, final double xi, final double
yi) {
- double tol = Math.scalb(tolerance, -it);
- if (tol >= 0.5) {
+ }
/*
- * If the point was inside the grid and the grid is
well-formed, we assume that iteration should have converged.
- * But during extrapolations since there is no authoritative
results, we consider that a more approximate result
- * is okay. In particular it does not make sense to require a
1E-7 accuracy (relative to cell size) if we don't
- * really know what the answer should be.
+ * At this point iteration converged. Store the result before
to continue with next point.
*/
- if (forward.grid.isCellInGrid(xi, yi) || tol > 2) {
- return Double.NaN; //
No more iteration - caller will throw an exception.
+ if (DEBUG) {
+ if (!Double.isNaN(xi) && !Double.isNaN(yi)) {
+ forward.grid.interpolateInCell(xi, yi, vector);
+ final double xf = xi + vector[0];
+ final double yf = yi + vector[1];
+ assert Math.abs(xf - x) <= tol : xf;
+ assert Math.abs(yf - y) <= tol : yf;
+ }
+ }
+ 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;
+ dstOff += inc;
+ srcOff += inc;
}
- return tol;
}
}
}
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 a983095..fd1e2c7 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
@@ -57,16 +57,16 @@ 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,
+ * Creates an {@link InterpolatedTransform} derived from a sinusoidal
formula.
+ * We do not really need {@code InterpolatedTransform} for sinusoidal
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);
+ private double[][] createSinusoidal(final double rotation) throws
TransformException {
+ final SinusoidalShiftGrid grid = new SinusoidalShiftGrid(rotation);
transform = new InterpolatedTransform(grid);
return grid.samplePoints();
}
@@ -112,12 +112,16 @@ public final strictfp class InterpolatedTransformTest
extends MathTransformTestC
@Test
public void testForwardTransform() throws TransformException {
isInverseTransformSupported = false;
// For focusing on a single aspect.
- final double[][] samplePoints = createQuadratic(-15);
+ final double[][] samplePoints = createSinusoidal(-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(Arrays.copyOf(samplePoints[0],
SinusoidalShiftGrid.FIRST_FRACTIONAL_COORDINATE),
+ Arrays.copyOf(samplePoints[1],
SinusoidalShiftGrid.FIRST_FRACTIONAL_COORDINATE));
+ /*
+ * For non-integer coordinates, we need to relax the tolerance
threshold because the linear interpolations
+ * computed by InterpolatedTransform do not give the same results than
the calculation done with cosine by
+ * SinudoisalShiftGrid. The result of tested point is about (81.96
22.89).
+ */
+ tolerance = 0.01;
verifyTransform(samplePoints[0], samplePoints[1]);
}
@@ -128,13 +132,20 @@ public final strictfp class InterpolatedTransformTest
extends MathTransformTestC
* @throws TransformException if an error occurred while transforming a
coordinate.
*/
@Test
- @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);
+ final double[][] samplePoints = createSinusoidal(-20);
transform = transform.inverse();
- tolerance = QuadraticShiftGrid.PRECISION;
+ tolerance = SinusoidalShiftGrid.PRECISION;
+ verifyTransform(Arrays.copyOf(samplePoints[1],
SinusoidalShiftGrid.FIRST_FRACTIONAL_COORDINATE),
+ Arrays.copyOf(samplePoints[0],
SinusoidalShiftGrid.FIRST_FRACTIONAL_COORDINATE));
+ /*
+ * For non-integer coordinates, we need to relax the tolerance
threshold because the linear interpolations
+ * computed by InterpolatedTransform do not give the same results than
the calculation done with cosine by
+ * SinudoisalShiftGrid.
+ */
+ tolerance = 0.01;
verifyTransform(samplePoints[1], samplePoints[0]);
}
@@ -147,17 +158,17 @@ public final strictfp class InterpolatedTransformTest
extends MathTransformTestC
@Test
@DependsOnMethod("testForwardTransform")
public void testDerivative() throws TransformException {
- final double[][] samplePoints = createQuadratic(-40);
- final double[] point = new double[QuadraticShiftGrid.DIMENSION];
// A single point from 'samplePoints'
+ final double[][] samplePoints = createSinusoidal(-40);
+ final double[] point = new double[SinusoidalShiftGrid.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.
- }
+ for (int i=0; i < samplePoints[0].length; i +=
SinusoidalShiftGrid.DIMENSION) {
+ /*
+ * The tolerance threshold must be relaxed for derivative at a
position having factional digits
+ * 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;
verifyDerivative(point);
/*
* Verify derivative at the same point but using inverse transform,
@@ -165,7 +176,7 @@ public final strictfp class InterpolatedTransformTest
extends MathTransformTestC
*/
if (isInverseTransformSupported) {
transform = transform.inverse();
- System.arraycopy(samplePoints[1], i, point, 0,
QuadraticShiftGrid.DIMENSION);
+ System.arraycopy(samplePoints[1], i, point, 0,
SinusoidalShiftGrid.DIMENSION);
verifyDerivative(point);
transform = transform.inverse(); // Back to forward
transform.
}
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/SinusoidalShiftGrid.java
similarity index 84%
rename from
core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/QuadraticShiftGrid.java
rename to
core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/transform/SinusoidalShiftGrid.java
index 6eb20df..c88ad6e 100644
---
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/SinusoidalShiftGrid.java
@@ -27,8 +27,8 @@ import org.apache.sis.parameter.Parameters;
/**
- * 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
+ * Dummy implementation of {@link DatumShiftGrid} containing translation
vectors that are computed by a sinusoidal function.
+ * This class has no computational interest compared to a direct
implementation of sinusoidal formulas, but is convenient
* for debugging {@link InterpolatedTransform} because of its predictable
behavior (easier to see with rotation of 0°).
*
* @author Martin Desruisseaux (Geomatys)
@@ -37,7 +37,7 @@ import org.apache.sis.parameter.Parameters;
* @module
*/
@SuppressWarnings("serial") // Not intended to be
serialized.
-final strictfp class QuadraticShiftGrid extends
DatumShiftGrid<Dimensionless,Dimensionless> {
+final strictfp class SinusoidalShiftGrid extends
DatumShiftGrid<Dimensionless,Dimensionless> {
/**
* Number of source and target dimensions of the grid.
*/
@@ -60,8 +60,9 @@ final strictfp class QuadraticShiftGrid extends
DatumShiftGrid<Dimensionless,Dim
/**
* Grid size in number of pixels. The translation vectors in the middle of
this grid will be (0,0).
+ * We simulate a grid spanning from 80°S to 80°N.
*/
- private static final int WIDTH = 200, HEIGHT = 2000;
+ private static final int WIDTH = 100, HEIGHT = 160;
/**
* An internal step performed during computation of translation vectors at
a given point.
@@ -80,11 +81,11 @@ final strictfp class QuadraticShiftGrid extends
DatumShiftGrid<Dimensionless,Dim
*
* @param rotation the rotation angle, in degrees.
*/
- QuadraticShiftGrid(final double rotation) {
+ SinusoidalShiftGrid(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);
+ offsets = AffineTransform.getTranslateInstance(0.5*WIDTH, 0.5*HEIGHT);
+ offsets.rotate(StrictMath.toRadians(rotation));
+ offsets.scale(-0.75, 0.95);
buffer = new double[DIMENSION];
}
@@ -105,8 +106,8 @@ final strictfp class QuadraticShiftGrid extends
DatumShiftGrid<Dimensionless,Dim
final double[][] samplePoints() {
final double[] sources = {
/*[0]*/ WIDTH/2, HEIGHT/2,
- /*[1]*/ 50, 1400,
- /*[2]*/ 175, 200,
+ /*[1]*/ 30, 120,
+ /*[2]*/ 75, 40,
/*[3]*/ 10.356, 30.642 //
FIRST_FRACTIONAL_COORDINATE must point here.
};
final double[] targets = new double[sources.length];
@@ -119,15 +120,14 @@ final strictfp class QuadraticShiftGrid extends
DatumShiftGrid<Dimensionless,Dim
* 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;
+ for (int i=0; i<sources.length;) {
+ double x = sources[i ] - WIDTH / 2;
+ double y = sources[i+1] - HEIGHT / 2;
+ x /= StrictMath.max(0.1, StrictMath.cos(StrictMath.toRadians(y)));
+ targets[i++] = x;
+ targets[i++] = y;
}
+ offsets.transform(targets, 0, targets, 0, targets.length / DIMENSION);
}
/**
