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: Enable InterpolatedTransform.Inverse method based on Jacobian matrix (but still needs more test).
Date Sun, 14 Apr 2019 17:36:42 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 4685195  Enable InterpolatedTransform.Inverse method based on Jacobian matrix (but
still needs more test).
4685195 is described below

commit 468519583b54f5499b0c11c8799a272a629ee966
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
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);
     }
 
     /**


Mime
View raw message