This is an automated email from the ASF dual-hosted git repository. desruisseaux pushed a commit to branch geoapi-4.0 in repository https://gitbox.apache.org/repos/asf/sis.git
commit 6caf39fbe6d503936073a0c3e8052ea8088d68dc Author: Martin Desruisseaux <[email protected]> AuthorDate: Tue Mar 23 17:02:24 2021 +0100 Refactor `GridGeometry.linearize(…)` method in a separated `DomainLinearizer` class. The intent is to make linearization more parametrable, e.g. with a scale factor. --- .../apache/sis/coverage/grid/DomainLinearizer.java | 205 +++++++++++++++++++++ .../org/apache/sis/coverage/grid/GridGeometry.java | 79 +++----- .../operation/builder/LinearTransformBuilder.java | 2 +- 3 files changed, 227 insertions(+), 59 deletions(-) diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DomainLinearizer.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DomainLinearizer.java new file mode 100644 index 0000000..ce90603 --- /dev/null +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/DomainLinearizer.java @@ -0,0 +1,205 @@ +/* + * Licensed to the Apache Software Foundation (ASF) under one or more + * contributor license agreements. See the NOTICE file distributed with + * this work for additional information regarding copyright ownership. + * The ASF licenses this file to You under the Apache License, Version 2.0 + * (the "License"); you may not use this file except in compliance with + * the License. You may obtain a copy of the License at + * + * http://www.apache.org/licenses/LICENSE-2.0 + * + * Unless required by applicable law or agreed to in writing, software + * distributed under the License is distributed on an "AS IS" BASIS, + * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. + * See the License for the specific language governing permissions and + * limitations under the License. + */ +package org.apache.sis.coverage.grid; + +import java.util.Arrays; +import org.opengis.util.FactoryException; +import org.opengis.referencing.datum.PixelInCell; +import org.opengis.referencing.operation.MathTransform; +import org.opengis.referencing.operation.TransformException; +import org.apache.sis.referencing.operation.builder.LinearTransformBuilder; +import org.apache.sis.referencing.operation.transform.LinearTransform; +import org.apache.sis.referencing.operation.transform.MathTransforms; +import org.apache.sis.geometry.GeneralEnvelope; +import org.apache.sis.geometry.Envelopes; +import org.apache.sis.util.ArgumentChecks; + + +/** + * Method for replacing a non-linear "grid to CRS" conversion by a linear conversion (affine transform). + * The {@code GridGeometry} class allows non-linear {@linkplain GridGeometry#getGridToCRS "grid to CRS"} + * conversions, but some {@code GridGeometry} usages are restricted to linear (affine) conversions. + * The {@code DomainLinearizer} class encapsulates the method used for replacing non-linear conversions + * by a linear approximation. + * + * <p>The same instance can be reused by invoking {@code apply(…)} methods for many {@link GridCoverage} + * or {@link GridGeometry} instances.</p> + * + * <h2>Limitations</h2> + * Current implementation is designed for two-dimensional grid geometries. + * Support for higher dimensions is not guaranteed. + * + * @author Martin Desruisseaux (Geomatys) + * @version 1.1 + * + * @see LinearTransformBuilder#approximate(MathTransform, Envelope) + * + * @since 1.1 + * @module + */ +public class DomainLinearizer { + /** + * Whether to force lower grid coordinates to (0,0,…). + * + * @see #getGridStartsAtZero() + */ + private boolean gridStartsAtZero; + + /** + * Scale factor applied on grid coordinates, or 1 if none. + */ + private double scale = 1; + + /** + * Creates a new linearizer. + */ + public DomainLinearizer() { + } + + /** + * Returns whether {@code GridExtent} should have their lower grid coordinates set to zero. + * If {@code true}, then {@code DomainLinearizer} will opportunistically apply translation + * on the {@linkplain GridGeometry#getGridToCRS "grid to CRS"} conversion in such a way that + * {@link GridExtent#getLow()} is 0 for all dimensions. + * + * @return whether to force lower grid coordinates to (0,0,…). + * + * @see GridExtent#startsAtZero() + */ + public boolean getGridStartsAtZero() { + return gridStartsAtZero; + } + + /** + * Sets whether {@code GridExtent} should have their lower grid coordinates set to zero. + * The default value is {@code false}. + * + * @param force whether to force lower grid coordinates to (0,0,…). + */ + public void setGridStartsAtZero(final boolean force) { + gridStartsAtZero = force; + } + + /** + * Returns the scale factor applied on coordinates in all dimensions. + * + * @return scale factor applied on coordinates in all dimensions, or 1 if none. + */ + public double getUniformScale() { + return scale; + } + + /** + * Sets the scale factor to apply on coordinates in all dimensions. + * The default value is 1. + * + * @param factor scale factor applied on coordinates in all dimensions, or 1 if none. + */ + public void setUniformScale(final double factor) { + ArgumentChecks.ensureStrictlyPositive("factor", factor); + scale = factor; + } + + /** + * Returns a grid coverage with a linear approximation of the <cite>grid to CRS</cite> conversion. + * The linear approximation is computed by {@link #apply(GridGeometry)}. If the <cite>grid to CRS</cite> + * conversion of the given coverage is already linear, then this method returns {@code coverage}. + * + * @param coverage the grid coverage in which to make the <cite>grid to CRS</cite> conversion linear. + * @return a grid coverage with a linear approximation of the <cite>grid to CRS</cite> conversion. + * @throws TransformException if some cell coordinates can not be computed. + */ + public GridCoverage apply(final GridCoverage coverage) throws TransformException { + ArgumentChecks.ensureNonNull("coverage", coverage); + final GridGeometry gg = coverage.getGridGeometry(); + final GridGeometry linearized = apply(gg); + if (gg.equals(linearized)) { + return coverage; + } + final GridCoverageProcessor processor = new GridCoverageProcessor(); + return processor.resample(coverage, linearized); + } + + /** + * Creates a grid geometry with a linear approximation of the <cite>grid to CRS</cite> conversion. + * The approximation is computed by <cite>Least Mean Squares</cite> method: the affine transform + * coefficients are chosen in way making the average value of (<var>position</var> − <var>linear + * approximation of position</var>)² as small as possible for all cells in given grid geometry. + * If the <cite>grid to CRS</cite> conversion of the given grid geometry is already linear, + * then this method returns {@code gg}. + * + * @param gg the grid geometry in which to make the <cite>grid to CRS</cite> conversion linear. + * @return a grid geometry with a linear approximation of the <cite>grid to CRS</cite> conversion. + * @throws TransformException if some cell coordinates can not be computed. + */ + public GridGeometry apply(final GridGeometry gg) throws TransformException { + ArgumentChecks.ensureNonNull("gg", gg); + if (gg.nonLinears != 0) try { + MathTransform gridToCRS = gg.requireGridToCRS(true); + GeneralEnvelope domain = gg.extent.toCRS(null, null, null); + MathTransform approximate = modify(LinearTransformBuilder.approximate(gridToCRS, domain)); + MathTransform gridToGrid = MathTransforms.concatenate(gridToCRS, approximate.inverse()); + domain = Envelopes.transform(gridToGrid, domain); + final int dimension = domain.getDimension(); + final long[] coordinates = new long[dimension * 2]; + final double[] shift = new double[dimension]; + for (int i=0; i<dimension; i++) { + long low = Math.round(domain.getMinimum(i)); + long high = Math.round(domain.getMaximum(i)); + high = Math.max(low, Math.decrementExact(high)); + if (gridStartsAtZero) { + high = Math.subtractExact(high, low); + shift[i] = low; + } else { + coordinates[i] = low; + } + coordinates[i + dimension] = high; + } + approximate = MathTransforms.concatenate(MathTransforms.translation(shift), approximate); + if (!approximate.equals(gridToCRS)) { + return new GridGeometry(new GridExtent(gg.extent, coordinates), PixelInCell.CELL_CENTER, + approximate, gg.envelope.getCoordinateReferenceSystem()); + } + } catch (FactoryException e) { + final Throwable cause = e.getCause(); + if (cause instanceof TransformException) { + throw (TransformException) cause; + } + throw new TransformException(e); + } + return gg; + } + + /** + * Callback for custom modification of linear approximation. This method is invoked by {@link #apply(GridGeometry)} + * after a linear "grid to CRS" approximation has been computed by the <cite>Least Mean Squares</cite> method. + * Subclasses can override this method for example in order to scale the conversion by some arbitrary factor. + * + * <div class="note"><b>Tip:</b> + * scales (if desired) should be applied <em>before</em> {@code gridToCRS}, i.e. on grid coordinates. + * Scales applied after the transform (i.e. on "real world" coordinates) may give unexpected results + * if the conversion contains a rotation.</div> + * + * @param gridToCRS an approximation of the "grid to CRS" conversion computed by {@code DomainLinearizer}. + * @return the approximation to use for creating a new {@link GridGeometry}. Should be linear. + */ + private MathTransform modify(final LinearTransform gridToCRS) { + final double[] factors = new double[gridToCRS.getTargetDimensions()]; + Arrays.fill(factors, scale); + return MathTransforms.concatenate(MathTransforms.scale(factors), gridToCRS); + } +} diff --git a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java index 768ab7e..cc50d6b 100644 --- a/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java +++ b/core/sis-feature/src/main/java/org/apache/sis/coverage/grid/GridGeometry.java @@ -51,7 +51,6 @@ import org.apache.sis.referencing.operation.matrix.Matrices; import org.apache.sis.referencing.operation.matrix.MatrixSIS; import org.apache.sis.referencing.operation.transform.MathTransforms; import org.apache.sis.referencing.operation.transform.PassThroughTransform; -import org.apache.sis.referencing.operation.builder.LinearTransformBuilder; import org.apache.sis.internal.referencing.DirectPositionView; import org.apache.sis.internal.referencing.TemporalAccessor; import org.apache.sis.internal.referencing.AxisDirections; @@ -100,15 +99,30 @@ import org.apache.sis.xml.NilReason; * <li>An {@linkplain #isConversionLinear indication of whether conversion for some axes is linear or not}.</li> * </ul> * - * The first three properties should be mandatory, but are allowed to be temporarily absent during - * grid coverage construction. Temporarily absent properties are allowed because they may be inferred - * from a wider context. For example a grid geometry know nothing about {@link RenderedImage}, - * but {@code GridCoverage2D} does and may use that information for providing a missing grid extent. + * The first three properties should be mandatory, + * but are allowed to be temporarily absent during grid coverage construction. + * Temporarily absent properties are allowed because they may be inferred from a wider context. + * For example a {@code GridGeometry} knows nothing about {@link RenderedImage}, + * but {@code GridCoverage2D} has this information and may use it for providing a missing grid extent. * By default, any request for an undefined property will throw an {@link IncompleteGridGeometryException}. * In order to check if a property is defined, use {@link #isDefined(int)}. * - * <p>{@code GridGeometry} instances are immutable and thread-safe. - * The same instance can be shared by different {@link GridCoverage} instances.</p> + * <h2>Non-linear referencing</h2> + * A key property is the {@linkplain #getGridToCRS(PixelInCell) "grid to CRS"} conversion, + * which defines how to map pixel coordinates to "real world" coordinates such as latitudes and longitudes. + * This relationship is often linear (an affine transform), but {@linkplain #isConversionLinear not necessarily}; + * {@code GridGeometry} accepts non-linear conversions as well. Non-linear conversions may occur with images + * using {@linkplain org.apache.sis.referencing.operation.builder.LocalizationGridBuilder localization grids}, + * but non-linear conversions should not be used for expressing map projections (projections should be specified + * in the {@linkplain #getCoordinateReferenceSystem() Coordinate Reference System} (CRS) instead). + * + * <p>Some applications can not handle non-linear "grid to CRS" conversions. + * For example encoding an image in a GeoTIFF file is much simpler if the "grid to CRS" conversion is linear. + * The {@link DomainLinearizer} class can be used for replacing non-linear conversions by linear approximations.</p> + * + * <h2>Multi-threading</h2> + * {@code GridGeometry} instances are immutable and thread-safe. + * The same instance can be shared by different {@link GridCoverage} instances. * * @author Martin Desruisseaux (IRD, Geomatys) * @version 1.1 @@ -1291,57 +1305,6 @@ public class GridGeometry implements LenientComparable, Serializable { } /** - * Creates a grid geometry with a linear approximation of the <cite>grid to CRS</cite> transform. - * The approximation is computed by <cite>Least Mean Squares</cite> method: the affine transform - * coefficients are chosen in way making the average value of (<var>position</var> − <var>linear - * approximation of position</var>)² as small as possible for all cells in this grid geometry. - * If the <cite>grid to CRS</cite> transform of this grid geometry is already linear, then this - * method returns {@code this}. - * - * @param forceLowerToZero whether to force lower grid coordinates to (0,0,…). - * @return a grid geometry with a linear approximation of the <cite>grid to CRS</cite> transform. - * @throws TransformException if some cell coordinates can not be computed. - * - * @see LinearTransformBuilder#approximate(MathTransform, Envelope) - * @see GridExtent#startsAtZero() - */ - public GridGeometry linearize(final boolean forceLowerToZero) throws TransformException { - if (nonLinears != 0 && extent != null) try { - GeneralEnvelope domain = extent.toCRS(null, null, null); - MathTransform approximate = LinearTransformBuilder.approximate(gridToCRS, domain); - MathTransform gridToGrid = MathTransforms.concatenate(gridToCRS, approximate.inverse()); - domain = Envelopes.transform(gridToGrid, domain); - final int dimension = domain.getDimension(); - final long[] coordinates = new long[dimension * 2]; - final double[] shift = new double[dimension]; - for (int i=0; i<dimension; i++) { - long low = Math.round(domain.getMinimum(i)); - long high = Math.round(domain.getMaximum(i)); - high = Math.max(low, Math.decrementExact(high)); - if (forceLowerToZero) { - high = Math.subtractExact(high, low); - shift[i] = low; - } else { - coordinates[i] = low; - } - coordinates[i + dimension] = high; - } - approximate = MathTransforms.concatenate(MathTransforms.translation(shift), approximate); - if (!approximate.equals(gridToCRS)) { - return new GridGeometry(new GridExtent(extent, coordinates), PixelInCell.CELL_CENTER, - approximate, envelope.getCoordinateReferenceSystem()); - } - } catch (FactoryException e) { - final Throwable cause = e.getCause(); - if (cause instanceof TransformException) { - throw (TransformException) cause; - } - throw new TransformException(e); - } - return this; - } - - /** * Creates a transform from cell coordinates in this grid to cell coordinates in the given grid. * The returned transform handles change of Coordinate Reference System and wraparound axes * (e.g. longitude axis crossing the ±180° meridian) if applicable. diff --git a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java index fc78b3e..aeca54a 100644 --- a/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java +++ b/core/sis-referencing/src/main/java/org/apache/sis/referencing/operation/builder/LinearTransformBuilder.java @@ -260,7 +260,7 @@ public class LinearTransformBuilder extends TransformBuilder { * @throws FactoryException if the transform approximation can not be computed. * * @see #setControlPoints(MathTransform) - * @see org.apache.sis.coverage.grid.GridGeometry#linearize(boolean) + * @see org.apache.sis.coverage.grid.DomainLinearizer * * @since 1.1 */
