sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] branch master updated: Apply the fix for NADCON grids crossing anti-meridian to NetCDF localization grids crossing the anti-meridian.
Date Thu, 20 Feb 2020 15:09:50 GMT
This is an automated email from the ASF dual-hosted git repository.

desruisseaux pushed a commit to branch master
in repository https://gitbox.apache.org/repos/asf/sis.git


The following commit(s) were added to refs/heads/master by this push:
     new 9de779a  Apply the fix for NADCON grids crossing anti-meridian to NetCDF localization grids crossing the anti-meridian.
9de779a is described below

commit 9de779a9108be4d4c739078b8a05e62eb1944256
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Thu Feb 20 16:08:33 2020 +0100

    Apply the fix for NADCON grids crossing anti-meridian to NetCDF localization grids crossing the anti-meridian.
---
 .../provider/DatumShiftGridCompressed.java         |  47 ++++----
 .../referencing/provider/DatumShiftGridFile.java   |  62 +++++------
 .../sis/referencing/datum/DatumShiftGrid.java      |  96 +++++++++-------
 .../operation/builder/LinearTransformBuilder.java  |   5 +
 .../operation/builder/LocalizationGridBuilder.java |  40 +++++--
 .../operation/builder/ResidualGrid.java            | 121 ++++++++++++++++++---
 .../operation/builder/package-info.java            |   2 +-
 .../sis/referencing/operation/builder/readme.html  |   6 +-
 .../operation/transform/LinearTransform.java       |   3 +-
 .../operation/transform/ProjectiveTransform.java   |  44 +++++---
 .../operation/builder/ResidualGridTest.java        |   8 +-
 11 files changed, 293 insertions(+), 141 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 139760e..0f6ec07 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
@@ -42,13 +42,7 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T
     /**
      * Serial number for inter-operability with different versions.
      */
-    private static final long serialVersionUID = 4847888093457104917L;
-
-    /**
-     * Maximal grid index along the <var>y</var> axis, inclusive.
-     * This is the number of grid cells minus 2.
-     */
-    private final int ymax;
+    private static final long serialVersionUID = 1889111858140209014L;
 
     /**
      * An "average" value for the offset in each dimension.
@@ -75,7 +69,6 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T
             final short[][] data, final double scale)
     {
         super(grid);
-        this.ymax     = getGridSize()[1] - 2;
         this.averages = averages;
         this.data     = data;
         this.scale    = scale;
@@ -181,7 +174,7 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T
      */
     @Override
     public double getCellValue(final int dim, final int gridX, final int gridY) {
-        return data[dim][gridX + gridY*nx] * scale + averages[dim];
+        return data[dim][gridX + gridY*scanlineStride] * scale + averages[dim];
     }
 
     /**
@@ -190,24 +183,26 @@ final class DatumShiftGridCompressed<C extends Quantity<C>, T extends Quantity<T
      */
     @Override
     public void interpolateInCell(double gridX, double gridY, double[] vector) {
-        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) {
-            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.
-            skipY |= (iy != ymax+1 || gridY != 0);      // Keep value 'false' if gridY == gridSize[1] - 1.
-            iy     = ymax;
-            gridY  = 1;
+        final int xmax = getGridSize(0) - 2;
+        final int ymax = getGridSize(1) - 2;
+        int ix = (int) gridX;                               // Really want rounding toward zero (not floor).
+        int iy = (int) gridY;
+        if (ix < 0 || ix > xmax || iy < 0 || iy > ymax) {
+            final double[] gridCoordinates = {gridX, gridY};
+            replaceOutsideGridCoordinates(gridCoordinates);
+            gridX = gridCoordinates[0];
+            gridY = gridCoordinates[1];
+            ix = Math.max(0, Math.min(xmax, (int) gridX));
+            iy = Math.max(0, Math.min(ymax, (int) gridY));
         }
-        final int p00 = nx*iy + ix;
-        final int p10 = nx + p00;
+        gridX -= ix;                                        // If was negative, will continue to be negative.
+        gridY -= iy;
+        boolean skipX = (gridX < 0); if (skipX) gridX = 0;
+        boolean skipY = (gridY < 0); if (skipY) gridY = 0;
+        if (gridX > 1)  {gridX = 1; skipX = true;}
+        if (gridY > 1)  {gridY = 1; skipY = true;}
+        final int p00 = scanlineStride * iy + ix;
+        final int p10 = scanlineStride + p00;
         final int n   = data.length;
         boolean derivative = (vector.length >= n + INTERPOLATED_DIMENSIONS * INTERPOLATED_DIMENSIONS);
         for (int dim = 0; dim < n; dim++) {
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
index 7ebb2d7..be13d1a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridFile.java
@@ -77,7 +77,7 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
     /**
      * Serial number for inter-operability with different versions.
      */
-    private static final long serialVersionUID = -4471670781277328193L;
+    private static final long serialVersionUID = -5801692909082130314L;
 
     /**
      * Cache of grids loaded so far. Those grids will be stored by soft references until the amount of
@@ -111,10 +111,11 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
     private final Path[] files;
 
     /**
-     * Number of grid cells along the <var>x</var> axis.
-     * This is <code>{@linkplain #getGridSize()}[0]</code> as a field for performance reasons.
+     * Number of cells between the start of adjacent rows in the grid. This is usually {@code getGridSize(0)},
+     * stored as a field for performance reasons. Value could be greater than {@code getGridSize(0)} if there
+     * is some elements to ignore at the end of each row.
      */
-    protected final int nx;
+    protected final int scanlineStride;
 
     /**
      * Number of cells that the grid would have if it was spanning 360° of longitude, or 0 if no wraparound
@@ -122,9 +123,9 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
      * an integer number of cells in 360°. This value is used for longitude values that are on the other side
      * of the ±180° meridian compared to the region where the grid is defined.
      *
-     * @see #replaceOutsideGridCoordinate(int, double)
+     * @see #replaceOutsideGridCoordinates(double[])
      */
-    private final double cycle;
+    private final double periodX;
 
     /**
      * The best translation accuracy that we can expect from this file.
@@ -182,13 +183,13 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
     {
         super(coordinateUnit, new AffineTransform2D(Δx, 0, 0, Δy, x0, y0).inverse(),
               new int[] {nx, ny}, isCellValueRatio, translationUnit);
-        this.descriptor = descriptor;
-        this.files      = files;
-        this.nx         = nx;
+        this.descriptor     = descriptor;
+        this.files          = files;
+        this.scanlineStride = nx;
         if (Units.isAngular(coordinateUnit)) {
-            cycle = Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / Math.abs(Δx));
+            periodX = Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / Math.abs(Δx));
         } else {
-            cycle = 0;
+            periodX = 0;
             /*
              * Note: non-angular source coordinates are currently never used in this package.
              * If it continue to be like that in the future, we should remove the check for
@@ -205,12 +206,12 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
      */
     protected DatumShiftGridFile(final DatumShiftGridFile<C,T> other) {
         super(other);
-        descriptor = other.descriptor;
-        files      = other.files;
-        nx         = other.nx;
-        accuracy   = other.accuracy;
-        subgrids   = other.subgrids;
-        cycle      = other.cycle;
+        descriptor     = other.descriptor;
+        files          = other.files;
+        scanlineStride = other.scanlineStride;
+        accuracy       = other.accuracy;
+        subgrids       = other.subgrids;
+        periodX        = other.periodX;
     }
 
     /**
@@ -230,11 +231,11 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
     {
         super(other.getCoordinateUnit(), gridToCRS.inverse(), new int[] {nx, ny},
               other.isCellValueRatio(), other.getTranslationUnit());
-        descriptor = other.descriptor;
-        files      = other.files;
-        this.nx    = nx;
-        cycle      = (other.cycle == 0) ? 0 :
-                Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE) / AffineTransforms2D.getScaleX0(gridToCRS));
+        scanlineStride = nx;
+        descriptor     = other.descriptor;
+        files          = other.files;
+        periodX        = (other.periodX == 0) ? 0 : Math.rint((Longitude.MAX_VALUE - Longitude.MIN_VALUE)
+                                                  / AffineTransforms2D.getScaleX0(gridToCRS));
         // Accuracy to be set by caller. Initial value needs to be zero.
     }
 
@@ -395,11 +396,10 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
      * that coordinate inside the range.
      */
     @Override
-    protected double replaceOutsideGridCoordinate(final int dimension, final double gridCoordinate) {
-        if (dimension == 0 && cycle != 0) {
-            return Math.IEEEremainder(gridCoordinate, cycle);
+    protected void replaceOutsideGridCoordinates(final double[] gridCoordinates) {
+        if (periodX != 0) {
+            gridCoordinates[0] = Math.IEEEremainder(gridCoordinates[0], periodX);
         }
-        return super.replaceOutsideGridCoordinate(dimension, gridCoordinate);
     }
 
     /**
@@ -495,7 +495,7 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
         /**
          * The translation values. {@code offsets.length} is the number of dimensions, and {@code offsets[dim].length}
          * shall be the same for all {@code dim} value. Component {@code dim} of the translation vector at coordinate
-         * {@code gridX}, {@code gridY} is {@code offsets[dim][gridX + gridY*nx]}.
+         * {@code gridX}, {@code gridY} is {@code offsets[dim][gridX + gridY*scanlineStride]}.
          */
         final float[][] offsets;
 
@@ -564,13 +564,13 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
          * from an ASCII file, or any other medium that format numbers in base 10.
          *
          * @param  dim    the dimension for which to get an average value.
-         * @param  gridX  the grid index along the <var>x</var> axis, from 0 inclusive to {@link #nx} exclusive.
-         * @param  gridY  the grid index along the <var>y</var> axis, from 0 inclusive to {@code  ny} exclusive.
+         * @param  gridX  the grid index along the <var>x</var> axis, from 0 inclusive to {@code nx} exclusive.
+         * @param  gridY  the grid index along the <var>y</var> axis, from 0 inclusive to {@code ny} exclusive.
          * @return the offset at the given dimension in the grid cell at the given index.
          */
         @Override
         public final double getCellValue(final int dim, final int gridX, final int gridY) {
-            return DecimalFunctions.floatToDouble(offsets[dim][gridX + gridY*nx]);
+            return DecimalFunctions.floatToDouble(offsets[dim][gridX + gridY*scanlineStride]);
         }
 
         /**
@@ -673,7 +673,7 @@ abstract class DatumShiftGridFile<C extends Quantity<C>, T extends Quantity<T>>
          */
         @Override
         public final double getCellValue(final int dim, final int gridX, final int gridY) {
-            return offsets[dim][gridX + gridY*nx];
+            return offsets[dim][gridX + gridY*scanlineStride];
         }
     }
 }
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 158540a..ab08860 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
@@ -124,7 +124,7 @@ import org.apache.sis.measure.Units;
  * in which case a longitude of 170° needs to be replaced by −190° before it can be processed by the grid.
  * The default {@code DatumShiftGrid} class does not apply longitude wraparound automatically
  * (it does not even know which axis, if any, is longitude),
- * but subclasses can add this support by overriding the {@link #replaceOutsideGridCoordinate(int, double)} method.
+ * but subclasses can add this support by overriding the {@link #replaceOutsideGridCoordinates(double[])} method.
  *
  * <h2>Sub-grids</h2>
  * Some datum shift grid files provide a grid valid on a wide region, refined with denser sub-grids in smaller regions.
@@ -283,6 +283,20 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
     }
 
     /**
+     * Returns the number of cells in the specified dimension.
+     * Invoking this method is equivalent to {@code getGridSize()[dimension]}.
+     *
+     * @param  dimension  the dimension for which to get the grid size.
+     * @return the number of grid cells in the specified dimension.
+     * @throws IndexOutOfBoundsException if the given index is out of bounds.
+     *
+     * @since 1.1
+     */
+    public int getGridSize(final int dimension) {
+        return gridSize[dimension];
+    }
+
+    /**
      * Returns the domain of validity of input coordinates that can be specified to the
      * {@link #interpolateAt interpolateAt(…)} method. Coordinates outside that domain
      * will still be accepted, but results may be extrapolations far from reality.
@@ -506,13 +520,13 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
         final int xmax = gridSize[0] - 2;
         final int ymax = gridSize[1] - 2;
         int ix = (int) gridX;                               // Really want rounding toward zero (not floor).
-        if (ix < 0 || ix > xmax) {
-            gridX = replaceOutsideGridCoordinate(0, gridX);
-            ix = Math.max(0, Math.min(xmax, (int) gridX));
-        }
         int iy = (int) gridY;
-        if (iy < 0 || iy > ymax) {
-            gridY = replaceOutsideGridCoordinate(1, gridY);
+        if (ix < 0 || ix > xmax || iy < 0 || iy > ymax) {
+            final double[] gridCoordinates = {gridX, gridY};
+            replaceOutsideGridCoordinates(gridCoordinates);
+            gridX = gridCoordinates[0];
+            gridY = gridCoordinates[1];
+            ix = Math.max(0, Math.min(xmax, (int) gridX));
             iy = Math.max(0, Math.min(ymax, (int) gridY));
         }
         gridX -= ix;                                        // If was negative, will continue to be negative.
@@ -608,13 +622,13 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
         final int xmax = gridSize[0] - 2;
         final int ymax = gridSize[1] - 2;
         int ix = (int) gridX;
-        if (ix < 0 || ix > xmax) {
-            gridX = replaceOutsideGridCoordinate(0, gridX);
-            ix = Math.max(0, Math.min(xmax, (int) gridX));
-        }
         int iy = (int) gridY;
-        if (iy < 0 || iy > ymax) {
-            gridY = replaceOutsideGridCoordinate(1, gridY);
+        if (ix < 0 || ix > xmax || iy < 0 || iy > ymax) {
+            final double[] gridCoordinates = {gridX, gridY};
+            replaceOutsideGridCoordinates(gridCoordinates);
+            gridX = gridCoordinates[0];
+            gridY = gridCoordinates[1];
+            ix = Math.max(0, Math.min(xmax, (int) gridX));
             iy = Math.max(0, Math.min(ymax, (int) gridY));
         }
         gridX -= ix;
@@ -745,17 +759,16 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
      * @since 1.0
      */
     public boolean isCellInGrid(double gridX, double gridY) {
-        double max = gridSize[1] - 1;
-        if (!(gridY >= 0 && gridY <= max)) {
-            gridY = replaceOutsideGridCoordinate(1, gridY);
-            if (!(gridY >= 0 && gridY <= max)) return false;
-        }
-        max = gridSize[0] - 1;
-        if (!(gridX >= 0 && gridX <= max)) {
-            gridX = replaceOutsideGridCoordinate(0, gridX);
-            if (!(gridX >= 0 && gridX <= max)) return false;
+        final double xmax = gridSize[0] - 1;
+        final double ymax = gridSize[1] - 1;
+        if (gridX >= 0 && gridX <= xmax && gridY >= 0 && gridY <= ymax) {
+            return true;
         }
-        return true;
+        final double[] gridCoordinates = {gridX, gridY};
+        replaceOutsideGridCoordinates(gridCoordinates);
+        gridX = gridCoordinates[0];
+        gridY = gridCoordinates[1];
+        return (gridX >= 0 && gridX <= xmax && gridY >= 0 && gridY <= ymax);
     }
 
     /**
@@ -765,35 +778,40 @@ public abstract class DatumShiftGrid<C extends Quantity<C>, T extends Quantity<T
      * longitude range, so a longitude of 170° needs to be converted to a longitude of −190° before it can be
      * processed by that grid.
      *
-     * <p>The default implementation returns {@code gridCoordinate} unchanged. Subclasses need to override this
-     * method if they want to handle longitude cycles. Note that the coordinate value is a grid index, not a
-     * longitude value. So the cycle to add or remove is the number of cells that the grid would have if it was
-     * spanning 360° of longitude.</p>
+     * <p>The default implementation does nothing. Subclasses need to override this method if they want to handle
+     * longitude wraparounds. Note that the coordinate values are grid indices, not longitude or latitude values.
+     * So the period to add or remove is the number of cells that the grid would have if it was spanning 360° of
+     * longitude.</p>
      *
      * <div class="note"><b>Example:</b>
-     * this method may be implemented as below:
+     * if longitude values are mapped to {@code gridX} coordinates (in dimension 0), and if a shift of 360° in
+     * longitude values is equivalent to a shift of {@code periodX} cells in the grid, then this method can be
+     * implemented as below:
      *
      * {@preformat java
+     *     private final double periodX = ...;      // Number of grid cells in 360° of longitude.
+     *
      *     &#64;Override
-     *     protected double replaceOutsideGridCoordinate(int dimension, double gridCoordinate) {
-     *         if (dimension == 0) {
-     *             return Math.IEEEremainder(gridCoordinate, cycle);
-     *         } else {
-     *             return super.replaceOutsideGridCoordinate(dimension, gridCoordinate);
-     *         }
+     *     protected void replaceOutsideGridCoordinates(double[] gridCoordinates) {
+     *         gridCoordinates[0] = Math.IEEEremainder(gridCoordinates[0], periodX);
      *     }
      * }</div>
      *
-     * @param  dimension       0 if the given coordinate is {@code gridX}, or 1 if it is {@code gridY}.
-     * @param  gridCoordinate  the grid coordinate which is outside the range of valid values.
-     * @return the grid coordinate to use. May still be outside the range of valid grid values.
+     * This method receives all grid coordinates in the {@code gridCoordinates} argument and can modify any
+     * of them, possibly many at once. The reason is because a shift of 360° of longitude (for example) may
+     * correspond to an offset in both {@code gridX} and {@code gridY} indices if the grid is inclined.
+     * The {@code gridCoordinates} array length is the number of grid dimensions,
+     * typically {@value #INTERPOLATED_DIMENSIONS}.
+     *
+     * @param  gridCoordinates  on input, the cell indices of the point which is outside the grid.
+     *         On output, the cell indices of an equivalent point inside the grid if possible.
+     *         Coordinate values are modified in-place.
      *
      * @see Math#IEEEremainder(double, double)
      *
      * @since 1.1
      */
-    protected double replaceOutsideGridCoordinate(final int dimension, final double gridCoordinate) {
-        return gridCoordinate;
+    protected void replaceOutsideGridCoordinates(double[] gridCoordinates) {
     }
 
     /**
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 412f7d7..d5849ab 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
@@ -1324,6 +1324,11 @@ search:         for (int j=domain(); --j >= 0;) {
      * Computes the matrix of the linear approximation. This is the implementation of {@link #create(MathTransformFactory)}
      * without the step creating the {@link LinearTransform} from a matrix. The {@link #correlations} field is set as a side
      * effect of this method call.
+     *
+     * <p>In current implementation, the transform represented by the returned matrix is always affine
+     * (i.e. the last row is fixed to [0 0 … 0 1]). If this is no longer the case in a future version,
+     * some codes may not work anymore. Search for {@code isAffine()} statements for locating codes
+     * that depend on affine transform assumption.</p>
      */
     @SuppressWarnings("serial")
     private MatrixSIS fit() throws FactoryException {
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 8952f30..c8046ab 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
@@ -84,7 +84,7 @@ import static org.apache.sis.referencing.operation.builder.ResidualGrid.SOURCE_D
  * See the <cite>Linearizers</cite> section in {@link LinearTransformBuilder} for more discussion.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see InterpolatedTransform
  * @see LinearTransform
@@ -139,6 +139,17 @@ public class LocalizationGridBuilder extends TransformBuilder {
     private MathTransform transform;
 
     /**
+     * If the coordinates in some dimensions are cyclic, their periods. Otherwise {@code null}.
+     * Values are in units of the target CRS. For longitude wraparounds, the period is typically 360°.
+     * Array length shall be {@code linear.getTargetDimensions()} and non-cyclic dimensions shall have
+     * a period of zero (not {@link Double#NaN}, because we will use this array as a displacement vector).
+     *
+     * @see #resolveWraparoundAxis(int, int, double)
+     * @see ResidualGrid#cellPeriods
+     */
+    private double[] periods;
+
+    /**
      * Creates a new, initially empty, builder for a localization grid of the given size.
      *
      * @param width   the number of columns in the grid of target positions.
@@ -549,7 +560,8 @@ public class LocalizationGridBuilder extends TransformBuilder {
      *                    The recommended direction is the direction of most stable values, typically 1 (rows) for longitudes.
      * @param  period     that wraparound range (typically 360° for longitudes). Must be strictly positive.
      * @return the range of coordinate values in the specified dimension after correction for wraparound values.
-     * @throws IllegalStateException if {@link #create(MathTransformFactory) create(…)} has already been invoked.
+     * @throws IllegalStateException if this method has already been invoked for the same dimension,
+     *         or if {@link #create(MathTransformFactory) create(…)} has already been invoked.
      *
      * @since 1.0
      */
@@ -558,6 +570,14 @@ public class LocalizationGridBuilder extends TransformBuilder {
         ArgumentChecks.ensureBetween("dimension", 0, linear.getTargetDimensions() - 1, dimension);
         ArgumentChecks.ensureBetween("direction", 0, linear.getSourceDimensions() - 1, direction);
         ArgumentChecks.ensureStrictlyPositive("period", period);
+        if (periods == null) {
+            periods = new double[linear.getTargetDimensions()];
+        }
+        if (periods[dimension] != 0) {
+            throw new IllegalStateException(Errors.format(
+                    Errors.Keys.ValueAlreadyDefined_1, Strings.bracket("periods", dimension)));
+        }
+        periods[dimension] = period;
         return linear.resolveWraparoundAxis(dimension, direction, period);
     }
 
@@ -670,16 +690,16 @@ public class LocalizationGridBuilder extends TransformBuilder {
                             residual[k++] = (float) dy;
                         }
                     }
+                    if (isLinear) {
+                        step = MathTransforms.concatenate(sourceToGrid, gridToCoord);
+                    } else {
+                        step = InterpolatedTransform.createGeodeticTransformation(nonNull(factory),
+                                new ResidualGrid(sourceToGrid, gridToCoord, width, height, residual,
+                                (gridPrecision > 0) ? gridPrecision : DEFAULT_PRECISION, periods));
+                    }
                 } catch (TransformException e) {
                     throw new FactoryException(e);                                          // Should never happen.
                 }
-                if (isLinear) {
-                    step = MathTransforms.concatenate(sourceToGrid, gridToCoord);
-                } else {
-                    step = InterpolatedTransform.createGeodeticTransformation(nonNull(factory),
-                            new ResidualGrid(sourceToGrid, gridToCoord, width, height, residual,
-                            (gridPrecision > 0) ? gridPrecision : DEFAULT_PRECISION));
-                }
             }
             /*
              * At this point we finished to compute the transformation to target coordinates.
@@ -693,7 +713,7 @@ public class LocalizationGridBuilder extends TransformBuilder {
                 throw new InvalidGeodeticParameterException(Resources.format(
                         Resources.Keys.NonInvertibleOperation_1, linear.linearizerID()), e);
             }
-            transform = step;
+            transform = step;                               // Set only after everything succeeded.
         }
         return transform;
     }
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 46565c3..4776b98 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
@@ -22,6 +22,7 @@ import javax.measure.quantity.Dimensionless;
 import org.opengis.parameter.ParameterDescriptor;
 import org.opengis.parameter.ParameterDescriptorGroup;
 import org.opengis.referencing.operation.Matrix;
+import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.parameter.Parameters;
 import org.apache.sis.parameter.ParameterBuilder;
 import org.apache.sis.referencing.datum.DatumShiftGrid;
@@ -42,7 +43,7 @@ import org.apache.sis.measure.Units;
  * The residuals after an affine approximation has been created for a set of matching control point pairs.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.8
  * @module
  */
@@ -50,7 +51,7 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless> {
     /**
      * For cross-version compatibility.
      */
-    private static final long serialVersionUID = 5207799661806374259L;
+    private static final long serialVersionUID = -3668228260650927123L;
 
     /**
      * Number of source dimensions of the residual grid.
@@ -103,10 +104,11 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless> {
     }
 
     /**
-     * Number of grid cells along the <var>x</var> axis.
-     * This is <code>{@linkplain #getGridSize()}[0]</code> as a field for performance reasons.
+     * Number of cells between the start of adjacent rows in the grid. This is usually {@code getGridSize(0)},
+     * stored as a field for performance reasons. Value could be greater than {@code getGridSize(0)} if there
+     * is some elements to ignore at the end of each row.
      */
-    private final int nx;
+    private final int scanlineStride;
 
     /**
      * The residual data, as translations to apply on the result of affine transform.
@@ -133,21 +135,70 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless> {
     private final double accuracy;
 
     /**
+     * If grid coordinates in some target dimensions are cyclic, the period in number of cells.
+     * For each scalar value in the {@link LocalizationGridBuilder#periods} array (in units of
+     * target CRS), the corresponding period in number of cells is a vector. For example a 360°
+     * shift in longitude does not necessarily correspond to an horizontal or vertical offset
+     * in grid indices; it may be a combination of both if the grid is inclined.
+     *
+     * <p>We should have as many vectors as non-zero values in {@link LocalizationGridBuilder#periods}.
+     * Each {@code periodVector} (in cell units) should be computed from a {@code periods} vector with
+     * exactly one non-zero value (in CRS units) for allowing shifts in different CRS dimensions to be
+     * applied independently. Consequently this field should actually be of type {@code double[][]}.
+     * But current version uses only one vector for avoiding the complexity of searching how to combine
+     * multiple vectors. It is okay for the usual case where only one CRS axis has wraparound range,
+     * but may need to be revisited in the future.</p>
+     *
+     * <p>This array is {@code null} if no period has been specified, or if a period has been specified
+     * but we can not convert it from CRS units to a constant number of cells.</p>
+     *
+     * @see LocalizationGridBuilder#periods
+     * @see #replaceOutsideGridCoordinates(double[])
+     */
+    private final double[] periodVector;
+
+    /**
      * Creates a new residual grid.
      *
      * @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 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.
+     * @param periods       if grid coordinates in some dimensions are cyclic, their periods in units of target CRS.
      */
     ResidualGrid(final LinearTransform sourceToGrid, final LinearTransform gridToTarget,
-            final int nx, final int ny, final float[] residuals, final double precision)
+            final int nx, final int ny, final float[] residuals, final double precision,
+            final double[] periods) throws TransformException
     {
         super(Units.UNITY, sourceToGrid, new int[] {nx, ny}, true, Units.UNITY);
-        this.gridToTarget = gridToTarget;
-        this.offsets      = residuals;
-        this.accuracy     = precision;
-        this.nx           = nx;
+        this.gridToTarget   = gridToTarget;
+        this.offsets        = residuals;
+        this.accuracy       = precision;
+        this.scanlineStride = nx;
+        double[] periodVector = null;
+        if (periods != null && gridToTarget.isAffine()) {
+            /*
+             * We require the transform to be affine because it makes the Jacobian independent of
+             * coordinate values. It allows us to replace a period in target CRS units by periods
+             * in grid units without having to take the coordinate values in account.
+             */
+            final MatrixSIS m = MatrixSIS.castOrCopy(gridToTarget.inverse().derivative(null));
+            periodVector = m.multiply(periods);
+            /*
+             * Above code is not really right. We should verify that `periods` contains exactly
+             * one non-zero element, and if not the case execute above code in a loop with one
+             * non-zero element by iteration, creating one independent vector each time.
+             * We don't do that for now because `replaceOutsideGridCoordinates(…)` is not yet
+             * capable to search the best combination of many vectors.
+             *
+             * With current implementation, if the `periods` array contains more than 1 non-zero value,
+             * then `replaceOutsideGridCoordinates(…)` always shift all wraparound dimensions together.
+             * For example if the first dimension has a period of 360° and the second dimension has a
+             * period of 12 months, then `replaceOutsideGridCoordinates(…)` will only shift by 360°
+             * AND 12 months together, never 360° only or 12 months only.
+             */
+        }
+        this.periodVector = periodVector;
     }
 
     /**
@@ -192,7 +243,48 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless> {
      */
     @Override
     public double getCellValue(int dim, int gridX, int gridY) {
-        return offsets[(gridX + gridY*nx) * SOURCE_DIMENSION + dim];
+        return offsets[(gridX + gridY*scanlineStride) * SOURCE_DIMENSION + dim];
+    }
+
+    /**
+     * Invoked when a {@code gridX} or {@code gridY} coordinate is outside the range of valid grid coordinates.
+     * If the coordinate outside the range is a longitude value and if we handle those values as cyclic, brings
+     * that coordinate inside the range.
+     */
+    @Override
+    protected void replaceOutsideGridCoordinates(final double[] gridCoordinates) {
+        if (periodVector != null) {
+            /*
+             * We will try to shift the point by an integral multiple of `periodVector`.
+             * Test along each dimension which range of multiples could bring the point
+             * inside the grid, and take the range intersection in all dimensions.
+             */
+            double min = Double.NEGATIVE_INFINITY;
+            double max = Double.POSITIVE_INFINITY;
+            for (int i=0; i<gridCoordinates.length; i++) {
+                final double period = periodVector[i];
+                double toLower = gridCoordinates[i];                // If we subtract a shift ≤ toLower, then coordinate ≥ lower bound (which is 0).
+                double toUpper = toLower - (getGridSize(i) - 1);    // If we subtract a shift ≥ toUpper, then coordinate ≤ upper bound (which is gridSize−1).
+                toLower = Math.floor(toLower / period);             // Convert to integral number of periods, reverse may be ≤ original value.
+                toUpper = Math.ceil (toUpper / period);             // Convert to integral number of periods, reverse may be ≥ original value.
+                if (toLower < max) max = toLower;                   // Must shift by no more than this value otherwise coordinate < lower bound.
+                if (toUpper > min) min = toUpper;                   // Must shift by no less than this value otherwise coordinate > upper bound.
+            }
+            if (min <= max) {
+                /*
+                 * If at least one multiple exists, take the multiple closest to zero
+                 * (i.e. apply the smallest displacement). Note that the range should
+                 * not include zero, otherwise this point would be inside the grid and
+                 * this method should not have been invoked.
+                 */
+                final double n = (min >= 0) ? min : max;
+                if (Double.isFinite(n)) {
+                    for (int i=0; i<gridCoordinates.length; i++) {
+                        gridCoordinates[i] -= periodVector[i] * n;
+                    }
+                }
+            }
+        }
     }
 
     /**
@@ -224,8 +316,8 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless> {
         @SuppressWarnings("CloneInNonCloneableClass")
         @Override public Matrix  clone()                            {return this;}
         @Override public boolean isIdentity()                       {return false;}
-        @Override public int     getNumCol()                        {return nx;}
-        @Override public int     getNumRow()                        {return getGridSize()[1];}
+        @Override public int     getNumCol()                        {return getGridSize(0);}
+        @Override public int     getNumRow()                        {return getGridSize(1);}
         @Override public Number  apply     (int[] p)                {return getElement(p[1], p[0]);}
         @Override public void    setElement(int y, int x, double v) {throw new UnsupportedOperationException();}
 
@@ -241,10 +333,9 @@ final class ResidualGrid extends DatumShiftGrid<Dimensionless,Dimensionless> {
          * in the table formatted for {@link ParameterDescriptorGroup} string representation.
          */
         @Override public String toString() {
-            final int[] size = getGridSize();
             return new StringBuilder(80).append('[')
                     .append(getElement(0, 0)).append(", …, ")
-                    .append(getElement(size[1] - 1, size[0] - 1))
+                    .append(getElement(getGridSize(1) - 1, getGridSize(0) - 1))
                     .append(']').toString();
         }
 
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
index 5c0774b..2ffb543 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/package-info.java
@@ -30,7 +30,7 @@
  * convenience.</p>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   0.5
  * @module
  */
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html
index bef92a0..6cd3ec1 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/readme.html
@@ -21,15 +21,17 @@
       typical linear regression method assumes that <var>x</var> values are exact and all errors are in <var>y</var> values.
       Applied to the <code>LocalizationGridBuilder</code> context, it means that linear regression method
       assumes that grid indices are exact and all errors are in geospatial coordinates.
-      This is a reasonable assumption if the linear regression is used directly.
+    </p><p>
+      The assumption that all errors are in geospatial coordinates is reasonable if the linear regression is used directly.
       But in <code>LocalizationGridBuilder</code> context, having the smallest errors on geospatial coordinates
       may not be so important because those errors are corrected by the residual grids during <em>forward</em> transformations.
       However during <em>inverse</em> transformations, it may be useful that grid indices estimated by the linear regression
       are as close as possible to the real grid indices in order to allow iterations to converge faster
       (such iterations exist only in inverse operations, not in forward operations).
       For that reason, <code>LocalizationGridBuilder</code> may want to minimize errors on grid indices instead than geospatial coordinates.
-      We can achieve this result by computing linear regression coefficients for an equation of the form
+      We could achieve this result by computing linear regression coefficients for an equation of the form
       <var>x</var> = <var>C₀</var> + <var>C₁</var> × <var>y</var>, then inverting that equation.
+    </p><p>
       This approach was attempted in a previous version and gave encouraging results, but we nevertheless reverted that change.
       One reason is that even if inverting the linear regression calculation allowed iterations to converge more often with curved
       localization grids, it was not sufficient anyway: we still had too many cases where inverse transformations did not converge.
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java
index 9cac821..3776331 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/LinearTransform.java
@@ -75,7 +75,8 @@ import org.opengis.referencing.operation.NoninvertibleTransformException;
 public interface LinearTransform extends MathTransform {
     /**
      * Returns {@code true} if this transform is affine.
-     * An affine transform preserves parallelism.
+     * An affine transform preserves parallelism and has the same
+     * {@linkplain AbstractMathTransform#derivative derivative} at every locations.
      *
      * @return {@code true} if this transform is affine.
      *
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
index a4efca3..7516f1f 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/transform/ProjectiveTransform.java
@@ -34,7 +34,7 @@ import org.apache.sis.util.ArgumentChecks;
  * lines in the source is preserved in the output.</p>
  *
  * @author  Martin Desruisseaux (IRD, Geomatys)
- * @version 0.7
+ * @version 1.1
  *
  * @see java.awt.geom.AffineTransform
  *
@@ -263,7 +263,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
         while (--numPts >= 0) {
             int mix = 0;
             for (int j=0; j<numRow; j++) {
-                double sum = elt[mix + srcDim];
+                double sum = elt[mix + srcDim];         // Initialize to translation term.
                 for (int i=0; i<srcDim; i++) {
                     final double e = elt[mix++];
                     if (e != 0) {
@@ -282,7 +282,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
             }
             final double w = buffer[dstDim];
             for (int j=0; j<dstDim; j++) {
-                // 'w' is equal to 1 if the transform is affine.
+                // `w` is equal to 1 if the transform is affine.
                 dstPts[dstOff + j] = buffer[j] / w;
             }
             srcOff += srcInc;
@@ -309,8 +309,8 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
     @Override
     public final void transform(float[] srcPts, int srcOff, final float[] dstPts, int dstOff, int numPts) {
         final int srcDim, dstDim;
-        int srcInc = srcDim = numCol-1;
-        int dstInc = dstDim = numRow-1;
+        int srcInc = srcDim = numCol - 1;
+        int dstInc = dstDim = numRow - 1;
         if (srcPts == dstPts) {
             switch (IterationStrategy.suggest(srcOff, srcDim, dstOff, dstDim, numPts)) {
                 case ASCENDING: {
@@ -337,7 +337,7 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
                 double sum = elt[mix + srcDim];
                 for (int i=0; i<srcDim; i++) {
                     final double e = elt[mix++];
-                    if (e != 0) { // See comment in transform(double[], ...)
+                    if (e != 0) {                                   // See comment in transform(double[], ...)
                         sum += srcPts[srcOff + i] * e;
                     }
                 }
@@ -364,8 +364,8 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
      */
     @Override
     public final void transform(final double[] srcPts, int srcOff, final float[] dstPts, int dstOff, int numPts) {
-        final int srcDim = numCol-1;
-        final int dstDim = numRow-1;
+        final int srcDim = numCol - 1;
+        final int dstDim = numRow - 1;
         final double[] buffer = new double[numRow];
         while (--numPts >= 0) {
             int mix = 0;
@@ -424,20 +424,38 @@ class ProjectiveTransform extends AbstractLinearTransform implements ExtendedPre
     }
 
     /**
-     * Gets the derivative of this transform at a point.
-     * For a matrix transform, the derivative is the same everywhere.
+     * Gets the derivative of this transform at a point. For an affine transform,
+     * the derivative is the same everywhere and the given point is ignored.
+     * For a perspective transform, the given point is used.
      *
-     * @param  point  ignored (can be {@code null}).
+     * @param  point  the coordinate point where to evaluate the derivative.
      */
     @Override
     public final Matrix derivative(final DirectPosition point) {
         final int srcDim = numCol - 1;
         final int dstDim = numRow - 1;
+        /*
+         * In the `transform(…)` method, all coordinate values are divided by a `w` coefficient
+         * which depends on the position. We need to reproduce that division here. Note that `w`
+         * coefficient is different than 1 only if the transform is non-affine.
+         */
+        int mix = dstDim * numCol;
+        double w = elt[mix + srcDim];                   // `w` is equal to 1 if the transform is affine.
+        for (int i=0; i<srcDim; i++) {
+            final double e = elt[mix++];
+            if (e != 0) {                               // For avoiding NullPointerException if affine.
+                w += point.getOrdinate(i) * e;
+            }
+        }
+        /*
+         * In the usual affine case (w=1), the derivative is a copy of this matrix
+         * with last row and last column omitted.
+         */
+        mix = 0;
         final MatrixSIS matrix = Matrices.createZero(dstDim, srcDim);
-        int mix = 0;
         for (int j=0; j<dstDim; j++) {
             for (int i=0; i<srcDim; i++) {
-                matrix.setElement(j, i, elt[mix++]);
+                matrix.setElement(j, i, elt[mix++] / w);
             }
             mix++;                                      // Skip translation column.
         }
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 800e646..a989968 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
@@ -31,7 +31,7 @@ import static org.apache.sis.test.Assert.*;
  * Tests {@link ResidualGrid}.
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 1.0
+ * @version 1.1
  * @since   1.0
  * @module
  */
@@ -45,13 +45,15 @@ public final strictfp class ResidualGridTest extends TestCase {
      * Creates a new test case with a 3×4 grid with 2 values in each cells.
      * Those two values are typically the horizontal components of translation vectors.
      * The grid has no "source to grid" or "grid to CRS" transformations.
+     *
+     * @throws TransformException if an error occurred while handling a wraparound axis.
      */
-    public ResidualGridTest() {
+    public ResidualGridTest() throws TransformException {
         grid = new ResidualGrid(MathTransforms.identity(2), MathTransforms.identity(2), 3, 4, new float[] {
                 0,2  ,  1,2  ,  2,1,
                 1,3  ,  2,2  ,  1,1,
                 0,4  ,  2,3  ,  3,2,
-                1,4  ,  3,3  ,  3,2}, 0.1);
+                1,4  ,  3,3  ,  3,2}, 0.1, null);
     }
 
     /**


Mime
View raw message