sis-commits mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From desruisse...@apache.org
Subject [sis] 03/05: When a point to transform is outside all sub-grids provided by a NTv2 file, take the nearest grid.
Date Sat, 15 Feb 2020 17:13:01 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 c484cfb720ec4c5d683f44e94a642764515f0f98
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Sat Feb 15 15:31:47 2020 +0100

    When a point to transform is outside all sub-grids provided by a NTv2 file, take the nearest
grid.
---
 .../referencing/j2d/IntervalRectangle.java         | 48 ++++++++++++++++-
 .../referencing/provider/DatumShiftGridGroup.java  | 63 +++++++++++-----------
 .../internal/referencing/provider/NTv2Test.java    |  4 +-
 3 files changed, 82 insertions(+), 33 deletions(-)

diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java
index a303887..fb589f4 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/j2d/IntervalRectangle.java
@@ -41,6 +41,12 @@ import org.apache.sis.internal.util.Strings;
  *       {@link Envelope2D#contains(Rectangle2D)} or {@link Envelope2D#intersects(Rectangle2D)}
methods.</li>
  * </ul>
  *
+ * This class provides the following additional methods which are not defined in {@link Rectangle2D}:
+ * <ul>
+ *   <li>{@link #containsInclusive(double, double)}</li>
+ *   <li>{@link #distance(double, double)}</li>
+ * </ul>
+ *
  * This class does <strong>not</strong> support by itself rectangles crossing
the anti-meridian of a geographic CRS.
  * However the {@link #getX()}, {@link #getY()}, {@link #getWidth()} and {@link #getHeight()}
methods are defined in
  * the straightforward way expected by {@link Envelope2D#intersects(Rectangle2D)} and similar
methods for computing
@@ -53,7 +59,7 @@ import org.apache.sis.internal.util.Strings;
  * recommended approach, but for Apache SIS private classes this is a way to reduce pressure
on garbage collector.</div>
  *
  * @author  Martin Desruisseaux (Geomatys)
- * @version 0.8
+ * @version 1.1
  * @since   0.8
  * @module
  */
@@ -370,6 +376,8 @@ public class IntervalRectangle extends Rectangle2D {
 
     /**
      * Tests if a specified coordinate is inside the boundary of this {@code Rectangle2D}.
+     * The maximal <var>x</var> and <var>y</var> values (i.e. the
right and bottom edges
+     * of this rectangle) are exclusive.
      *
      * @param  x  the <var>x</var> coordinates to test.
      * @param  y  the <var>y</var> coordinates to test.
@@ -381,6 +389,44 @@ public class IntervalRectangle extends Rectangle2D {
     }
 
     /**
+     * Tests if a specified coordinate is inside the boundary of this {@code Rectangle2D}
+     * with all edges inclusive.
+     *
+     * @param  x  the <var>x</var> coordinates to test.
+     * @param  y  the <var>y</var> coordinates to test.
+     * @return {@code true} if the specified coordinates are inside this rectangle, all edges
included.
+     */
+    public final boolean containsInclusive(final double x, final double y) {
+        return (x >= xmin && x <= xmax && y >= ymin && y
<= ymax);
+    }
+
+    /**
+     * Returns the minimal distance between a point and this rectangle.
+     * If the point is inside the rectangle or on the edge, then this method returns 0.
+     *
+     * @param  x  the <var>x</var> coordinates to test.
+     * @param  y  the <var>y</var> coordinates to test.
+     * @return minimal distance, or 0 if the point is inside this rectangle.
+     */
+    public final double distance(final double x, final double y) {
+        int outcode = 0;
+        double dx = java.lang.Double.POSITIVE_INFINITY;
+        double dy = java.lang.Double.POSITIVE_INFINITY;
+        double d;
+        if ((d = (x - xmax)) >= 0)           {dx = d; outcode =  1;}
+        if ((d = (y - ymax)) >= 0)           {dy = d; outcode |= 2;}
+        if ((d = (xmin - x)) >= 0 && d < dx) {dx = d; outcode |= 1;}
+        if ((d = (ymin - y)) >= 0 && d < dy) {dy = d; outcode |= 2;}
+        switch (outcode) {
+            case 1:  return dx;                                     // Only x coordinate
is outside.
+            case 2:  return dy;                                     // Only y coordinate
is outside.
+            case 3:  return Math.hypot(dx, dy);                     // Both coordinates are
outside.
+            case 0:  assert containsInclusive(x, y); return 0;      // No coordinate is outside.
+            default: throw new AssertionError(outcode);
+        }
+    }
+
+    /**
      * Determines where the specified coordinates lie with respect to this rectangle.
      * This method computes a binary OR of the appropriate mask values indicating,
      * for each side of this {@code Rectangle2D}, whether or not the specified coordinates
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
index 966bbbf..a009a92 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/provider/DatumShiftGridGroup.java
@@ -29,6 +29,7 @@ import org.opengis.util.FactoryException;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
 import org.apache.sis.referencing.operation.transform.InterpolatedTransform;
 import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
+import org.apache.sis.internal.referencing.j2d.IntervalRectangle;
 import org.apache.sis.internal.referencing.j2d.TileOrganizer;
 import org.apache.sis.internal.referencing.j2d.Tile;
 import org.apache.sis.internal.referencing.Resources;
@@ -65,16 +66,13 @@ final class DatumShiftGridGroup<C extends Quantity<C>, T extends
Quantity<T>> ex
      * All values in this class are integers, but nevertheless stored as {@code double} for
avoiding to cast them every
      * time {@link DatumShiftGridGroup#interpolateInCell(double, double, double[])} is executed.
      */
-    private static final class Region {
-        /** Grid bounds in units of the grid having finest resolution. */
-        private final double xmin, xmax, ymin, ymax;
-
+    private static final class Region extends IntervalRectangle {
         /** Subsampling compared to the grid having finest resolution. */
         private final double sx, sy;
 
         /** Creates a new instance from the given {@link TileOrganizer} result. */
         Region(final Tile tile) throws IOException {
-            final Rectangle r = tile.getAbsoluteRegion();
+            final Rectangle r = tile.getAbsoluteRegion();       // In units of the grid having
finest resolution.
             final Dimension s = tile.getSubsampling();
             xmin = r.getMinX();
             xmax = r.getMaxX();
@@ -84,11 +82,6 @@ final class DatumShiftGridGroup<C extends Quantity<C>, T extends
Quantity<T>> ex
             sy   = s.height;
         }
 
-        /** Tests whether the given coordinates are included in this region. */
-        final boolean contains(final double x, final double y) {
-            return x >= xmin && x <= xmax && y >= ymin &&
y <= ymax;
-        }
-
         /** Converts a coordinate from the parent grid to this grid. */
         final double x(final double p) {return (p - xmin) / sx;}
         final double y(final double p) {return (p - ymin) / sy;}
@@ -234,6 +227,9 @@ final class DatumShiftGridGroup<C extends Quantity<C>, T extends
Quantity<T>> ex
      * but should never be invoked. The {@link InterpolatedTransform} class will rather invoke
the
      * {@code interpolateInCell(…)} method for efficiency.
      *
+     * <p>Caller must ensure that all arguments given to this method are in their expected
ranges.
+     * The behavior of this method is undefined if any argument value is out-of-range.</p>
+     *
      * @param  dim    the dimension of the translation vector component to get.
      * @param  gridX  the grid index on the <var>x</var> axis, from 0 inclusive
to {@code gridSize[0]} exclusive.
      * @param  gridY  the grid index on the <var>y</var> axis, from 0 inclusive
to {@code gridSize[1]} exclusive.
@@ -243,7 +239,7 @@ final class DatumShiftGridGroup<C extends Quantity<C>, T extends
Quantity<T>> ex
     public double getCellValue(final int dim, final int gridX, final int gridY) {
         for (int i=0; i<regions.length; i++) {
             final Region r = regions[i];
-            if (r.contains(gridX, gridY)) {
+            if (r.containsInclusive(gridX, gridY)) {
                 double shift = subgrids[i].getCellValue(dim,
                         Math.toIntExact(Math.round(r.x(gridX))),
                         Math.toIntExact(Math.round(r.y(gridY))));
@@ -258,14 +254,20 @@ final class DatumShiftGridGroup<C extends Quantity<C>, T extends
Quantity<T>> ex
                 return shift;
             }
         }
-        throw new IndexOutOfBoundsException();
+        /*
+         * May be in the valid range of this DatumShiftGridGroup but not in the range of
a subgrid.
+         * This situation may happen if there is holes in the data coverage provided by subgrids.
+         */
+        throw new IllegalArgumentException();
     }
 
     /**
-     * Interpolates the translation to apply for the given two-dimensional grid indices.
The result is stored
-     * in the given {@code vector} array. This method is invoked only as a fallback if the
transform has not
-     * been able to use directly one of the child transforms. Consequently this implementation
does not need
-     * to be very fast.
+     * Interpolates the translation to apply for the given two-dimensional grid indices.
+     * This method is a fallback used when {@code SpecializableTransform} has not been able
to use directly
+     * one of the child transforms — it should happen only in a minority of cases, so performance
is not the
+     * priority here. The given point is assumed outside all sub-grids (otherwise {@code
SpecializableTransform}
+     * would not be invoking this fallback). Consequently searching a sub-grid containing
the given point is not
+     * sufficient; we have to search for the nearest grid even if the point is outside.
      *
      * @param  gridX   first grid coordinate of the point for which to get the translation.
      * @param  gridY   second grid coordinate of the point for which to get the translation.
@@ -273,22 +275,23 @@ final class DatumShiftGridGroup<C extends Quantity<C>, T extends
Quantity<T>> ex
      */
     @Override
     public void interpolateInCell(final double gridX, final double gridY, final double[]
vector) {
-        for (int i=0; i<regions.length; i++) {
+        int ni = 0;
+        Region nearest = regions[ni];
+        double distance = nearest.distance(gridX, gridY);
+        for (int i=1; i<regions.length; i++) {
             final Region r = regions[i];
-            if (r.contains(gridX, gridY)) {
-                subgrids[i].interpolateInCell(r.x(gridX), r.y(gridY), vector);
-                if (isCellValueRatio()) {
-                    for (int dim=0; dim < INTERPOLATED_DIMENSIONS; dim++) {
-                        vector[dim] *= r.relativeCellSize(dim);
-                    }
-                }
-                return;
+            final double d = r.distance(gridX, gridY);
+            if (d < distance) {
+                distance = d;
+                nearest  = r;
+                ni       = i;
+            }
+        }
+        subgrids[ni].interpolateInCell(nearest.x(gridX), nearest.y(gridY), vector);
+        if (isCellValueRatio()) {
+            for (int dim=0; dim < INTERPOLATED_DIMENSIONS; dim++) {
+                vector[dim] *= nearest.relativeCellSize(dim);
             }
         }
-        /*
-         * The following method call will (indirectly) invokes the above `getCellValue(…)`
method.
-         * It can be used as a way to test that method.
-         */
-        super.interpolateInCell(gridX, gridY, vector);
     }
 }
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java
b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java
index 3efe47f..27581e4 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/provider/NTv2Test.java
@@ -224,7 +224,7 @@ public final strictfp class NTv2Test extends DatumShiftTestCase {
          * should not be invoked in normal usage, so a direct invocation is the only way
to test it.
          */
         final double[] result = new double[] {
-            position[0] - grid.getCellValue(0, gridX, gridY) * cellSize,
+            position[0] - grid.getCellValue(0, gridX, gridY) * cellSize,    // Positive translation
is toward west.
             position[1] + grid.getCellValue(1, gridX, gridY) * cellSize
         };
         assertArrayEquals("getCellValue", expected, result, 0.001);
@@ -233,7 +233,7 @@ public final strictfp class NTv2Test extends DatumShiftTestCase {
          * when `SpecializableTransform` has not been able to find the most appropriate grid.
          */
         grid.interpolateInCell(indices[0], indices[1], result);
-        result[0] = position[0] - result[0] * cellSize;
+        result[0] = position[0] - result[0] * cellSize;                     // Positive translation
is toward west.
         result[1] = position[1] + result[1] * cellSize;
         assertArrayEquals("interpolateInCell", expected, result, Formulas.ANGULAR_TOLERANCE
* DEGREES_TO_SECONDS);
     }


Mime
View raw message