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: Initial support of wraparound axes during raster resampling.
Date Mon, 07 Sep 2020 12:12:40 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 cc12164  Initial support of wraparound axes during raster resampling.
cc12164 is described below

commit cc12164a1be36fdfa7d1ba736460b22ca65c6445
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Mon Sep 7 14:12:13 2020 +0200

    Initial support of wraparound axes during raster resampling.
---
 .../coverage/grid/CoordinateOperationFinder.java   | 232 +++++++++++++++++++--
 .../sis/coverage/grid/ResampledGridCoverage.java   |   8 +-
 .../java/org/apache/sis/image/PlanarImage.java     |   2 +-
 .../coverage/grid/ResampledGridCoverageTest.java   |  39 +++-
 4 files changed, 255 insertions(+), 26 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 b5fdaee..35e111f 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
@@ -22,21 +22,28 @@ import javax.measure.Quantity;
 import javax.measure.quantity.Length;
 import org.opengis.geometry.Envelope;
 import org.opengis.util.FactoryException;
+import org.opengis.geometry.DirectPosition;
 import org.opengis.referencing.datum.PixelInCell;
+import org.opengis.referencing.cs.RangeMeaning;
+import org.opengis.referencing.cs.CoordinateSystem;
 import org.opengis.referencing.crs.CoordinateReferenceSystem;
 import org.opengis.referencing.operation.CoordinateOperation;
 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;
 import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox;
 import org.apache.sis.geometry.Envelopes;
+import org.apache.sis.geometry.GeneralDirectPosition;
 import org.apache.sis.image.ImageProcessor;
 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.referencing.operation.transform.MathTransforms;
 
 
 /**
@@ -69,7 +76,9 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     private final GridGeometry source, target;
 
     /**
-     * The target coordinate values, computed only if needed.
+     * 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.
      */
     private double[] coordinates;
 
@@ -79,6 +88,35 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     private CoordinateOperation operation;
 
     /**
+     * The {@link #operation} transform together with {@link WraparoundTransform} if needed.
+     * The wraparound is used for handling images crossing the anti-meridian.
+     */
+    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.
+     */
+    private MathTransform inverseOp;
+
+    /**
+     * 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
+     * because we use it only for transforming envelopes and for the test that determines
the value of
+     * this {@code isWraparoundNeeded} flag.
+     */
+    private boolean isWraparoundNeeded;
+
+    /**
+     * 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.
+     *
+     * @see org.apache.sis.coverage.grid.ImageRenderer#mtFactory
+     */
+    private final MathTransformFactory mtFactory;
+
+    /**
      * Creates a new finder.
      *
      * @param  source  the grid geometry which is the source of the coordinate operation
to find.
@@ -87,6 +125,8 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     CoordinateOperationFinder(final GridGeometry source, final GridGeometry target) {
         this.source = source;
         this.target = target;
+        mtFactory = DefaultFactories.forBuildin(MathTransformFactory.class,
+                org.apache.sis.referencing.operation.transform.DefaultMathTransformFactory.class).caching(false);
     }
 
     /**
@@ -119,6 +159,11 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
      * 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>
+     *
      * @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.
      * @throws FactoryException if no operation can be found between the source and target
CRS.
@@ -127,17 +172,26 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
      */
     final MathTransform gridToCRS(final PixelInCell anchor) throws FactoryException, TransformException
{
         /*
-         * If `coordinates` is non-null, it means that the `get()` method has been invoked
in previous call
-         * to this `cornerToCRS(…)` method, which implies that `operation` depends on the
`anchor` value.
+         * 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.
          */
         this.anchor = anchor;
-        if (operation == null || coordinates != null) {
-            operation   = null;
+        if (coordinates != null) {
             coordinates = null;
-            final CoordinateReferenceSystem sourceCRS = getSourceCRS();
-            final CoordinateReferenceSystem targetCRS = target.isDefined(GridGeometry.CRS)
?
-                                                        target.getCoordinateReferenceSystem()
: sourceCRS;
+            operation   = null;
+            forwardOp   = null;
+            inverseOp   = null;
+        }
+        /*
+         * Get the operation performing the change of CRS. If more than one operation is
defined for
+         * the given pair of CRS, select the most appropriate operation for the area of interest.
+         *
+         * TODO: specify also the desired resolution, computed from target grid geometry.
+         *       It will require more direct use of `CoordinateOperationContext`.
+         *       As a side effect, we could remove `CONSTANT_COORDINATES` hack.
+         */
+        if (operation == null) {
             final Envelope sourceEnvelope = source.envelope;
             final Envelope targetEnvelope = target.envelope;
             try {
@@ -145,13 +199,20 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
                 if (sourceEnvelope != null && targetEnvelope != null) {
                     operation = Envelopes.findOperation(sourceEnvelope, targetEnvelope);
                 }
-                if (operation == null && sourceCRS != null) {
-                    DefaultGeographicBoundingBox areaOfInterest = null;
-                    if (sourceEnvelope != null || targetEnvelope != null) {
-                        areaOfInterest = new DefaultGeographicBoundingBox();
-                        areaOfInterest.setBounds(targetEnvelope != null ? targetEnvelope
: sourceEnvelope);
+                if (operation == null && target.isDefined(GridGeometry.CRS)) {
+                    final CoordinateReferenceSystem sourceCRS = getSourceCRS();
+                    if (sourceCRS != null) {
+                        /*
+                         * Unconditionally create operation even if CRS are the same. A non-null
operation trig
+                         * the check for wraparound axes, which is necessary even if the
transform is identity.
+                         */
+                        DefaultGeographicBoundingBox areaOfInterest = null;
+                        if (sourceEnvelope != null || targetEnvelope != null) {
+                            areaOfInterest = new DefaultGeographicBoundingBox();
+                            areaOfInterest.setBounds(targetEnvelope != null ? targetEnvelope
: sourceEnvelope);
+                        }
+                        operation = CRS.findOperation(sourceCRS, target.getCoordinateReferenceSystem(),
areaOfInterest);
                     }
-                    operation = CRS.findOperation(sourceCRS, targetCRS, areaOfInterest);
                 }
             } finally {
                 CoordinateOperations.CONSTANT_COORDINATES.remove();
@@ -161,19 +222,148 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
          * The following line may throw IncompleteGridGeometryException, which is desired
because
          * if that transform is missing, we can not continue (we have no way to guess it).
          */
-        MathTransform tr = source.getGridToCRS(anchor);
-        if (operation != null) {
-            tr = MathTransforms.concatenate(tr, operation.getMathTransform());
+        final MathTransform tr = source.getGridToCRS(anchor);
+        if (operation == null) {
+            return tr;
+        }
+        /*
+         * At this point we have a "grid → source CRS" transform. Append a "source CRS
→ target CRS" transform,
+         * which may be identity. A wraparound may be applied for keeping target coordinates
inside the expected
+         * target domain.
+         */
+        if (forwardOp == null) {
+            forwardOp = operation.getMathTransform();
+            final CoordinateSystem cs = operation.getTargetCRS().getCoordinateSystem();
+wraparound: if (mayRequireWraparound(cs)) {
+                DirectPosition median = median(target);
+                if (median == null) {
+                    median = median(source);
+                    if (median == null) break wraparound;
+                    median = forwardOp.transform(median, null);
+                }
+                forwardOp = WraparoundTransform.forDomainOfUse(mtFactory, forwardOp, cs,
median);
+            }
+        }
+        return mtFactory.createConcatenatedTransform(tr, forwardOp);
+    }
+
+    /**
+     * Returns the inverse of the transform returned by last call to {@link #gridToCRS(PixelInCell)}.
+     * 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>
+     */
+    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;
+            }
+            /*
+             * Potentially append a wraparound if we determined (either in this method invocation
+             * or in a previous invocation of this method) that it is necessary.
+             */
+            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);
+                    }
+                }
+            }
         }
-        return tr;
+        return mtFactory.createConcatenatedTransform(inverseOp, tr);
     }
 
     /**
-     * Invoked when the target CRS have some dimensions that the source CRS does not have.
+     * Returns {@code true} if a transform to the specified target coordinate system may
require handling
+     * of wraparound axes. A {@code true} return value does not mean that wraparounds actually
happen;
+     * it only means that more expensive checks will be required.
+     */
+    private static boolean mayRequireWraparound(final CoordinateSystem cs) {
+        final int dimension = cs.getDimension();
+        for (int i=0; i<dimension; i++) {
+            if (RangeMeaning.WRAPAROUND.equals(cs.getAxis(i).getRangeMeaning())) {
+                return true;
+            }
+        }
+        return false;
+    }
+
+    /**
+     * Returns the point of interest converted to the Coordinate Reference System.
+     * If the grid does not define a point of interest or does not define a CRS,
+     * then this method returns {@code null}.
+     */
+    private static DirectPosition median(final GridGeometry grid) throws TransformException
{
+        if (grid.isDefined(GridGeometry.EXTENT | GridGeometry.GRID_TO_CRS)) {
+            final double[] poi = grid.getExtent().getPointOfInterest();
+            if (poi != null) {
+                final MathTransform tr = grid.getGridToCRS(PixelInCell.CELL_CENTER);
+                final GeneralDirectPosition median = new GeneralDirectPosition(tr.getTargetDimensions());
+                tr.transform(poi, 0, median.coordinates, 0, 1);
+                return median;
+            }
+        }
+        return null;
+    }
+
+    /**
+     * Invoked when the target CRS has some dimensions that the source CRS does not have.
      * For example this is invoked during the conversion from (<var>x</var>,
<var>y</var>)
      * coordinates to (<var>x</var>, <var>y</var>, <var>t</var>).
If constant values can
      * be given to the missing dimensions, than those values are returned. Otherwise this
      * method returns {@code null}.
+     *
+     * <p>The returned array has a length equals to the number of dimensions in the
target CRS.
+     * Only coordinates in dimensions without source (<var>t</var> in above example)
will be used.
+     * All other coordinate values will be ignored.</p>
+     *
+     * @see org.apache.sis.referencing.operation.CoordinateOperationContext#getConstantCoordinates()
      */
     @Override
     @SuppressWarnings("ReturnOfCollectionOrArrayField")
@@ -209,7 +399,7 @@ final class CoordinateOperationFinder implements Supplier<double[]>
{
     /**
      * Configures the accuracy hints on the given processor.
      */
-    final void setAccuracy(final ImageProcessor processor) {
+    final void setAccuracyOf(final ImageProcessor processor) {
         final double accuracy = CRS.getLinearAccuracy(operation);
         if (accuracy > 0) {
             Length qm = Quantities.create(accuracy, Units.METRE);
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
index 4ead751..901f12b 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/ResampledGridCoverage.java
@@ -115,7 +115,7 @@ final class ResampledGridCoverage extends GridCoverage {
         }
         processor = processor.clone();
         processor.setFillValues(fillValues);
-        changeOfCRS.setAccuracy(processor);
+        changeOfCRS.setAccuracyOf(processor);
         imageProcessor = GridCoverageProcessor.unique(processor);
     }
 
@@ -232,7 +232,9 @@ final class ResampledGridCoverage extends GridCoverage {
          * transform is missing, we can not continue (we have no way to guess it).
          */
         final MathTransform sourceCornerToCRS = changeOfCRS.gridToCRS(PixelInCell.CELL_CORNER);
+        final MathTransform crsToSourceCorner = changeOfCRS.inverse(sourceCornerToCRS);
         final MathTransform sourceCenterToCRS = changeOfCRS.gridToCRS(PixelInCell.CELL_CENTER);
+        final MathTransform crsToSourceCenter = changeOfCRS.inverse(null);
         /*
          * Compute the transform from target grid to target CRS. This transform may be unspecified,
          * in which case we need to compute a default transform trying to preserve resolution
at the
@@ -398,8 +400,8 @@ final class ResampledGridCoverage extends GridCoverage {
          */
         final MathTransform targetCornerToCRS = resampled.getGridToCRS(PixelInCell.CELL_CORNER);
         return new ResampledGridCoverage(source, resampled,
-                MathTransforms.concatenate(targetCornerToCRS, sourceCornerToCRS.inverse()),
-                MathTransforms.concatenate(targetCenterToCRS, sourceCenterToCRS.inverse()),
+                MathTransforms.concatenate(targetCornerToCRS, crsToSourceCorner),
+                MathTransforms.concatenate(targetCenterToCRS, crsToSourceCenter),
                 changeOfCRS, processor).specialize(isGeometryExplicit);
     }
 
diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java b/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java
index 276729a..c5483a4 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/image/PlanarImage.java
@@ -574,7 +574,7 @@ public abstract class PlanarImage implements RenderedImage {
     @Override
     public String toString() {
         final StringBuilder buffer = new StringBuilder(100).append(Classes.getShortClassName(this))
-                .append('[').append(getWidth()).append(" × ").append(getHeight()).append("
pixels");
+                .append("[(").append(getWidth()).append(" × ").append(getHeight()).append(")
pixels");
         final SampleModel sm = getSampleModel();
         if (sm != null) {
             buffer.append(" × ").append(sm.getNumBands()).append(" bands");
diff --git a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java
b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java
index 64a9a2e..cc30a7b 100644
--- a/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java
+++ b/core/sis-feature/src/test/java/org/apache/sis/coverage/grid/ResampledGridCoverageTest.java
@@ -18,6 +18,7 @@ package org.apache.sis.coverage.grid;
 
 import java.util.Arrays;
 import java.util.Random;
+import java.util.stream.IntStream;
 import java.awt.Color;
 import java.awt.Rectangle;
 import java.awt.image.BufferedImage;
@@ -39,6 +40,7 @@ import org.apache.sis.internal.referencing.Formulas;
 import org.apache.sis.internal.referencing.j2d.AffineTransform2D;
 import org.apache.sis.referencing.CommonCRS;
 import org.apache.sis.referencing.crs.HardCodedCRS;
+import org.apache.sis.referencing.cs.AxesConvention;
 import org.apache.sis.referencing.operation.HardCodedConversions;
 import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.operation.matrix.MatrixSIS;
@@ -559,7 +561,7 @@ public final strictfp class ResampledGridCoverageTest extends TestCase
{
         final GridGeometry targetGrid = new GridGeometry(new GridExtent(4, 4), CELL_CENTER,
gridToCRS, crs);
         final GridCoverage source     = new GridCoverage2D(sourceGrid, null, image);
         final GridCoverage target     = resample(source, targetGrid);
-        assertValuesEqual(target.render(null).getData(), 0, new int[][] {
+        assertValuesEqual(target.render(null), 0, new double[][] {
             {10, 12, 0, 0},
             {16, 14, 0, 0},
             { 0,  0, 0, 0},
@@ -568,6 +570,41 @@ public final strictfp class ResampledGridCoverageTest extends TestCase
{
     }
 
     /**
+     * Tests resampling of an image associated to a coordinate system using the 0 to 360°
range of longitude.
+     * The image crosses the 180° longitude.
+     *
+     * @throws TransformException if some coordinates can not be transformed to the target
grid geometry.
+     */
+    @Test
+    public void testLongitude360() throws TransformException {
+        final int width = 32;
+        final int height = 3;
+        final GridGeometry source = new GridGeometry(
+                new GridExtent(width, height), CELL_CENTER,
+                new AffineTransform2D(1, 0, 0, -1, 164, 0),
+                HardCodedCRS.WGS84.forConvention(AxesConvention.POSITIVE_RANGE));
+        /*
+         * Above grid extent is [164 … 195]° in longitude. The part that exceed 180°
is equivalent to
+         * a [-180 … -165]° range. The extent below requests only a part of it, namely
[-173 … -167]°.
+         * The first pixel of resampled image is the 23th pixel of original image.
+         */
+        final GridGeometry target = new GridGeometry(
+                new GridExtent(7, height), CELL_CENTER,
+                new AffineTransform2D(1, 0, 0, -1, -173, 0),
+                HardCodedCRS.WGS84);
+
+        final BufferedImage image = new BufferedImage(width, height, BufferedImage.TYPE_BYTE_GRAY);
+        image.getRaster().setPixels(0, 0, width, height, IntStream.range(100, 100 + width*height).toArray());
+        final GridCoverage2D coverage = new GridCoverage2D(source, null, image);
+        final GridCoverage resampled = resample(coverage, target);
+        assertValuesEqual(resampled.render(null), 0, new double[][] {
+            {123, 124, 125, 126, 127, 128, 129},
+            {155, 156, 157, 158, 159, 160, 161},
+            {187, 188, 189, 190, 191, 192, 193}
+        });
+    }
+
+    /**
      * Returns an image with only the queries part of the given image.
      * This is an helper tools which can be invoked during debugging
      * session in IDE capable to display images.


Mime
View raw message