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.
|