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: Make wraparound checks more robust by comparing "transform without wraparound" with "transform with wraparound".
Date Wed, 09 Sep 2020 10:37:50 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 625263a  Make wraparound checks more robust by comparing "transform without wraparound"
with "transform with wraparound".
625263a is described below

commit 625263abdc5f84d13b9a680013f626b47dea498b
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Wed Sep 9 12:32:16 2020 +0200

    Make wraparound checks more robust by comparing "transform without wraparound" with "transform
with wraparound".
---
 .../coverage/grid/CoordinateOperationFinder.java   | 258 +++++++++++++++------
 .../internal/referencing/WraparoundTransform.java  | 115 ++++++++-
 .../referencing/WraparoundTransformTest.java       |  20 ++
 .../sis/storage/geotiff/GridGeometryBuilder.java   |   2 +-
 4 files changed, 319 insertions(+), 76 deletions(-)

diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
index 35e111f..ea4f59c 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/CoordinateOperationFinder.java
@@ -32,7 +32,6 @@ import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.util.collection.BackingStoreException;
-import org.apache.sis.internal.util.Numerics;
 import org.apache.sis.internal.system.DefaultFactories;
 import org.apache.sis.internal.referencing.CoordinateOperations;
 import org.apache.sis.internal.referencing.WraparoundTransform;
@@ -44,19 +43,34 @@ import org.apache.sis.measure.Quantities;
 import org.apache.sis.measure.Units;
 import org.apache.sis.util.ArraysExt;
 import org.apache.sis.referencing.CRS;
+import org.apache.sis.util.logging.Logging;
+import org.apache.sis.internal.system.Modules;
 
 
 /**
- * Finds a transform from points expressed in the CRS of a source coverage to points in the
CRS of a target coverage.
+ * Finds a transform from grid cells in a source coverage to geospatial positions in the
CRS of a target coverage.
  * This class differs from {@link CRS#findOperation CRS#findOperation(…)} because of the
gridded aspect of inputs.
  * {@linkplain GridGeometry grid geometries} give more information about how referencing
is applied on datasets.
- * With them, we can detect dimensions where target coordinates are constrained to constant
values
- * because the {@linkplain GridExtent#getSize(int) grid size} is only one cell in those dimensions.
- * This is an important difference because it allows us to find operations normally impossible,
- * where we still can produce an operation to a target CRS even if some dimensions have no
corresponding source CRS.
+ * With them, this class provides three additional benefits:
  *
- * <p><b>Note:</b> this class does not provide complete chain of transformation
from grid to grid.
- * It provides only the operation from the CRS of source to the CRS of destination.</p>
+ * <ul class="verbose">
+ *   <li>Detect dimensions where target coordinates are constrained to constant values
because
+ *     the {@linkplain GridExtent#getSize(int) grid size} is only one cell in those dimensions.
+ *     This is an important because it makes possible to find operations normally impossible:
+ *     we can still produce an operation to a target CRS even if some dimensions have no
corresponding
+ *     source CRS.</li>
+ *
+ *   <li>Detect how to handle wraparound axes. For example if a raster spans 150°
to 200° of longitude,
+ *     this class understands that -170° of longitude should be translated to 190° for
thar particular
+ *     raster. This will work even if the minimum and maximum values declared in the longitude
axis do
+ *     not match that range.</li>
+ *
+ *   <li>Use the area of interest and grid resolution for refining the coordinate operation
between two CRS.</li>
+ * </ul>
+ *
+ * <b>Note:</b> this class does not provide the complete chain of operations
from grid to grid.
+ * It provides only the operation from <em>cell indices</em> in source grid to
<em>coordinates in the CRS</em>
+ * of destination grid. Callers must add the last step (conversion from target CRS to cell
indices) themselves.
  *
  * @author  Alexis Manin (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
@@ -78,28 +92,53 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     /**
      * The target coordinate values, computed only if needed. This is computed by {@link
#get()}, which is
      * itself invoked indirectly by {@link org.apache.sis.referencing.operation.CoordinateOperationFinder}.
-     * The value is cached in case {@code get()} is invoked multiple times during the same
finder execution.
+     * The value is cached in case {@code get()} is invoked many times during the same finder
execution.
+     *
+     * @see #get()
      */
     private double[] coordinates;
 
     /**
      * The coordinate operation from source to target CRS, computed when first needed.
+     *
+     * @see #gridToCRS(PixelInCell)
      */
     private CoordinateOperation operation;
 
     /**
      * The {@link #operation} transform together with {@link WraparoundTransform} if needed.
      * The wraparound is used for handling images crossing the anti-meridian.
+     *
+     * @see #gridToCRS(PixelInCell)
      */
     private MathTransform forwardOp;
 
     /**
      * Inverse of {@link #operation} transform together with {@link WraparoundTransform}
if needed.
      * The wraparound is used for handling images crossing the anti-meridian.
+     *
+     * <p>Contrarily to other {@link MathTransform} fields, determination of this {@code
inverseOp}
+     * transform should make an effort for checking if wraparound is really needed. This
is because
+     * {@code inverseOp} will be used much more extensively (for every pixels) than other
transforms.</p>
+     *
+     * @see #isWraparoundNeeded
+     * @see #isWraparoundApplied
+     * @see #inverse(MathTransform)
      */
     private MathTransform inverseOp;
 
     /**
+     * Transform from the target CRS to the source grid, with {@link WraparoundTransform}
applied if needed.
+     * This is the concatenation of {@code inverseOp} with inverse of {@link #source} "grid
to CRS", except
+     * for more conservative application of wraparound (i.e. determination of this transform
should not do
+     * the simplification effort mentioned in {@link #inverseOp}).
+     *
+     * @see #inverse(MathTransform)
+     * @see #applyWraparound(MathTransform)
+     */
+    private MathTransform crsToGrid;
+
+    /**
      * Whether {@link #inverseOp} needs to include a {@link WraparoundTransform} step. We
do this check
      * only for {@link #inverseOp} because it is the transform which will be executed for
every pixels.
      * By contrast, {@link #forwardOp} will systematically contains a {@link WraparoundTransform}
step
@@ -109,6 +148,20 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     private boolean isWraparoundNeeded;
 
     /**
+     * Whether {@link WraparoundTransform} has been applied on {@link #inverseOp}. This field
complements
+     * {@link #isWraparoundNeeded} because a delay may exist between the time we detected
that wraparound
+     * is needed and the time we added the necessary operation steps.
+     *
+     * <p>Note that despite this field name, a {@code true} value does not mean that
{@link #inverseOp}
+     * and {@link #crsToGrid} really contains some {@link WraparoundTransform} steps. It
only means that
+     * {@link WraparoundTransform#forDomainOfUse WraparoundTransform.forDomainOfUse(…)}
has been invoked.
+     * That method may have decided to not insert any wraparound steps.</p>
+     *
+     * @see #applyWraparound(MathTransform)
+     */
+    private boolean isWraparoundApplied;
+
+    /**
      * The factory to use for {@link MathTransform} creations. For now this is fixed to SIS
implementation.
      * But it may become a configurable reference in a future version if useful.
      *
@@ -155,26 +208,33 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     }
 
     /**
-     * Computes the transform from the grid coordinates of the source to geospatial coordinates
of the target.
+     * Computes the transform from “grid coordinates of the source” to “geospatial
coordinates of the target”.
      * It may be the identity operation. We try to take envelopes in account because the
operation choice may
      * depend on the geographic area.
      *
-     * <p>The transform returned by this method always apply wraparounds on every axes
having wraparound range.
-     * This method does not check if wraparounds actually happen. This is okay because this
transform is used
-     * only for transforming envelopes; it is not used for transforming pixel coordinates.
By contrast pixel
-     * transform (computed by {@link #inverse(MathTransform)}) will need to perform more
extensive check.</p>
+     * <p>The transform returned by this method applies wraparound checks systematically
on every axes having
+     * wraparound range. This method does not verify whether those checks are needed (i.e.
whether wraparound
+     * can possibly happen). This is okay because this transform is used only for transforming
envelopes;
+     * it is not used for transforming pixel coordinates.<p>
+     *
+     * <h4>Implementation note</h4>
+     * After invocation of this method, the following fields are valid:
+     * <ul>
+     *   <li>{@link #operation} — cached for {@link #inverse(MathTransform)} usage.</li>
+     *   <li>{@link #forwardOp} — cached for next invocation of this {@code gridToCRS(…)}
method.</li>
+     * </ul>
      *
      * @param  anchor  whether the operation is between cell centers or cell corners.
-     * @return operation from source CRS to target CRS, or {@code null} if CRS are unspecified
or equivalent.
+     * @return operation from source grid indices to target geospatial coordinates.
      * @throws FactoryException if no operation can be found between the source and target
CRS.
      * @throws TransformException if some coordinates can not be transformed to the specified
target.
      * @throws IncompleteGridGeometryException if required CRS or a "grid to CRS" information
is missing.
      */
     final MathTransform gridToCRS(final PixelInCell anchor) throws FactoryException, TransformException
{
         /*
-         * If `coordinates` is non-null, it means that the `get()` method has been invoked
during previous
-         * call to this `gridToCRS(…)` method, which implies that `operation` depends on
the `anchor` value.
-         * In such case we need to discard the previous `operation` value and recompute it.
+         * If `coordinates` is non-null, it means that `CoordinateOperationContext.getConstantCoordinates()`
+         * has been invoked during previous `gridToCRS(…)` execution, which implies that
`operation` depends
+         * on the`anchor` value. In such case we will need to recompute all fields that depend
on `operation`.
          */
         this.anchor = anchor;
         if (coordinates != null) {
@@ -182,6 +242,8 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
             operation   = null;
             forwardOp   = null;
             inverseOp   = null;
+            crsToGrid   = null;
+            // Do not clear `isWraparoundNeeded`; its value is still valid.
         }
         /*
          * Get the operation performing the change of CRS. If more than one operation is
defined for
@@ -207,13 +269,18 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
                          * the check for wraparound axes, which is necessary even if the
transform is identity.
                          */
                         DefaultGeographicBoundingBox areaOfInterest = null;
-                        if (sourceEnvelope != null || targetEnvelope != null) {
+                        if (sourceEnvelope != null || targetEnvelope != null) try {
                             areaOfInterest = new DefaultGeographicBoundingBox();
                             areaOfInterest.setBounds(targetEnvelope != null ? targetEnvelope
: sourceEnvelope);
+                        } catch (TransformException e) {
+                            areaOfInterest = null;
+                            recoverableException("gridToCRS", e);
                         }
                         operation = CRS.findOperation(sourceCRS, target.getCoordinateReferenceSystem(),
areaOfInterest);
                     }
                 }
+            } catch (BackingStoreException e) {                         // May be thrown
by getConstantCoordinates().
+                throw e.unwrapOrRethrow(TransformException.class);
             } finally {
                 CoordinateOperations.CONSTANT_COORDINATES.remove();
             }
@@ -252,71 +319,106 @@ wraparound: if (mayRequireWraparound(cs)) {
      * This is equivalent to invoking {@link MathTransform#inverse()} on the transform returned
by
      * {@link #gridToCRS(PixelInCell)}, except in the way wraparounds are handled.
      *
-     * <p>The {@code gridToCRS} argument is the value returned by last call to {@link
#gridToCRS(PixelInCell)}.
-     * That transform is used for testing whether wraparound is needed for calculation of
source coordinates.
-     * This method performs a more extensive check than {@code gridToCRS(…)} because
-     * the transform returned by {@code inverse(…)} will be applied on every pixels of
destination image.
-     * The argument should be non-null only the first time that this {@code inverse(…)}
method is invoked,
-     * preferably with a transform mapping {@link PixelInCell#CELL_CORNER}.
-     * Subsequent invocations will use the {@link #isWraparoundNeeded} status determined
by first invocation.</p>
+     * <p>The {@code gridToCRS} argument controls whether to perform a more extensive
check of wraparound occurrence.
+     * If {@code null}, the result of last check is reused. If non-null, then it shall be
the value returned by last
+     * call to {@link #gridToCRS(PixelInCell)}, preferably with a transform mapping {@link
PixelInCell#CELL_CORNER}.
+     * The argument should be non-null the first time that {@code inverse(…)} is invoked
and {@code null} next time.</p>
+     *
+     * @param  gridToCRS  result of previous call to {@link #gridToCRS(PixelInCell)}, or
{@code null}
+     *                    for reusing the {@link #isWraparoundNeeded} result of previous
invocation.
+     * @return operation from target geospatial coordinates to source grid indices.
+     * @throws FactoryException if some transform steps can not be concatenated.
+     * @throws TransformException if some coordinates can not be transformed.
      */
     final MathTransform inverse(final MathTransform gridToCRS) throws FactoryException, TransformException
{
         final MathTransform tr = source.getGridToCRS(anchor).inverse();
         if (operation == null) {
             return tr;
         }
-        if (inverseOp == null) {
-            inverseOp = operation.getMathTransform().inverse();
-check:      if (gridToCRS != null) {
-                /*
-                 * Transform all corners of source extent to the destination CRS, then back
to source grid coordinates.
-                 * We do not concatenate the forward and inverse transforms because we do
not want MathTransformFactory
-                 * to simplify the transformation chain (e.g. replacing "Mercator → Inverse
of Mercator" by an identity
-                 * transform), because such simplification would erase wraparound effects.
-                 */
-                final MathTransform crsToGrid = mtFactory.createConcatenatedTransform(inverseOp,
tr);
-                final GridExtent extent = source.getExtent();
-                final int dimension = extent.getDimension();
-                final double[] src = new double[dimension];
-                final double[] tgt = new double[Math.max(dimension, gridToCRS.getTargetDimensions())];
-                for (long maskOfUppers = Numerics.bitmask(dimension); --maskOfUppers != 0;)
{
-                    long bit = 1;
-                    for (int i=0; i<dimension; i++) {
-                        src[i] = (maskOfUppers & bit) != 0 ? extent.getHigh(i) + 1d :
extent.getLow(i);
-                        bit <<= 1;
-                    }
-                    gridToCRS.transform(src, 0, tgt, 0, 1);
-                    crsToGrid.transform(tgt, 0, tgt, 0, 1);
-                    for (int i=0; i<dimension; i++) {
-                        /*
-                         * Do not consider NaN as a need for wraparound because NaN values
occur when
-                         * the operation between the two CRS reduce the number of dimensions.
-                         *
-                         * TODO: we may require a more robust check for NaN values.
-                         */
-                        if (Math.abs(tgt[i] - src[i]) > 1) {
-                            isWraparoundNeeded = true;
-                            break check;
-                        }
-                    }
-                }
-                return crsToGrid;
-            }
+        if (inverseOp != null) {
+            // Here, `inverseOp` contains the wraparound step if needed.
+            return mtFactory.createConcatenatedTransform(inverseOp, tr);
+        }
+        /*
+         * Need to compute transform with wraparound checks, but contrarily to `gridToCRS(…)`
we do not want
+         * `WraparoundTransform` to be systematically inserted. This is for performance reasons,
because the
+         * transform returned by this method will be applied on every pixels of destination
image.  We start
+         * with transforms without wraparound, and modify those fields later depending on
whether wraparound
+         * appears to be necessary or not.
+         */
+        final MathTransform inverseNoWrap   = operation.getMathTransform().inverse();
+        final MathTransform crsToGridNoWrap = mtFactory.createConcatenatedTransform(inverseNoWrap,
tr);
+        inverseOp = inverseNoWrap;
+        crsToGrid = crsToGridNoWrap;
+        isWraparoundApplied = false;
+check:  if (gridToCRS != null) try {
             /*
-             * Potentially append a wraparound if we determined (either in this method invocation
-             * or in a previous invocation of this method) that it is necessary.
+             * We will do a more extensive check by converting all corners of source grid
to the target CRS,
+             * then convert back to the source grid and see if coordinates match. Only if
coordinates do not
+             * match, `WraparoundTransform.isNeeded(…)` will request a `crsToGrid` transform
which includes
+             * wraparound steps in order to check if it improves the results. By using a
`Supplier`,
+             * we avoid creating `WraparoundTransform` in the common case where it is not
needed.
              */
-            if (isWraparoundNeeded) {
-                final CoordinateSystem cs = operation.getSourceCRS().getCoordinateSystem();
-                if (mayRequireWraparound(cs)) {
-                    final DirectPosition median = median(source);
-                    if (median != null) {
-                        inverseOp = WraparoundTransform.forDomainOfUse(mtFactory, inverseOp,
cs, median);
-                    }
+            final Supplier<MathTransform> withWraparound = () -> {
+                try {
+                    return applyWraparound(tr);
+                } catch (FactoryException | TransformException e) {
+                    throw new BackingStoreException(e);
+                }
+            };
+            final GridExtent extent = source.getExtent();
+            final boolean isCorner = (anchor == PixelInCell.CELL_CORNER);
+            isWraparoundNeeded = WraparoundTransform.isNeeded(gridToCRS, crsToGridNoWrap,
withWraparound, (dim) -> {
+                long cc;                                          // Coordinate of a corner.
+                if (dim < 0) {
+                    cc = extent.getLow(~dim);
+                } else {
+                    cc = extent.getHigh(dim);                     // Inclusive.
+                    if (isCorner && cc != Long.MAX_VALUE) cc++;   // Make exclusive.
+                }
+                return cc;
+            });
+        } catch (BackingStoreException e) {
+            throw e.unwrapOrRethrow(FactoryException.class);
+            // TransformException is handled in `WraparoundTransform.isNeeded(…)` instead
of here.
+        }
+        /*
+         * At this point we determined whether wraparound is needed. The `inverseOp` and
`crsToGrid` fields
+         * may have been updated as part of this determination process, but not necessarily.
If wraparounds
+         * are needed, reuse available calculation results (completing them if needed). 
Otherwise rollback
+         * any change that `applyWraparound(…)` may have done to `inverseOp` and `crsToGrid`.
+         */
+        if (isWraparoundNeeded) {
+            return applyWraparound(tr);
+        } else {
+            inverseOp = inverseNoWrap;
+            crsToGrid = crsToGridNoWrap;
+            return crsToGridNoWrap;
+        }
+    }
+
+    /**
+     * If not already done, inserts {@link WraparoundTransform} steps into {@link #inverseOp}
and {@link #crsToGrid}
+     * transforms. The transform from geospatial target coordinates to source grid indices
is returned for convenience.
+     *
+     * @param  tr  value of {@code source.getGridToCRS(anchor).inverse()}.
+     * @return transform from geospatial target coordinates to source grid indices.
+     * @throws FactoryException if some transform steps can not be concatenated.
+     * @throws TransformException if some coordinates can not be transformed.
+     */
+    private MathTransform applyWraparound(final MathTransform tr) throws FactoryException,
TransformException {
+        if (!isWraparoundApplied) {
+            isWraparoundApplied = true;
+            final CoordinateSystem cs = operation.getSourceCRS().getCoordinateSystem();
+            if (mayRequireWraparound(cs)) {
+                final DirectPosition median = median(source);
+                if (median != null) {
+                    inverseOp = WraparoundTransform.forDomainOfUse(mtFactory, inverseOp,
cs, median);
+                    crsToGrid = mtFactory.createConcatenatedTransform(inverseOp, tr);
                 }
             }
         }
-        return mtFactory.createConcatenatedTransform(inverseOp, tr);
+        return crsToGrid;
     }
 
     /**
@@ -416,4 +518,14 @@ check:      if (gridToCRS != null) {
             processor.setPositionalAccuracyHints(hints);                        // Null elements
will be ignored.
         }
     }
+
+    /**
+     * Invoked when an ignorable exception occurred.
+     *
+     * @param  caller  the method where the exception occurred.
+     * @param  e       the ignorable exception.
+     */
+    private static void recoverableException(final String caller, final Exception e) {
+        Logging.recoverableException(Logging.getLogger(Modules.RASTER), CoordinateOperationFinder.class,
caller, e);
+    }
 }
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
index bb323d3..9917d92 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/internal/referencing/WraparoundTransform.java
@@ -17,6 +17,8 @@
 package org.apache.sis.internal.referencing;
 
 import java.util.List;
+import java.util.function.Supplier;
+import java.util.function.IntToDoubleFunction;
 import org.opengis.util.FactoryException;
 import org.opengis.geometry.DirectPosition;
 import org.opengis.referencing.cs.RangeMeaning;
@@ -27,6 +29,7 @@ import org.opengis.referencing.operation.MathTransform;
 import org.opengis.referencing.operation.MathTransformFactory;
 import org.opengis.referencing.operation.Matrix;
 import org.opengis.referencing.operation.NoninvertibleTransformException;
+import org.opengis.referencing.operation.TransformException;
 import org.apache.sis.referencing.factory.InvalidGeodeticParameterException;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
@@ -39,6 +42,7 @@ import org.apache.sis.util.ArgumentChecks;
 import org.apache.sis.util.ComparisonMode;
 import org.apache.sis.util.resources.Errors;
 import org.apache.sis.util.logging.Logging;
+import org.apache.sis.util.collection.BackingStoreException;
 
 
 /**
@@ -69,6 +73,38 @@ import org.apache.sis.util.logging.Logging;
  */
 public final class WraparoundTransform extends AbstractMathTransform {
     /**
+     * Number of dimensions of the first cached {@link WraparoundTransform} instance.
+     */
+    private static final int FIRST_CACHED_DIMENSION = 2;
+
+    /**
+     * Frequently used {@link WraparoundTransform} instances.
+     * The {@link #dimension} value ranges from {@value #FIRST_CACHED_DIMENSION} to 4 inclusive,
+     * and the {@link #wraparoundDimension} value alternates between 0 and 1.
+     */
+    private static final WraparoundTransform[] CACHE = new WraparoundTransform[6];
+    static {
+        for (int i=0; i<CACHE.length; i++) {
+            final int dimension = (i >>> 1) + FIRST_CACHED_DIMENSION;
+            CACHE[i] = new WraparoundTransform(dimension, i & 1);
+        }
+    }
+
+    /**
+     * Returns a transform with a wraparound behavior in the given dimension.
+     * Input and output values in the wraparound dimension shall be normalized in the [0
… 1] range.
+     */
+    static WraparoundTransform create(final int dimension, final int wraparoundDimension)
{
+        if ((wraparoundDimension & ~1) == 0) {
+            final int i = ((dimension - FIRST_CACHED_DIMENSION) << 1) | wraparoundDimension;
+            if (i >= 0 && i < CACHE.length) {
+                return CACHE[i];
+            }
+        }
+        return new WraparoundTransform(dimension, wraparoundDimension);
+    }
+
+    /**
      * The dimension of source and target coordinates.
      */
     private final int dimension;
@@ -95,7 +131,7 @@ public final class WraparoundTransform extends AbstractMathTransform {
     private WraparoundTransform redim(final int n) {
         if (n == dimension) return this;
         if (n >= wraparoundDimension) return null;
-        return new WraparoundTransform(n, wraparoundDimension);
+        return create(n, wraparoundDimension);
     }
 
     /**
@@ -158,6 +194,8 @@ public final class WraparoundTransform extends AbstractMathTransform {
 
     /**
      * Creates a transform with a "wrap around" behavior in the given dimension.
+     * The wraparound is implemented by a concatenation of affine transform before
+     * and after the {@link WraparoundTransform} instance.
      *
      * @param  factory              the factory to use for creating math transforms.
      * @param  dimension            the number of source and target dimensions.
@@ -183,7 +221,7 @@ public final class WraparoundTransform extends AbstractMathTransform {
             m.setElement(wraparoundDimension, wraparoundDimension, span);
             m.setElement(wraparoundDimension, dimension, Double.isNaN(median) ? minimum :
median - span/2);
             final MathTransform denormalize = factory.createAffineTransform(m);
-            final WraparoundTransform wraparound = new WraparoundTransform(dimension, wraparoundDimension);
+            final WraparoundTransform wraparound = create(dimension, wraparoundDimension);
             try {
                 return factory.createConcatenatedTransform(denormalize.inverse(),
                        factory.createConcatenatedTransform(wraparound, denormalize));
@@ -196,6 +234,79 @@ public final class WraparoundTransform extends AbstractMathTransform
{
     }
 
     /**
+     * Verifies whether wraparound is needed for a "CRS to grid" transform.
+     * This method converts coordinates of all corners of a grid to an arbitrary CRS, then
back to the grid.
+     * The forward and backward transforms are not exactly the inverse of each other: the
forward transform
+     * applies wraparounds while the backward transform does not. By checking whether or
not the roundtrip
+     * result is equal to the original coordinates, we determine if wraparound is necessary
or not.
+     *
+     * <p>The coordinates to test are provided by {@code gridExtent} using the following
convention.
+     * For function argument {@code int} <var>i</var>:</p>
+     * <ul>
+     *   <li>If <var>i</var> ≥ 0, return the upper value at dimension
{@code i}.</li>
+     *   <li>If <var>i</var> &lt; 0, return the lower value at dimension
{@code ~i} (tild operator, not minus).</li>
+     * </ul>
+     *
+     * @param  gridToCRS       a transform to an arbitrary CRS with handling of wraparound
axes.
+     * @param  crsToGrid       inverse of {@code gridToCRS} but without handling of wraparound
axes.
+     * @param  withWraparound  provider of a transform equivalent to {@code crsToGrid} but
with wraparound applied.
+     * @param  gridExtent      provider of grid corner coordinates (see above javadoc).
+     * @return whether wraparound transform seems needed.
+     * @throws TransformException if an error occurred while transforming coordinates.
+     */
+    public static boolean isNeeded(final MathTransform gridToCRS, final MathTransform crsToGrid,
+            final Supplier<MathTransform> withWraparound, final IntToDoubleFunction
gridExtent)
+            throws TransformException
+    {
+        MathTransform  watr = null;                     // Transform with wraparound (fetched
only if needed).
+        final int dimension = gridToCRS.getSourceDimensions();
+        final double[]  src = new double[dimension];    // Coordinates of a corner.
+        final double[]  nwp = new double[dimension];    // Roundtrip result without wraparound.
+        final double[]  tgt = new double[Math.max(gridToCRS.getTargetDimensions(), dimension)];
+        long  maskOfUppers  = Numerics.bitmask(dimension);
+        long  maskOfChanges = 0;
+        do {
+            maskOfChanges ^= --maskOfUppers;            // Source coordinates that changed
since last iteration.
+            do {
+                final int i = Long.numberOfTrailingZeros(maskOfChanges);
+                final long bit = 1L << i;
+                src[i] = gridExtent.applyAsDouble((maskOfUppers & bit) == 0 ? ~i : i);
+                maskOfChanges &= ~bit;
+            } while (maskOfChanges != 0);               // Iterate only over the coordinates
that changed.
+            maskOfChanges = maskOfUppers;
+            /*
+             * Transform corner from the source extent to the destination CRS, then back
to source grid coordinates.
+             * We do not concatenate the forward and inverse transforms because we do not
want MathTransformFactory
+             * to simplify the transformation chain (e.g. replacing "Mercator → Inverse
of Mercator" by an identity
+             * transform), because such simplification would erase wraparound effects.
+             */
+            boolean watrApplied = false;
+            gridToCRS.transform(src, 0, tgt, 0, 1);
+            crsToGrid.transform(tgt, 0, nwp, 0, 1);
+            for (int i=0; i<dimension; i++) {
+                final double error = Math.abs(nwp[i] - src[i]);
+                if (!(error <= 1)) {                                // Use `!` for catching
NaN.
+                    if (!watrApplied) try {
+                        watrApplied = true;
+                        if (watr == null) watr = withWraparound.get();
+                        watr.transform(tgt, 0, tgt, 0, 1);
+                    } catch (BackingStoreException e) {
+                        throw e.unwrapOrRethrow(TransformException.class);
+                    }
+                    /*
+                     * Do not consider NaN in `tgt` as a need for wraparound because NaN
values
+                     * occur when an operation between two CRS reduces the number of dimensions.
+                     */
+                    if (Math.abs(tgt[i] - src[i]) < (error <= Double.MAX_VALUE ? error
: 1)) {
+                        return true;
+                    }
+                }
+            }
+        } while (maskOfUppers != 0);
+        return false;
+    }
+
+    /**
      * Gets the dimension of input points.
      *
      * @return the dimension of input points.
diff --git a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java
b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java
index f66d662..ec8b366 100644
--- a/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java
+++ b/core/sis-referencing/src/test/java/org/apache/sis/internal/referencing/WraparoundTransformTest.java
@@ -45,6 +45,26 @@ import static org.opengis.test.Assert.*;
  */
 public final strictfp class WraparoundTransformTest extends TestCase {
     /**
+     * Tests {@link WraparoundTransform} cache.
+     */
+    @Test
+    public void testCache() {
+        final WraparoundTransform t1, t2, t3, t4;
+        assertSame   (WraparoundTransform.create(3, 0), t1 = WraparoundTransform.create(3,
0));
+        assertNotSame(WraparoundTransform.create(3, 0), t2 = WraparoundTransform.create(3,
1));
+        assertNotSame(WraparoundTransform.create(3, 0), t3 = WraparoundTransform.create(2,
0));
+        assertNotSame(WraparoundTransform.create(3, 0), t4 = WraparoundTransform.create(3,
2));
+        assertEquals(3, t1.getSourceDimensions());
+        assertEquals(3, t2.getSourceDimensions());
+        assertEquals(2, t3.getSourceDimensions());
+        assertEquals(3, t4.getSourceDimensions());
+        assertEquals(0, t1.wraparoundDimension);
+        assertEquals(1, t2.wraparoundDimension);
+        assertEquals(0, t3.wraparoundDimension);
+        assertEquals(2, t4.wraparoundDimension);
+    }
+
+    /**
      * Tests wraparound on one axis.
      *
      * @throws FactoryException if the transform can not be created.
diff --git a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
index 6c92d9e..b780694 100644
--- a/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
+++ b/storage/sis-geotiff/src/main/java/org/apache/sis/storage/geotiff/GridGeometryBuilder.java
@@ -46,7 +46,7 @@ import org.apache.sis.math.Vector;
  * The coordinate reference system part is built by {@link CRSBuilder}.
  *
  * <h2>Pixel center versus pixel corner</h2>
- * The policy about whether the conversion map pixel corner or pixel center if GeoTIFF files
does not seem
+ * The policy about whether the conversion maps pixel corner or pixel center in GeoTIFF files
does not seem
  * totally clear. But the practice at least with GDAL seems to consider the following as
equivalent:
  *
  * {@preformat text


Mime
View raw message