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: Compute a default GridExtent during resampling operation if it was not specified by user.
Date Sun, 29 Mar 2020 21:22:45 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 e2f3f2a  Compute a default GridExtent during resampling operation if it was not specified
by user.
e2f3f2a is described below

commit e2f3f2aecc2b6f3f6aeca4726d5dd1d8ddedce80
Author: Martin Desruisseaux <martin.desruisseaux@geomatys.com>
AuthorDate: Sun Mar 29 23:21:31 2020 +0200

    Compute a default GridExtent during resampling operation if it was not specified by user.
---
 .../apache/sis/coverage/grid/GridCoverage2D.java   |   2 +-
 .../sis/coverage/grid/ResampledGridCoverage.java   | 172 ++++++++++++++++++---
 .../org/apache/sis/coverage/grid/package-info.java |   1 +
 .../java/org/apache/sis/image/package-info.java    |   1 +
 .../sis/referencing/operation/matrix/Matrix1.java  |   6 +-
 .../sis/referencing/operation/matrix/Matrix2.java  |   8 +-
 .../referencing/operation/matrix/MatrixSIS.java    |   8 +-
 7 files changed, 170 insertions(+), 28 deletions(-)

diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java
b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java
index 4889928..a89bdb7 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridCoverage2D.java
@@ -102,7 +102,7 @@ public class GridCoverage2D extends GridCoverage {
     /**
      * Minimal number of dimension for this coverage.
      */
-    private static final int MIN_DIMENSION = 2;
+    static final int MIN_DIMENSION = 2;
 
     /**
      * The sample values stored as a {@code RenderedImage}.
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 68a1f40..160bce5 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
@@ -17,6 +17,7 @@
 package org.apache.sis.coverage.grid;
 
 import java.util.List;
+import java.util.Arrays;
 import java.util.Optional;
 import java.awt.Rectangle;
 import java.awt.image.RenderedImage;
@@ -31,15 +32,16 @@ import org.apache.sis.geometry.Envelopes;
 import org.apache.sis.image.Interpolation;
 import org.apache.sis.image.ResampledImage;
 import org.apache.sis.coverage.SampleDimension;
+import org.apache.sis.geometry.GeneralEnvelope;
+import org.apache.sis.internal.util.DoubleDouble;
 import org.apache.sis.internal.coverage.j2d.ImageUtilities;
+import org.apache.sis.internal.referencing.DirectPositionView;
 import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox;
 import org.apache.sis.referencing.operation.transform.MathTransforms;
+import org.apache.sis.referencing.operation.matrix.MatrixSIS;
+import org.apache.sis.referencing.operation.matrix.Matrices;
 import org.apache.sis.referencing.CRS;
 
-import static java.lang.Math.toIntExact;
-import static java.lang.Math.negateExact;
-import static java.lang.Math.subtractExact;
-
 
 /**
  * A multi-dimensional grid coverage where each two-dimensional slice is the resampling
@@ -47,6 +49,7 @@ import static java.lang.Math.subtractExact;
  * stored in a {@link GridCoverage2D}.
  *
  * @author  Martin Desruisseaux (Geomatys)
+ * @author  Johann Sorel (Geomatys)
  * @version 1.1
  * @since   1.1
  * @module
@@ -124,7 +127,7 @@ final class ResampledGridCoverage extends GridCoverage {
      */
     private GridCoverage specialize() {
         final GridExtent extent = gridGeometry.getExtent();
-        if (extent.getSubDimension() > BIDIMENSIONAL) {
+        if (extent.getDimension() < GridCoverage2D.MIN_DIMENSION || extent.getSubDimension()
> BIDIMENSIONAL) {
             return this;
         }
         return new GridCoverage2D(source, gridGeometry, extent, render(null), xDimension,
yDimension);
@@ -142,6 +145,9 @@ final class ResampledGridCoverage extends GridCoverage {
     static GridCoverage create(final GridCoverage source, GridGeometry target, final Interpolation
interpolation)
             throws FactoryException, TransformException
     {
+        final CoordinateReferenceSystem sourceCRS = source.getCoordinateReferenceSystem();
+        final CoordinateReferenceSystem targetCRS = target.isDefined(GridGeometry.CRS) ?
+                                                    target.getCoordinateReferenceSystem()
: sourceCRS;
         /*
          * Get the coordinate operation from source CRS to target CRS. It may be the identity
operation,
          * or null only if there is not enough information for determining the operation.
We try to take
@@ -153,9 +159,6 @@ final class ResampledGridCoverage extends GridCoverage {
             changeOfCRS = Envelopes.findOperation(sourceGG.getEnvelope(), target.getEnvelope());
         }
         if (changeOfCRS == null) try {
-            final CoordinateReferenceSystem sourceCRS = source.getCoordinateReferenceSystem();
-            final CoordinateReferenceSystem targetCRS = target.isDefined(GridGeometry.CRS)
?
-                                                        target.getCoordinateReferenceSystem()
: sourceCRS;
             DefaultGeographicBoundingBox areaOfInterest = null;
             if (sourceGG.isDefined(GridGeometry.ENVELOPE)) {
                 areaOfInterest = new DefaultGeographicBoundingBox();
@@ -171,21 +174,124 @@ final class ResampledGridCoverage extends GridCoverage {
          * 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 toTarget = sourceGG.getGridToCRS(PixelInCell.CELL_CENTER);
+        MathTransform sourceGridToCRS = sourceGG.getGridToCRS(PixelInCell.CELL_CENTER);
         if (changeOfCRS != null) {
-            toTarget = MathTransforms.concatenate(toTarget, changeOfCRS.getMathTransform());
+            sourceGridToCRS = MathTransforms.concatenate(sourceGridToCRS, changeOfCRS.getMathTransform());
         }
         /*
          * 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
          * point of interest.
          */
-        MathTransform toSource = null;
+        GridExtent targetExtent = target.isDefined(GridGeometry.EXTENT) ? target.getExtent()
: null;
+        final MathTransform targetGridToCRS;
         if (target.isDefined(GridGeometry.GRID_TO_CRS)) {
-            toSource = target.getGridToCRS(PixelInCell.CELL_CENTER);
+            targetGridToCRS = target.getGridToCRS(PixelInCell.CELL_CENTER);
+            if (targetExtent == null) {
+                targetExtent = targetExtent(sourceGG, changeOfCRS, targetGridToCRS.inverse());
+            }
         } else {
-            throw new UnsupportedOperationException("Automatic computation of gridToCRS not
yet implemented.");
-            // TODO: complete the target GridGeometry.
+            /*
+             * The first column below gives the displacement in target CRS when moving is
the source grid by one cell
+             * toward right, and the second column gives the displacement when moving one
cell toward up (positive y).
+             * More columns may exist in 3D, 4D, etc. cases.
+             */
+            final MatrixSIS vectors = MatrixSIS.castOrCopy(sourceGridToCRS.derivative(
+                    new DirectPositionView.Double(sourceGG.getExtent().getPointOfInterest())));
+            /*
+             * We will retain only the magnitudes of those vectors, in order to have directions
parallel with target
+             * grid axes. There is one magnitude value for each target CRS dimension. If
there is more target grid
+             * dimensions than magnitude values (unusual, but not forbidden), some grid dimensions
will be ignored
+             * provided that their size is 1 (otherwise a SubspaceNotSpecifiedException is
thrown).
+             */
+            final double[]  magnitudes = vectors.normalizeColumns();          // Length is
dimension  of source grid.
+            final int       crsDim     = vectors.getNumRow();                 // Number of
dimensions of target CRS.
+            final int       gridDim    = target.getDimension();               // Number of
dimensions of target grid.
+            final int       mappedDim  = Math.min(magnitudes.length, Math.min(crsDim, gridDim));
+            final MatrixSIS crsToGrid  = Matrices.createZero(gridDim + 1, crsDim + 1);
+            final int[]     dimSelect  = (gridDim > crsDim && targetExtent !=
null) ?
+                                         targetExtent.getSubspaceDimensions(crsDim) : null;
+            /*
+             * The goal below is to build a target "gridToCRS" which perform the same axis
swapping and flipping than
+             * the source "gridToCRS". For example if the source "gridToCRS" was flipping
y axis, then we want target
+             * "gridToCRS" to also flips that axis, unless the transformation from "source
CRS" to "target CRS" flips
+             * that axis, in which case the result in target "gridToCRS" would be to not
flip again:
+             *
+             *     (source gridToCRS)    →    (source CRS to target CRS)    →    (target
gridToCRS)⁻¹
+             *         flip y axis             no axis direction change              flip
y axis
+             *   or    flip y axis                  flip y axis                   no additional
flip
+             *
+             * For each column, the row index of the greatest absolute value is taken as
the target dimension where
+             * to set vector magnitude in target "crsToGrid" matrix. That way, if transformation
from source CRS to
+             * target CRS does not flip or swap axis, target `gridToCRS` matrix should looks
like source `gridToCRS`
+             * matrix (i.e. zero and non-zero coefficients at the same places). If two vectors
have their greatest
+             * value on the same row, the largest value win (since the `vectors` matrix has
been normalized to unit
+             * vectors, values in different columns are comparable).
+             */
+            for (;;) {
+                double max = 0;
+                int tgDim = -1;                         // Grid dimension of maximal value.
+                int tcDim = -1;                         // CRS dimension of maximal value.
+                for (int i=0; i<mappedDim; i++) {
+                    // `ci` differs from `i` only if the source grid has "too much" dimensions.
+                    final int ci = (dimSelect != null) ? dimSelect[i] : i;
+                    for (int j=0; j<crsDim; j++) {
+                        final double m = Math.abs(vectors.getElement(j,ci));
+                        if (m > max) {
+                            max   = m;
+                            tcDim = j;
+                            tgDim = ci;
+                        }
+                    }
+                }
+                if (tgDim < 0) break;                   // No more non-zero value in the
`vectors` matrix.
+                for (int j=0; j<crsDim; j++) {
+                    vectors.setElement(j, tgDim, 0);    // For preventing this column to
be selected again.
+                }
+                crsToGrid.setElement(tgDim, tcDim, magnitudes[tgDim]);
+            }
+            crsToGrid.setElement(gridDim, crsDim, 1);
+            /*
+             * At this point we got a first estimation of "target CRS to grid" transform,
without translation terms.
+             * Apply the complete transform chain on source extent; this will give us a tentative
target extent.
+             * This tentative extent will be compared with desired target.
+             */
+            final GridExtent tentative = targetExtent(sourceGG, changeOfCRS, MathTransforms.linear(crsToGrid));
+            if (targetExtent == null) {
+                // Create an extent of same size but with lower coordinates set to 0.
+                final long[] coordinates = new long[gridDim * 2];
+                for (int i=0; i<gridDim; i++) {
+                    coordinates[i + gridDim] = tentative.getSize(i);
+                }
+                targetExtent = new GridExtent(tentative, coordinates);
+            }
+            /*
+             * At this point we have the desired target extent and the extent that we actually
got by applying
+             * full "source to target" transform. Compute the scale and offset differences
between target and
+             * actual extents, then adjust matrix coefficients for compensating those differences.
+             */
+            final DoubleDouble scale  = new DoubleDouble();
+            final DoubleDouble offset = new DoubleDouble();
+            final DoubleDouble tmp    = new DoubleDouble();
+            for (int j=0; j<gridDim; j++) {
+                tmp.set(targetExtent.getSize(j));
+                scale.set(tentative.getSize(j));
+                scale.inverseDivide(tmp);
+
+                tmp.set(targetExtent.getLow(j));
+                offset.set(-tentative.getLow(j));
+                offset.multiply(scale);
+                offset.add(tmp);
+                crsToGrid.convertAfter(j, scale, offset);
+            }
+            targetGridToCRS = MathTransforms.linear(crsToGrid.inverse());
+        }
+        /*
+         * At this point all target grid geometry components are non-null.
+         * Build the final target GridGeometry if any components were missing.
+         */
+        if (!target.isDefined(GridGeometry.EXTENT | GridGeometry.GRID_TO_CRS | GridGeometry.CRS))
{
+            target = new GridGeometry(targetExtent, PixelInCell.CELL_CENTER, targetGridToCRS,
targetCRS);
         }
         if (sourceGG.equals(target)) {
             return source;
@@ -193,11 +299,35 @@ final class ResampledGridCoverage extends GridCoverage {
         /*
          * Complete the "target to source" transform.
          */
-        toSource = MathTransforms.concatenate(toSource, toTarget.inverse());
+        final MathTransform toSource = MathTransforms.concatenate(targetGridToCRS, sourceGridToCRS.inverse());
         return new ResampledGridCoverage(source, target, toSource, interpolation).specialize();
     }
 
     /**
+     * Computes a target grid extent by transforming the source grid extent. This method
expects a transform
+     * mapping cell centers, but will adjust for getting a result as if the transform was
mapping cell corners.
+     *
+     * @param  source       the source grid geometry.
+     * @param  changeOfCRS  the coordinate operation from source CRS to target CRS, or {@code
null}.
+     * @param  crsToGrid    the transform from target CRS to target grid, mapping cell <em>centers</em>.
+     * @return target grid extent.
+     */
+    private static GridExtent targetExtent(final GridGeometry source, final CoordinateOperation
changeOfCRS,
+            final MathTransform crsToGrid) throws TransformException
+    {
+        MathTransform extentMapper = source.getGridToCRS(PixelInCell.CELL_CORNER);
+        if (changeOfCRS != null) {
+            extentMapper = MathTransforms.concatenate(extentMapper, changeOfCRS.getMathTransform());
+        }
+        extentMapper = MathTransforms.concatenate(extentMapper, crsToGrid);
+        final GeneralEnvelope bounds = source.getExtent().toCRS(extentMapper, extentMapper,
null);
+        final double[] vector = new double[bounds.getDimension()];
+        Arrays.fill(vector, 0.5);
+        bounds.translate(vector);       // Convert cell centers to cell corners.
+        return new GridExtent(bounds, GridRoundingMode.NEAREST, null, null, null);
+    }
+
+    /**
      * Returns a two-dimensional slice of resampled grid data as a rendered image.
      */
     @Override
@@ -208,20 +338,20 @@ final class ResampledGridCoverage extends GridCoverage {
         final RenderedImage image        = source.render(null);           // TODO: compute
slice.
         final GridExtent    sourceExtent = source.getGridGeometry().getExtent();
         final GridExtent    targetExtent = gridGeometry.getExtent();
-        final Rectangle     bounds       = new Rectangle(toIntExact(targetExtent.getSize(xDimension)),
-                                                         toIntExact(targetExtent.getSize(yDimension)));
+        final Rectangle     bounds       = new Rectangle(Math.toIntExact(targetExtent.getSize(xDimension)),
+                                                         Math.toIntExact(targetExtent.getSize(yDimension)));
         /*
          * `this.toSource` is a transform from source cell coordinates to target cell coordinates.
          * We need a transform from source pixel coordinates to target pixel coordinates
(in images).
          * An offset may exist between cell coordinates and pixel coordinates.
          */
         final MathTransform pixelsToTransform = MathTransforms.translation(
-                subtractExact(sourceExtent.getLow(xDimension), image.getMinX()),
-                subtractExact(sourceExtent.getLow(yDimension), image.getMinY()));
+                targetExtent.getLow(xDimension),
+                targetExtent.getLow(yDimension));
 
         final MathTransform transformToPixels = MathTransforms.translation(
-                negateExact(targetExtent.getLow(xDimension)),
-                negateExact(targetExtent.getLow(yDimension)));
+                Math.subtractExact(image.getMinX(), sourceExtent.getLow(xDimension)),
+                Math.subtractExact(image.getMinY(), sourceExtent.getLow(yDimension)));
 
         final MathTransform toImage = MathTransforms.concatenate(pixelsToTransform, toSource,
transformToPixels);
         /*
diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/package-info.java
b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/package-info.java
index 29fb023..7e382c0 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/package-info.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/package-info.java
@@ -21,6 +21,7 @@
  * In the three-dimensional case, the cells are voxels.
  *
  * @author  Martin Desruisseaux (Geomatys)
+ * @author  Johann Sorel (Geomatys)
  * @version 1.1
  * @since   1.0
  * @module
diff --git a/core/sis-feature/src/main/java/org/apache/sis/image/package-info.java b/core/sis-feature/src/main/java/org/apache/sis/image/package-info.java
index 142faab..04fc7a6 100644
--- a/core/sis-feature/src/main/java/org/apache/sis/image/package-info.java
+++ b/core/sis-feature/src/main/java/org/apache/sis/image/package-info.java
@@ -34,6 +34,7 @@
  *
  * @author  Rémi Maréchal (Geomatys)
  * @author  Martin Desruisseaux (Geomatys)
+ * @author  Johann Sorel (Geomatys)
  * @version 1.1
  * @since   1.0
  * @module
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
index ff93401..49bdf2b 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix1.java
@@ -252,10 +252,14 @@ public class Matrix1 extends MatrixSIS {
     /**
      * Normalizes all columns in-place.
      * For a 1×1 matrix, this method just sets unconditionally the {@link #m00} value to
1.
+     *
+     * @return the magnitude of the column, which is the absolute value of {@link #m00}.
      */
     @Override
-    public void normalizeColumns() {
+    public double[] normalizeColumns() {
+        final double[] magnitudes = new double[] {Math.abs(m00)};
         m00 = 1;
+        return magnitudes;
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
index c3f455a..a717530 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/Matrix2.java
@@ -276,10 +276,12 @@ public class Matrix2 extends MatrixSIS {
      * {@inheritDoc}
      */
     @Override
-    public void normalizeColumns() {
+    public double[] normalizeColumns() {
+        final double[] magnitudes = new double[2];
         double m;
-        m = Math.hypot(m00, m10); m00 /= m; m10 /= m;
-        m = Math.hypot(m01, m11); m01 /= m; m11 /= m;
+        magnitudes[0] = m = Math.hypot(m00, m10); m00 /= m; m10 /= m;
+        magnitudes[1] = m = Math.hypot(m01, m11); m01 /= m; m11 /= m;
+        return magnitudes;
     }
 
     /**
diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
index a39ca3c..a34112a 100644
--- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
+++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/matrix/MatrixSIS.java
@@ -55,7 +55,7 @@ import org.apache.sis.util.resources.Errors;
  * </ul>
  *
  * @author  Martin Desruisseaux (IRD, Geomatys)
- * @version 1.0
+ * @version 1.1
  *
  * @see Matrices
  *
@@ -422,11 +422,13 @@ public abstract class MatrixSIS implements Matrix, LenientComparable,
Cloneable,
      * coordinate in the source space is increased by one. Invoking this method turns those
vectors
      * into unitary vectors, which is useful for forming the basis of a new coordinate system.</p>
      *
+     * @return the magnitude for each column. The length of this array is the number of columns.
      * @throws UnsupportedOperationException if this matrix is unmodifiable.
      */
-    public void normalizeColumns() {
+    public double[] normalizeColumns() {
         final int numRow = getNumRow();
         final int numCol = getNumCol();
+        final double[] magnitudes = new double[numCol];
         final DoubleDouble sum = new DoubleDouble();
         final DoubleDouble dot = new DoubleDouble();
         final DoubleDouble tmp = new DoubleDouble();
@@ -444,7 +446,9 @@ public abstract class MatrixSIS implements Matrix, LenientComparable,
Cloneable,
                 dot.inverseDivide(tmp);
                 set(j, i, dot);
             }
+            magnitudes[i] = sum.doubleValue();
         }
+        return magnitudes;
     }
 
     /**


Mime
View raw message