sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 03/03: Faster computation of LocalizationGrid, and fix a bug in the WKT representation of grid values as matrices.
Date Wed, 06 Feb 2019 20:12:19 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 8c8409ccade9b377b1910b9fd515623c98e10182
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Feb 6 21:11:01 2019 +0100

    Faster computation of LocalizationGrid, and fix a bug in the WKT representation of grid
values as matrices.
---
 .../operation/builder/LinearTransformBuilder.java  | 17 +++--
 .../operation/builder/LocalizationGridBuilder.java | 80 +++++++++++++++-------
 .../operation/builder/ResidualGrid.java            | 65 +++++-------------
 .../operation/builder/ResidualGridTest.java        |  2 +-
 4 files changed, 86 insertions(+), 78 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
index c5e5c90..972f768 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java
@@ -897,13 +897,22 @@ search:         for (int j=domain(); --j >= 0;) {
      * More straightforward version of {@link #getControlPoint(int[])} for the case where
this
      * {@code LinearTransformBuilder} is known to have been built for grid source coordinates.
      * This method is for {@link LocalizationGridBuilder#create(MathTransformFactory)} internal
usage.
+     *
+     * @param  source  the source coordinates.
+     * @param  target  where to store target coordinates for a full row.
      */
-    final void getControlPoint2D(final int[] source, final double[] target) {
+    final void getControlRow(final int[] source, final double[] target) {
         assert gridSize != null;
-        final int index = flatIndex(source);
+        final int start  = flatIndex(source);
+        final int stop   = start + gridSize[0];
         final int tgtDim = targets.length;
-        for (int i=0; i<tgtDim; i++) {
-            target[i] = targets[i][index];
+        for (int j=0; j<tgtDim; j++) {
+            int index = j;
+            final double[] row = targets[j];
+            for (int i=start; i<stop; i++) {
+                target[index] = row[i];
+                index += tgtDim;
+            }
         }
     }
 
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
index 60a2583..f940417 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LocalizationGridBuilder.java
@@ -27,19 +27,20 @@ import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
 import org.apache.sis.referencing.operation.transform.LinearTransform;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
-import org.apache.sis.referencing.operation.matrix.MatrixSIS;
 import org.apache.sis.referencing.operation.matrix.Matrix3;
 import org.apache.sis.referencing.datum.DatumShiftGrid;
 import org.apache.sis.internal.referencing.Resources;
-import org.apache.sis.geometry.DirectPosition2D;
 import org.apache.sis.geometry.GeneralEnvelope;
 import org.apache.sis.geometry.Envelopes;
 import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.measure.NumberRange;
 import org.apache.sis.math.MathFunctions;
+import org.apache.sis.math.Statistics;
 import org.apache.sis.math.Vector;
 
+import static org.apache.sis.referencing.operation.builder.ResidualGrid.SOURCE_DIMENSION;
+
 
 /**
  * Creates an "almost linear" transform mapping the given source points to the given target
points.
@@ -86,7 +87,7 @@ public class LocalizationGridBuilder extends TransformBuilder {
      * A temporary array for two-dimensional source coordinates.
      * Used for reducing object allocations.
      */
-    private final int[] tmp = new int[2];
+    private final int[] tmp = new int[SOURCE_DIMENSION];
 
     /**
      * Conversions from source real-world coordinates to grid indices before interpolation.
@@ -118,7 +119,7 @@ public class LocalizationGridBuilder extends TransformBuilder {
      */
     public LocalizationGridBuilder(final int width, final int height) {
         linear = new LinearTransformBuilder(width, height);
-        sourceToGrid = MathTransforms.identity(2);
+        sourceToGrid = MathTransforms.identity(SOURCE_DIMENSION);
     }
 
     /**
@@ -175,14 +176,14 @@ public class LocalizationGridBuilder extends TransformBuilder {
     public LocalizationGridBuilder(final LinearTransformBuilder localizations) {
         ArgumentChecks.ensureNonNull("localizations", localizations);
         int n = localizations.getGridDimensions();
-        if (n == 2) {
+        if (n == SOURCE_DIMENSION) {
             linear = localizations;
-            sourceToGrid = MathTransforms.identity(2);
+            sourceToGrid = MathTransforms.identity(SOURCE_DIMENSION);
         } else {
             if (n < 0) {
                 final Vector[] sources = localizations.sources();
                 n = sources.length;
-                if (n == 2) {
+                if (n == SOURCE_DIMENSION) {
                     final Matrix fromGrid = new Matrix3();
                     final int width  = infer(sources[0], fromGrid, 0);
                     final int height = infer(sources[1], fromGrid, 1);
@@ -196,7 +197,8 @@ public class LocalizationGridBuilder extends TransformBuilder {
                     return;
                 }
             }
-            throw new IllegalArgumentException(Resources.format(Resources.Keys.MismatchedTransformDimension_3,
0, 2, n));
+            throw new IllegalArgumentException(Resources.format(
+                    Resources.Keys.MismatchedTransformDimension_3, 0, SOURCE_DIMENSION, n));
         }
     }
 
@@ -240,7 +242,7 @@ public class LocalizationGridBuilder extends TransformBuilder {
          * for preventing useless allocation of large arrays - https://issues.apache.org/jira/browse/SIS-407
          */
         fromGrid.setElement(dim, dim, inc);
-        fromGrid.setElement(dim,   2, min);
+        fromGrid.setElement(dim, SOURCE_DIMENSION, min);
         final double n = span / inc;
         if (n >= 0.5 && n < source.size() - 0.5) {          // Compare as 'double'
in case the value is large.
             return ((int) Math.round(n)) + 1;
@@ -318,16 +320,16 @@ public class LocalizationGridBuilder extends TransformBuilder {
         ArgumentChecks.ensureNonNull("sourceToGrid", sourceToGrid);
         int isTarget = 0;
         int dim = sourceToGrid.getSourceDimensions();
-        if (dim >= 2) {
+        if (dim >= SOURCE_DIMENSION) {
             isTarget = 1;
             dim = sourceToGrid.getTargetDimensions();
-            if (dim == 2) {
+            if (dim == SOURCE_DIMENSION) {
                 this.sourceToGrid = sourceToGrid;
                 return;
             }
         }
         throw new MismatchedDimensionException(Resources.format(
-                Resources.Keys.MismatchedTransformDimension_3, isTarget, 2, dim));
+                Resources.Keys.MismatchedTransformDimension_3, isTarget, SOURCE_DIMENSION,
dim));
     }
 
     /**
@@ -449,9 +451,8 @@ public class LocalizationGridBuilder extends TransformBuilder {
         }
         final int      width    = linear.gridSize(0);
         final int      height   = linear.gridSize(1);
-        final int      tgtDim   = gridToCoord.getTargetDimensions();
-        final double[] residual = new double[tgtDim * linear.gridLength];
-        final double[] point    = new double[tgtDim + 1];
+        final double[] residual = new double[SOURCE_DIMENSION * linear.gridLength];
+        final double[] grid     = new double[SOURCE_DIMENSION * width];
         double gridPrecision    = precision;
         try {
             /*
@@ -477,19 +478,15 @@ public class LocalizationGridBuilder extends TransformBuilder {
              * than the desired precision,  then the returned MathTransform will need to
apply corrections
              * after linear transforms. Those corrections will be done by InterpolatedTransform.
              */
-            final MatrixSIS coordToGrid = MatrixSIS.castOrCopy(gridToCoord.inverse().getMatrix());
-            final DirectPosition2D src = new DirectPosition2D();
-            point[tgtDim] = 1;
+            final MathTransform coordToGrid = gridToCoord.inverse();
             for (int k=0,y=0; y<height; y++) {
-                src.y  = y;
+                tmp[0] = 0;
                 tmp[1] = y;
+                linear.getControlRow(tmp, grid);                                    // Expected
positions.
+                coordToGrid.transform(grid, 0, residual, k, width);                 // As
grid coordinate.
                 for (int x=0; x<width; x++) {
-                    src.x  = x;
-                    tmp[0] = x;
-                    linear.getControlPoint2D(tmp, point);                           // Expected
position.
-                    double[] grid = coordToGrid.multiply(point);                    // As
grid coordinate.
-                    isLinear &= (residual[k++] = grid[0] - x) <= gridPrecision;
-                    isLinear &= (residual[k++] = grid[1] - y) <= gridPrecision;
+                    isLinear &= (residual[k++] -= x) <= gridPrecision;
+                    isLinear &= (residual[k++] -= y) <= gridPrecision;
                 }
             }
         } catch (TransformException e) {
@@ -499,7 +496,38 @@ public class LocalizationGridBuilder extends TransformBuilder {
             return MathTransforms.concatenate(sourceToGrid, gridToCoord);
         }
         return InterpolatedTransform.createGeodeticTransformation(nonNull(factory),
-                new ResidualGrid(sourceToGrid, gridToCoord, width, height, tgtDim, residual,
+                new ResidualGrid(sourceToGrid, gridToCoord, width, height, residual,
                 (gridPrecision > 0) ? gridPrecision : DEFAULT_PRECISION));
     }
+
+    /**
+     * Returns statistics of differences between values calculated by the given transform
and actual values.
+     * The given math transform is typically the transform computed by {@link #create(MathTransformFactory)},
+     * but not necessarily.
+     *
+     * @param  mt  the transform to test.
+     * @return statistics of difference between computed values and expected value.
+     * @throws TransformException if an error occurred while transforming a coordinate.
+     *
+     * @since 1.0
+     */
+    public Statistics error(final MathTransform mt) throws TransformException {
+        final Statistics s = new Statistics(null);
+        final int width  = linear.gridSize(0);
+        final int height = linear.gridSize(1);
+        final int tgtDim = Math.max(mt.getTargetDimensions(), SOURCE_DIMENSION);
+        final double[] t = new double[tgtDim];
+        for (int y=0; y<height; y++) {
+            for (int x=0; x<width; x++) {
+                t[0] = tmp[0] = x;
+                t[1] = tmp[1] = y;
+                mt.transform(t, 0, t, 0, 1);
+                final double[] expected = linear.getControlPoint(tmp);
+                for (int i=0; i<tgtDim; i++) {
+                    s.accept(t[i] - expected[i]);
+                }
+            }
+        }
+        return s;
+    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
index 6c53cc3..b6de64a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/ResidualGrid.java
@@ -49,13 +49,17 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
     /**
      * For cross-version compatibility.
      */
-    private static final long serialVersionUID = -6661539177698674636L;
+    private static final long serialVersionUID = 1445697681304159019L;
+
+    /**
+     * Number of source dimensions of the residual grid.
+     */
+    static final int SOURCE_DIMENSION = 2;
 
     /**
      * The parameter descriptors for the "Localization grid" operation.
-     * Current implementation accepts a maximum of 3 dimensions; if the residual grid contains
more dimensions,
-     * the extra dimensions will not be shown in the parameters. This is not a committed
set of parameters and
-     * they may change in any future SIS version. We define them mostly for {@code toString()}
implementation.
+     * Current implementation is fixed to 2 dimensions. This is not a committed set of parameters
and they
+     * may change in any future SIS version. We define them mostly for {@code toString()}
implementation.
      */
     private static final ParameterDescriptorGroup PARAMETERS;
     static {
@@ -64,10 +68,8 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
         final ParameterDescriptor<?>[] grids = new ParameterDescriptor[] {
             builder.addName(Constants.NUM_ROW).createBounded(Integer.class, 2, null, null),
             builder.addName(Constants.NUM_COL).createBounded(Integer.class, 2, null, null),
-            builder.addName("num_dim").createBounded(Integer.class, 1, null, 2),
             builder.addName("grid_x").create(Matrix.class, null),
-            builder.addName("grid_y").setRequired(false).create(Matrix.class, null),
-            builder.addName("grid_z").create(Matrix.class, null)
+            builder.addName("grid_y").create(Matrix.class, null)
         };
         PARAMETERS = builder.addName("Localization grid").createGroup(grids);
     }
@@ -89,17 +91,8 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
         final int[] size = getGridSize();
         parameters.parameter(Constants.NUM_ROW).setValue(size[1]);
         parameters.parameter(Constants.NUM_COL).setValue(size[0]);
-        parameters.parameter("num_dim").setValue(numDim);
-        String name = "grid_x";
-        for (int dim = 0; dim < numDim; dim++) {
-            final Matrix m = new Data(dim, denormalization);
-            parameters.parameter(name).setValue(m);
-            switch (dim) {
-                case 0: name = "grid_y"; break;         // Next parameter to set after "grid_x".
-                case 1: name = "grid_z"; break;         // Next parameter to set after "grid_y".
-                default: return;                        // Current implementation does not
show more than 3 dimensions.
-            }
-        }
+        parameters.parameter("grid_x").setValue(new Data(0, denormalization));
+        parameters.parameter("grid_y").setValue(new Data(1, denormalization));
     }
 
     /**
@@ -109,11 +102,6 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
     private final double[] offsets;
 
     /**
-     * Number of dimension of target coordinates.
-     */
-    private final int numDim;
-
-    /**
      * Conversion from grid coordinates to the final "real world" coordinates.
      *
      * @see #gridToTarget()
@@ -125,16 +113,14 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
      *
      * @param sourceToGrid  conversion from the "real world" source coordinates to grid indices
including fractional parts.
      * @param gridToTarget  conversion from grid coordinates to the final "real world" coordinates.
-     * @param numDim        number of dimension of target coordinates.
      * @param residuals     the residual data, as translations to apply on the result of
affine transform.
      * @param precision     desired precision of inverse transformations in unit of grid
cells.
      */
     ResidualGrid(final LinearTransform sourceToGrid, final LinearTransform gridToTarget,
-            final int nx, final int ny, final int numDim, final double[] residuals, final
double precision)
+            final int nx, final int ny, final double[] residuals, final double precision)
     {
         super(Units.UNITY, Units.UNITY, true, sourceToGrid, nx, ny, PARAMETERS);
         this.gridToTarget = gridToTarget;
-        this.numDim       = numDim;
         this.offsets      = residuals;
         this.accuracy     = precision;
     }
@@ -146,7 +132,6 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
     private ResidualGrid(final ResidualGrid other, final double[] data) {
         super(other);
         gridToTarget = other.gridToTarget;
-        numDim       = other.numDim;
         accuracy     = other.accuracy;
         offsets      = data;
     }
@@ -181,7 +166,7 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
      */
     @Override
     public int getTranslationDimensions() {
-        return numDim;
+        return SOURCE_DIMENSION;
     }
 
     /**
@@ -196,10 +181,11 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
 
     /**
      * 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) {
-        return offsets[(gridX + gridY*nx) * numDim + dim];
+        return offsets[(gridX + gridY*nx) * SOURCE_DIMENSION + dim];
     }
 
     /**
@@ -228,12 +214,9 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
 
         /** Computes the matrix element in the given row and column. */
         @Override public double  getElement(final int y, final int x) {
-            int i = affine.length;
-            double sum = affine[--i];
-            while (--i >= 0) {
-                sum += affine[i] * getCellValue(i, x, y);       // TODO: use Math.fma with
JDK9.
-            }
-            return sum;
+            return affine[0] * (x + getCellValue(0, x, y)) +                // TODO: use
Math.fma with JDK9.
+                   affine[1] * (y + getCellValue(1, x, y)) +
+                   affine[2];
         }
 
         /**
@@ -292,16 +275,4 @@ final class ResidualGrid extends DatumShiftGridFile<Dimensionless,Dimensionless>
             return "Matrix";
         }
     }
-
-    /**
-     * Returns {@code true} if the given object is a grid containing the same data than this
grid.
-     */
-    @Override
-    public boolean equals(final Object other) {
-        if (super.equals(other)) {
-            // Offset array has been compared by the parent class.
-            return numDim == ((ResidualGrid) other).numDim;
-        }
-        return false;
-    }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
index b15d4f0..9030ac0 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/referencing/operation/builder/ResidualGridTest.java
@@ -47,7 +47,7 @@ public final strictfp class ResidualGridTest extends TestCase {
      * The grid has no "source to grid" or "grid to CRS" transformations.
      */
     public ResidualGridTest() {
-        grid = new ResidualGrid(MathTransforms.identity(2), MathTransforms.identity(2), 3,
4, 2, new double[] {
+        grid = new ResidualGrid(MathTransforms.identity(2), MathTransforms.identity(2), 3,
4, new double[] {
                 0,2  ,  1,2  ,  2,1,
                 1,3  ,  2,2  ,  1,1,
                 0,4  ,  2,3  ,  3,2,


Mime
View raw message