| /* |
| * 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.referencing.operation.builder; |
| |
| import java.util.Map; |
| import java.util.Locale; |
| import java.util.Optional; |
| import java.io.IOException; |
| import java.io.UncheckedIOException; |
| import org.opengis.util.FactoryException; |
| import org.opengis.geometry.Envelope; |
| import org.opengis.geometry.MismatchedDimensionException; |
| import org.opengis.referencing.operation.Matrix; |
| import org.opengis.referencing.operation.MathTransform; |
| import org.opengis.referencing.operation.MathTransformFactory; |
| import org.opengis.referencing.operation.TransformException; |
| import org.opengis.referencing.operation.NoninvertibleTransformException; |
| import org.apache.sis.referencing.factory.InvalidGeodeticParameterException; |
| import org.apache.sis.referencing.operation.transform.InterpolatedTransform; |
| import org.apache.sis.referencing.operation.transform.LinearTransform; |
| import org.apache.sis.referencing.operation.transform.MathTransforms; |
| import org.apache.sis.referencing.operation.matrix.Matrix3; |
| import org.apache.sis.referencing.datum.DatumShiftGrid; |
| import org.apache.sis.internal.referencing.Resources; |
| import org.apache.sis.geometry.GeneralEnvelope; |
| import org.apache.sis.geometry.Envelopes; |
| import org.apache.sis.internal.util.Numerics; |
| import org.apache.sis.internal.util.Strings; |
| import org.apache.sis.util.resources.Vocabulary; |
| import org.apache.sis.util.resources.Errors; |
| import org.apache.sis.util.ArgumentChecks; |
| import org.apache.sis.measure.NumberRange; |
| import org.apache.sis.math.MathFunctions; |
| import org.apache.sis.math.Statistics; |
| import org.apache.sis.math.StatisticsFormat; |
| import org.apache.sis.math.Vector; |
| |
| import static org.apache.sis.referencing.operation.builder.ResidualGrid.SOURCE_DIMENSION; |
| |
| |
| /** |
| * Creates an "almost linear" transform mapping the given source points to the given target points. |
| * The transform is backed by a <cite>grid of localization</cite>, a two-dimensional array of coordinate points. |
| * Grid size is {@code width} × {@code height} and input coordinates are (<var>i</var>,<var>j</var>) indices in the grid, |
| * where <var>i</var> must be in the [0…{@code width}-1] range and <var>j</var> in the [0…{@code height}-1] range inclusive. |
| * Output coordinates are the values stored in the grid of localization at the specified index. |
| * After a {@code LocalizationGridBuilder} instance has been fully populated (i.e. real world coordinates have been |
| * specified for all grid cells), a transformation from grid coordinates to "real world" coordinates can be obtained |
| * with the {@link #create(MathTransformFactory)} method. If this transform is close enough to an affine transform, |
| * then an instance of {@link LinearTransform} is returned. |
| * Otherwise, a transform backed by the localization grid is returned. |
| * |
| * <p>This builder performs the following steps:</p> |
| * <ol> |
| * <li>Compute a linear approximation of the transformation using {@link LinearTransformBuilder}.</li> |
| * <li>Compute {@link DatumShiftGrid} with the residuals.</li> |
| * <li>Create a {@link InterpolatedTransform} with the above shift grid.</li> |
| * <li>If a {@linkplain LinearTransformBuilder#linearizer() linearizer has been applied}, |
| * concatenate the inverse transform of that linearizer.</li> |
| * </ol> |
| * |
| * Builders are not thread-safe. Builders can be used only once; |
| * points can not be added or modified after {@link #create(MathTransformFactory)} has been invoked. |
| * |
| * <h2>Linearizers</h2> |
| * If the localization grid is not close enough to a linear transform, {@link InterpolatedTransform} may not converge. |
| * To improve the speed and reliability of the transform, a non-linear step can be {@linkplain #addLinearizers specified}. |
| * Many candidates can be specified in case the exact form of that non-linear step is unknown; |
| * {@code LocalizationGridBuilder} will select the non-linear step that provides the best improvement, if any. |
| * See the <cite>Linearizers</cite> section in {@link LinearTransformBuilder} for more discussion. |
| * |
| * @author Martin Desruisseaux (Geomatys) |
| * @version 1.1 |
| * |
| * @see InterpolatedTransform |
| * @see LinearTransform |
| * @see DatumShiftGrid |
| * |
| * @since 0.8 |
| * @module |
| */ |
| public class LocalizationGridBuilder extends TransformBuilder { |
| /** |
| * Tolerance threshold for comparing pixel coordinates relative to integer values. |
| */ |
| private static final double EPS = Numerics.COMPARISON_THRESHOLD; |
| |
| /** |
| * The transform for the linear part. |
| * Always created with a grid size specified to the constructor. |
| */ |
| private final LinearTransformBuilder linear; |
| |
| /** |
| * A temporary array for two-dimensional source coordinates. |
| * Used for reducing object allocations. |
| */ |
| private final int[] tmp = new int[SOURCE_DIMENSION]; |
| |
| /** |
| * Conversions from source real-world coordinates to grid indices before interpolation. |
| * If there is no such conversion to apply, then this is the identity transform. |
| */ |
| private LinearTransform sourceToGrid; |
| |
| /** |
| * The desired precision of inverse transformations in unit of source coordinates, or 0 in unspecified. |
| * If no {@link #sourceToGrid} transform has been specified, then this is in unit of grid cells |
| * (i.e. a value of 1 is the size of one grid cell). |
| * |
| * @see #setDesiredPrecision(double) |
| */ |
| private double precision; |
| |
| /** |
| * Arbitrary default {@link #precision} value, in unit of grid cells. Used if no explicit inverse transform |
| * precision has been specified. The {@code sourceToGrid} transform shall not be applied on this value. |
| * This default precision may change in any future SIS version. |
| */ |
| private static final double DEFAULT_PRECISION = 1E-7; |
| |
| /** |
| * The transform created by {@link #create(MathTransformFactory)}. |
| */ |
| private MathTransform transform; |
| |
| /** |
| * If the coordinates in some dimensions are cyclic, their periods. Otherwise {@code null}. |
| * Values are in units of the target CRS. For longitude wraparounds, the period is typically 360°. |
| * Array length shall be {@code linear.getTargetDimensions()} and non-cyclic dimensions shall have |
| * a period of zero (not {@link Double#NaN}, because we will use this array as a displacement vector). |
| * |
| * @see #resolveWraparoundAxis(int, int, double) |
| * @see ResidualGrid#periodVector |
| */ |
| private double[] periods; |
| |
| /** |
| * Creates a new, initially empty, builder for a localization grid of the given size. |
| * |
| * @param width the number of columns in the grid of target positions. |
| * @param height the number of rows in the grid of target positions. |
| */ |
| public LocalizationGridBuilder(final int width, final int height) { |
| linear = new LinearTransformBuilder(width, height); |
| sourceToGrid = MathTransforms.identity(SOURCE_DIMENSION); |
| } |
| |
| /** |
| * Creates a new, initially empty, builder for a localization grid of a size inferred from the given points. |
| * This constructor uses the given vectors for computing a grid size and the following initial conversion: |
| * |
| * <blockquote>({@code sourceX}, {@code sourceY}) → ({@code gridX}, {@code gridY})</blockquote> |
| * |
| * Above conversion can be obtained by {@link #getSourceToGrid()}. |
| * |
| * <p>Values in the given vectors should be integers, but this constructor is tolerant to non-integer values |
| * if they have a constant offset (typically 0.5) relative to integer values. The two vectors do not need to |
| * have the same length (i.e. {@code sourceX[i]} are not necessarily related to {@code sourceY[i]}).</p> |
| * |
| * @param sourceX all possible <var>x</var> inputs before conversion to grid coordinates. |
| * @param sourceY all possible <var>y</var> inputs before conversion to grid coordinates. |
| * @throws ArithmeticException if this constructor can not infer a reasonable grid size from the given vectors. |
| */ |
| public LocalizationGridBuilder(final Vector sourceX, final Vector sourceY) { |
| final Matrix fromGrid = new Matrix3(); |
| final int width = infer(sourceX, fromGrid, 0); |
| final int height = infer(sourceY, fromGrid, 1); |
| linear = new LinearTransformBuilder(width, height); |
| try { |
| sourceToGrid = MathTransforms.linear(fromGrid).inverse(); |
| } catch (NoninvertibleTransformException e) { |
| // Should not happen because infer(…) verified that the coefficients are okay. |
| throw (ArithmeticException) new ArithmeticException(e.getLocalizedMessage()).initCause(e); |
| } |
| } |
| |
| /** |
| * Creates a new builder for a localization grid inferred from the given provider of control points. |
| * The {@linkplain LinearTransformBuilder#getSourceDimensions() number of source dimensions} in the |
| * given {@code localizations} argument shall be 2. The {@code localization} can be used in two ways: |
| * |
| * <ul class="verbose"> |
| * <li>If the {@code localizations} instance has been |
| * {@linkplain LinearTransformBuilder#LinearTransformBuilder(int...) created with a fixed grid size}, |
| * then that instance is used as-is — it is not copied. It is okay to specify an empty instance and |
| * to provide control points later by calls to {@link #setControlPoint(int, int, double...)}.</li> |
| * <li>If the {@code localizations} instance has been |
| * {@linkplain LinearTransformBuilder#LinearTransformBuilder() created for a grid of unknown size}, |
| * then this constructor tries to infer a grid size by inspection of the control points present in |
| * {@code localizations} at the time this constructor is invoked. Changes in {@code localizations} |
| * after construction will not be reflected in this new builder.</li> |
| * </ul> |
| * |
| * @param localizations the provider of control points for which to create a localization grid. |
| * @throws ArithmeticException if this constructor can not infer a reasonable grid size from the given localizations. |
| * |
| * @since 1.0 |
| */ |
| public LocalizationGridBuilder(final LinearTransformBuilder localizations) { |
| ArgumentChecks.ensureNonNull("localizations", localizations); |
| int n = localizations.getGridDimensions(); |
| if (n == SOURCE_DIMENSION) { |
| linear = localizations; |
| sourceToGrid = MathTransforms.identity(SOURCE_DIMENSION); |
| } else { |
| if (n < 0) { |
| final Vector[] sources = localizations.sources(); |
| n = sources.length; |
| if (n == SOURCE_DIMENSION) { |
| final Matrix fromGrid = new Matrix3(); |
| final int width = infer(sources[0], fromGrid, 0); |
| final int height = infer(sources[1], fromGrid, 1); |
| linear = new LinearTransformBuilder(width, height); |
| linear.setControlPoints(localizations.getControlPoints()); |
| try { |
| sourceToGrid = MathTransforms.linear(fromGrid).inverse(); |
| } catch (NoninvertibleTransformException e) { |
| throw (ArithmeticException) new ArithmeticException(e.getLocalizedMessage()).initCause(e); |
| } |
| linear.setLinearizers(localizations); |
| return; |
| } |
| } |
| throw new IllegalArgumentException(Resources.format( |
| Resources.Keys.MismatchedTransformDimension_3, 0, SOURCE_DIMENSION, n)); |
| } |
| } |
| |
| /** |
| * Infers a grid size by searching for the greatest common divisor (GCD) for values in the given vector. |
| * The vector values should be integers, but this method is tolerant to constant offsets (typically 0.5). |
| * The GCD is taken as a "grid to source" scale factor and the minimal value as the translation term. |
| * Those two values are stored in the {@code dim} row of the given matrix. |
| * |
| * @param source the vector of values for which to get the GCD and minimum value. |
| * @param fromGrid matrix where to store the minimum value and the GCD. |
| * @param dim index of the matrix row to update. |
| * @return grid size. |
| */ |
| private static int infer(final Vector source, final Matrix fromGrid, final int dim) { |
| final NumberRange<?> range = source.range(); |
| final double min = range.getMinDouble(true); |
| final double span = range.getMaxDouble(true) - min; |
| final Number increment = source.increment(EPS * span); |
| double inc; |
| if (increment != null) { |
| inc = increment.doubleValue(); |
| } else { |
| inc = span; |
| final int size = source.size(); |
| for (int i=0; i<size; i++) { |
| double v = source.doubleValue(i) - min; |
| if (Math.abs(v % inc) > EPS) { |
| do { |
| final double r = (inc % v); // Both 'inc' and 'v' are positive, so 'r' will be positive too. |
| inc = v; |
| v = r; |
| } while (Math.abs(v) > EPS); |
| } |
| } |
| } |
| /* |
| * Compute the size from the increment that we found. If the size is larger than the vector length, |
| * consider as too large. The rational is that attempt to create a localization grid with that size |
| * would fail anyway, because it would contain holes where no value is defined. A limit is important |
| * for preventing useless allocation of large arrays - https://issues.apache.org/jira/browse/SIS-407 |
| */ |
| fromGrid.setElement(dim, dim, inc); |
| fromGrid.setElement(dim, SOURCE_DIMENSION, min); |
| final double n = span / inc; |
| if (n >= 0.5 && n < source.size() - 0.5) { // Compare as 'double' in case the value is large. |
| return ((int) Math.round(n)) + 1; |
| } |
| throw new ArithmeticException(Resources.format(Resources.Keys.CanNotInferGridSizeFromValues_1, range)); |
| } |
| |
| /** |
| * Throws {@link IllegalStateException} if this builder can not be modified anymore. |
| */ |
| private void ensureModifiable() throws IllegalStateException { |
| if (!linear.isModifiable()) { |
| throw new IllegalStateException(Errors.format(Errors.Keys.UnmodifiableObject_1, LocalizationGridBuilder.class)); |
| } |
| } |
| |
| /** |
| * Sets the desired precision of <em>inverse</em> transformations, in units of source coordinates. |
| * If a conversion from "real world" to grid coordinates {@linkplain #setSourceToGrid has been specified}, |
| * then the given precision is in "real world" units. Otherwise the precision is in units of grid cells |
| * (i.e. a value of 1 is the size of one grid cell). |
| * |
| * <div class="note"><b>Note:</b> |
| * there is no method for setting the desired target precision because forward transformations <em>precision</em> |
| * (not to be confused with <em>accuracy</em>) are limited only by rounding errors. Of course the accuracy of both |
| * forward and inverse transformations still limited by the accuracy of given control points and the grid resolution. |
| * </div> |
| * |
| * @param precision desired precision of the results of inverse transformations. |
| * @throws IllegalStateException if {@link #create(MathTransformFactory) create(…)} has already been invoked. |
| * |
| * @see DatumShiftGrid#getCellPrecision() |
| */ |
| public void setDesiredPrecision(final double precision) { |
| ensureModifiable(); |
| ArgumentChecks.ensureStrictlyPositive("precision", precision); |
| this.precision = precision; |
| } |
| |
| /** |
| * Returns the desired precision of <em>inverse</em> transformations, in units of source coordinates. |
| * This is the precision sets by the last call to {@link #setDesiredPrecision(double)}. |
| * |
| * @return desired precision of the results of inverse transformations. |
| */ |
| public double getDesiredPrecision() { |
| return precision; |
| } |
| |
| /** |
| * Defines relationship between "real-world" source coordinates and grid coordinates. |
| * The given transform is usually two-dimensional, in which case conversions from (<var>x</var>,<var>y</var>) |
| * source coordinates to ({@code gridX}, {@code gridY}) indices can be done with the following formulas: |
| * <ul> |
| * <li><var>gridX</var> = (<var>x</var> - <var>x₀</var>) / <var>Δx</var></li> |
| * <li><var>gridY</var> = (<var>y</var> - <var>y₀</var>) / <var>Δy</var></li> |
| * </ul> |
| * |
| * where: |
| * <ul> |
| * <li>(<var>x₀</var>, <var>y₀</var>) is the coordinate of the center of the cell at grid index (0,0).</li> |
| * <li><var>Δx</var> and <var>Δy</var> are the distances between two cells on the <var>x</var> and <var>y</var> |
| * axes respectively, in the same unit of measurement than the one documented in the |
| * {@link #setDesiredPrecision(double)} method.</li> |
| * </ul> |
| * |
| * The {@code coordinateToGrid} transform for the above formulas can be represented by the following matrix: |
| * |
| * {@preformat math |
| * ┌ ┐ |
| * │ 1/Δx 0 -x₀/Δx │ |
| * │ 0 1/Δy -y₀/Δy │ |
| * │ 0 0 1 │ |
| * └ ┘ |
| * } |
| * |
| * If this method is never invoked, then the default conversion is identity. |
| * If a {@linkplain #setDesiredPrecision(double) desired precision} has been specified before this method call, |
| * it is caller's responsibility to convert that value to new source units if needed. |
| * |
| * @param sourceToGrid conversion from the "real world" source coordinates to grid indices including fractional parts. |
| * @throws IllegalStateException if {@link #create(MathTransformFactory) create(…)} has already been invoked. |
| * |
| * @see DatumShiftGrid#getCoordinateToGrid() |
| */ |
| public void setSourceToGrid(final LinearTransform sourceToGrid) { |
| ensureModifiable(); |
| ArgumentChecks.ensureNonNull("sourceToGrid", sourceToGrid); |
| int isTarget = 0; |
| int dim = sourceToGrid.getSourceDimensions(); |
| if (dim >= SOURCE_DIMENSION) { |
| isTarget = 1; |
| dim = sourceToGrid.getTargetDimensions(); |
| if (dim == SOURCE_DIMENSION) { |
| this.sourceToGrid = sourceToGrid; |
| return; |
| } |
| } |
| throw new MismatchedDimensionException(Resources.format( |
| Resources.Keys.MismatchedTransformDimension_3, isTarget, SOURCE_DIMENSION, dim)); |
| } |
| |
| /** |
| * Returns the current relationship between "real-world" source coordinates and grid coordinates. |
| * This is the value set by the last call to {@link #setSourceToGrid(LinearTransform)}. |
| * If that setter method has never been invoked, then this is an automatically computed transform |
| * if the grid coordinates {@linkplain #LocalizationGridBuilder(Vector, Vector) have been specified |
| * to the constructor}, or the identity transform {@linkplain #LocalizationGridBuilder(int, int) otherwise}. |
| * |
| * @return the current relationship between "real-world" source coordinates and grid coordinates. |
| */ |
| public LinearTransform getSourceToGrid() { |
| return sourceToGrid; |
| } |
| |
| /** |
| * Returns the envelope of source coordinates. The {@code fullArea} argument control whether |
| * the returned envelope shall encompass full surface of every cells or only their centers: |
| * <ul> |
| * <li>If {@code true}, then the returned envelope encompasses full cell surfaces, |
| * from lower border to upper border. In other words, the returned envelope encompasses all |
| * {@linkplain org.opengis.referencing.datum.PixelInCell#CELL_CORNER cell corners}.</li> |
| * <li>If {@code false}, then the returned envelope encompasses only |
| * {@linkplain org.opengis.referencing.datum.PixelInCell#CELL_CENTER cell centers}, inclusive.</li> |
| * </ul> |
| * |
| * This is the envelope of the grid domain (i.e. the ranges of valid {@code gridX} and {@code gridY} argument |
| * values in calls to {@code get/setControlPoint(…)} methods) transformed as below: |
| * <ol> |
| * <li>expanded by ½ cell on each side if {@code fullArea} is {@code true}</li> |
| * <li>transformed by the inverse of {@linkplain #getSourceToGrid() source to grid} transform.</li> |
| * </ol> |
| * |
| * @param fullArea whether the the envelope shall encompass the full cell surfaces instead than only their centers. |
| * @return the envelope of grid points, from lower corner to upper corner. |
| * @throws IllegalStateException if the grid points are not yet known. |
| * @throws TransformException if the envelope can not be calculated. |
| * |
| * @see LinearTransformBuilder#getSourceEnvelope() |
| * |
| * @since 1.0 |
| */ |
| public Envelope getSourceEnvelope(final boolean fullArea) throws TransformException { |
| Envelope envelope = linear.getSourceEnvelope(); |
| if (fullArea) { |
| for (int i = envelope.getDimension(); --i >= 0;) { |
| final GeneralEnvelope ge = GeneralEnvelope.castOrCopy(envelope); |
| ge.setRange(i, ge.getLower(i) - 0.5, |
| ge.getUpper(i) + 0.5); |
| envelope = ge; |
| } |
| } |
| return Envelopes.transform(sourceToGrid.inverse(), envelope); |
| } |
| |
| /** |
| * Sets all control points. The length of given vectors must be equal to the total number of cells in the grid. |
| * The first vector provides the <var>x</var> coordinates; the second vector provides the <var>y</var> coordinates, |
| * <i>etc.</i>. Coordinates are stored in row-major order (column index varies faster, followed by row index). |
| * |
| * @param coordinates coordinates in each target dimensions, stored in row-major order. |
| * @throws IllegalStateException if {@link #create(MathTransformFactory) create(…)} has already been invoked. |
| * |
| * @since 1.0 |
| */ |
| public void setControlPoints(final Vector... coordinates) { |
| ensureModifiable(); |
| ArgumentChecks.ensureNonNull("coordinates", coordinates); |
| linear.setControlPoints(coordinates); |
| } |
| |
| /** |
| * Sets a single matching control point pair. Source position is assumed precise and target position is assumed uncertain. |
| * If the given source position was already associated with another target position, then the old target position is discarded. |
| * |
| * <p>If a {@linkplain #getSourceToGrid() source to grid} conversion exists, it shall have been applied |
| * by the caller for computing the ({@code gridX}, {@code gridY}) coordinates given to this method.</p> |
| * |
| * @param gridX the column index in the grid where to store the given target position. |
| * @param gridY the row index in the grid where to store the given target position. |
| * @param target the target coordinates, assumed uncertain. |
| * @throws IllegalStateException if {@link #create(MathTransformFactory) create(…)} has already been invoked. |
| * @throws IllegalArgumentException if the {@code x} or {@code y} coordinate value is out of grid range. |
| * @throws MismatchedDimensionException if the target position does not have the expected number of dimensions. |
| */ |
| public void setControlPoint(final int gridX, final int gridY, final double... target) { |
| ensureModifiable(); |
| tmp[0] = gridX; |
| tmp[1] = gridY; |
| linear.setControlPoint(tmp, target); |
| } |
| |
| /** |
| * Returns a single target coordinate for the given source coordinate, or {@code null} if none. |
| * If {@linkplain #addLinearizers linearizers} have been specified and {@link #create create(…)} |
| * has already been invoked, then the control points may be projected using one of the linearizers. |
| * |
| * @param gridX the column index in the grid where to read the target position. |
| * @param gridY the row index in the grid where to read the target position. |
| * @return the target coordinates associated to the given source, or {@code null} if none. |
| * @throws IllegalArgumentException if the {@code x} or {@code y} coordinate value is out of grid range. |
| */ |
| public double[] getControlPoint(final int gridX, final int gridY) { |
| tmp[0] = gridX; |
| tmp[1] = gridY; |
| return linear.getControlPoint(tmp); |
| } |
| |
| /** |
| * Returns a row of coordinate values in the given dimension. |
| * The returned vector is a view; changes in the returned vector will be reflected in this builder. |
| * |
| * @param dimension the target dimension for which to get coordinate values. |
| * @param row index of the row to get. |
| * @return coordinate values of the specified row in the specified dimension. |
| * |
| * @since 1.0 |
| */ |
| public Vector getRow(final int dimension, final int row) { |
| tmp[0] = 0; |
| tmp[1] = row; |
| return linear.getTransect(dimension, tmp, 0); |
| } |
| |
| /** |
| * Returns a column of coordinate values in the given dimension. |
| * The returned vector is a view; changes in the returned vector will be reflected in this builder. |
| * |
| * @param dimension the target dimension for which to get coordinate values. |
| * @param column index of the column to get. |
| * @return coordinate values of the specified column in the specified dimension. |
| * |
| * @since 1.0 |
| */ |
| public Vector getColumn(final int dimension, final int column) { |
| tmp[0] = column; |
| tmp[1] = 0; |
| return linear.getTransect(dimension, tmp, 1); |
| } |
| |
| /** |
| * Tries to remove discontinuities in coordinates values caused by anti-meridian crossing. |
| * This method can be invoked when the localization grid may cross the anti-meridian, |
| * where longitude values may suddenly jump from +180° to -180° or conversely. |
| * This method walks through the coordinate values of the given dimension (typically the longitudes dimension) |
| * in the given direction (grid rows or grid columns). |
| * If a difference greater than {@code period/2} (typically 180°) is found between two consecutive values, |
| * then a multiple of {@code period} (typically 360°) is added or subtracted in order to make a value as close |
| * as possible from its previous value. |
| * |
| * <p>This method needs a direction to be specified:</p> |
| * <ul> |
| * <li>Direction 0 means that each value is compared with the value in the previous column, |
| * except the value in the first column which is compared to the value in previous row.</li> |
| * <li>Direction 1 means that each value is compared with the value in the previous row, |
| * except the value in the first row which is compared to the value in previous column.</li> |
| * </ul> |
| * The recommended value is the direction of most stable values. Typically, longitude values increase with column indices |
| * and are almost constant when increasing row indices. In such case, the recommended direction is 1 for comparing each |
| * value with the value in previous row, since that value should be closer than the value in previous column. |
| * |
| * <div class="note"><b>Example:</b> |
| * for a grid of (<var>longitude</var>, <var>latitude</var>) values in decimal degrees where longitude values |
| * vary (increase or decrease) with increasing column indices and latitude values vary (increase or decrease) |
| * with increasing row indices, the the following method should be invoked for protecting the grid against |
| * discontinuities on anti-meridian: |
| * |
| * {@preformat java |
| * grid.resolveWraparoundAxis(0, 1, 360); |
| * } |
| * </div> |
| * |
| * @param dimension the dimension to process. |
| * This is 0 for longitude dimension in a (<var>longitudes</var>, <var>latitudes</var>) grid. |
| * @param direction the direction to walk through: 0 for columns or 1 for rows. |
| * The recommended direction is the direction of most stable values, typically 1 (rows) for longitudes. |
| * @param period that wraparound range (typically 360° for longitudes). Must be strictly positive. |
| * @return the range of coordinate values in the specified dimension after correction for wraparound values. |
| * @throws IllegalStateException if this method has already been invoked for the same dimension, |
| * or if {@link #create(MathTransformFactory) create(…)} has already been invoked. |
| * |
| * @since 1.0 |
| */ |
| public NumberRange<Double> resolveWraparoundAxis(final int dimension, final int direction, final double period) { |
| ensureModifiable(); |
| ArgumentChecks.ensureBetween("dimension", 0, linear.getTargetDimensions() - 1, dimension); |
| ArgumentChecks.ensureBetween("direction", 0, linear.getSourceDimensions() - 1, direction); |
| ArgumentChecks.ensureStrictlyPositive("period", period); |
| if (periods == null) { |
| periods = new double[linear.getTargetDimensions()]; |
| } |
| if (periods[dimension] != 0) { |
| throw new IllegalStateException(Errors.format( |
| Errors.Keys.ValueAlreadyDefined_1, Strings.bracket("periods", dimension))); |
| } |
| periods[dimension] = period; |
| return linear.resolveWraparoundAxis(dimension, direction, period); |
| } |
| |
| /** |
| * Adds transforms to potentially apply on target coordinates before to compute the transform. |
| * This method can be invoked if the departure from a linear transform is too large, resulting |
| * in {@link InterpolatedTransform} to fail with "no convergence error" messages. |
| * If linearizers have been specified, then the {@link #create(MathTransformFactory)} method |
| * will try to apply each transform on target coordinates and check which one results in the |
| * best correlation coefficients. It may be none. |
| * |
| * <p>The linearizers are specified as {@link MathTransform}s from current target coordinates |
| * to other spaces where <cite>sources to new targets</cite> transforms may be more linear. |
| * The keys in the map are arbitrary identifiers used in {@link #toString()} for debugging purpose. |
| * The {@code dimensions} argument specifies which target dimensions to project and can be null or omitted |
| * if the projections shall be applied on all target coordinates. It is possible to invoke this method many |
| * times with different {@code dimensions} argument values.</p> |
| * |
| * @param projections projections from current target coordinates to other spaces which may result in more linear transforms. |
| * @param dimensions the target dimensions to project, or null or omitted for projecting all target dimensions. |
| * If non-null and non-empty, then all transforms in the {@code projections} map shall have a |
| * number of source and target dimensions equals to the length of this array. |
| * @throws IllegalStateException if {@link #create(MathTransformFactory) create(…)} has already been invoked. |
| * |
| * @see LinearTransformBuilder#addLinearizers(Map, int...) |
| * |
| * @since 1.0 |
| */ |
| public void addLinearizers(final Map<String,MathTransform> projections, int... dimensions) { |
| ensureModifiable(); |
| linear.addLinearizers(projections, dimensions); |
| } |
| |
| /** |
| * Creates a transform from the source points to the target points. |
| * This method assumes that source points are precise and all uncertainty is in the target points. |
| * If this transform is close enough to an affine transform, then an instance of {@link LinearTransform} is returned. |
| * |
| * <p>If this method is invoked more than once, the previously created transform instance is returned.</p> |
| * |
| * @param factory the factory to use for creating the transform, or {@code null} for the default factory. |
| * The {@link MathTransformFactory#createAffineTransform(Matrix)} method of that factory |
| * shall return {@link LinearTransform} instances. |
| * @return the transform from source to target points. |
| * @throws FactoryException if the transform can not be created, |
| * for example because the target points have not be specified. |
| */ |
| @Override |
| public MathTransform create(final MathTransformFactory factory) throws FactoryException { |
| if (transform == null) { |
| final LinearTransform gridToCoord = linear.create(factory); |
| /* |
| * Make a first check about whether the result of above LinearTransformBuilder.create() call |
| * can be considered a good fit. If true, then we may return the linear transform directly. |
| */ |
| boolean isExact = true; |
| boolean isLinear = true; |
| for (final double c : linear.correlation()) { |
| isExact &= (c == 1); |
| if (!(c >= 0.9999)) { // Empirical threshold (may need to be revisited). |
| isLinear = false; |
| break; |
| } |
| } |
| MathTransform step; |
| if (isExact) { |
| step = MathTransforms.concatenate(sourceToGrid, gridToCoord); |
| } else { |
| final int width = linear.gridSize(0); |
| final int height = linear.gridSize(1); |
| final float[] residual = new float [SOURCE_DIMENSION * linear.gridLength]; |
| final double[] grid = new double[SOURCE_DIMENSION * width]; |
| double gridPrecision = precision; |
| try { |
| /* |
| * If the user specified a precision, we need to convert it from source units to grid units. |
| * We convert each dimension separately, then retain the largest magnitude of vector results. |
| */ |
| if (gridPrecision > 0 && !sourceToGrid.isIdentity()) { |
| final double[] vector = new double[sourceToGrid.getSourceDimensions()]; |
| final double[] offset = new double[sourceToGrid.getTargetDimensions()]; |
| double converted = 0; |
| for (int i=0; i<vector.length; i++) { |
| vector[i] = precision; |
| sourceToGrid.deltaTransform(vector, 0, offset, 0, 1); |
| final double length = MathFunctions.magnitude(offset); |
| if (length > converted) converted = length; |
| vector[i] = 0; |
| } |
| gridPrecision = converted; |
| } |
| /* |
| * Compute the residuals, i.e. the differences between the coordinates that we get by a linear |
| * transformation and the coordinates that we want to get. If at least one residual is greater |
| * than the desired precision, then the returned MathTransform will need to apply corrections |
| * after linear transforms. Those corrections will be done by InterpolatedTransform. |
| */ |
| final MathTransform coordToGrid = gridToCoord.inverse(); |
| for (int k=0,y=0; y<height; y++) { |
| tmp[0] = 0; |
| tmp[1] = y; |
| linear.getControlRow(tmp, grid); // Expected positions. |
| coordToGrid.transform(grid, 0, grid, 0, width); // As grid coordinate. |
| for (int i=0,x=0; x<width; x++) { |
| final double dx = grid[i++] - x; |
| final double dy = grid[i++] - y; |
| isLinear &= (dx <= gridPrecision); |
| isLinear &= (dy <= gridPrecision); |
| residual[k++] = (float) dx; |
| residual[k++] = (float) dy; |
| } |
| } |
| if (isLinear) { |
| step = MathTransforms.concatenate(sourceToGrid, gridToCoord); |
| } else { |
| step = InterpolatedTransform.createGeodeticTransformation(nonNull(factory), |
| new ResidualGrid(sourceToGrid, gridToCoord, width, height, residual, |
| (gridPrecision > 0) ? gridPrecision : DEFAULT_PRECISION, periods)); |
| } |
| } catch (TransformException e) { |
| throw new FactoryException(e); // Should never happen. |
| } |
| } |
| /* |
| * At this point we finished to compute the transformation to target coordinates. |
| * If those target coordinates have been modified in order to make that step more |
| * linear, apply the inverse transformation after the step. |
| */ |
| final Optional<MathTransform> linearizer = linear.linearizer(); |
| if (linearizer.isPresent()) try { |
| step = factory.createConcatenatedTransform(step, linearizer.get().inverse()); |
| } catch (NoninvertibleTransformException e) { |
| throw new InvalidGeodeticParameterException(Resources.format( |
| Resources.Keys.NonInvertibleOperation_1, linear.linearizerID()), e); |
| } |
| transform = step; // Set only after everything succeeded. |
| } |
| return transform; |
| } |
| |
| /** |
| * Returns statistics of differences between values calculated by the given transform and actual values. |
| * The given math transform is typically the transform computed by {@link #create(MathTransformFactory)}, |
| * but not necessarily. The returned statistics are: |
| * |
| * <ol class="verbose"> |
| * <li>One {@code Statistics} instance for each target dimension, containing statistics about the differences between |
| * coordinates computed by the given transform and expected coordinates. For each (<var>i</var>,<var>j</var>) indices |
| * in this grid, the indices are transformed by a call to {@code mt.transform(…)} and the result is compared with the |
| * coordinates given by <code>{@linkplain #getControlPoint(int, int) getControlPoint}(i,j)</code>. |
| * Those statistics are identified by labels like “P → x” and “P → y” where <var>P</var> stands for pixel coordinates.</li> |
| * <li>One {@code Statistics} instance for each source dimension, containing statistics about the differences between |
| * coordinates computed by the <em>inverse</em> of the transform and expected coordinates. |
| * For each (<var>x</var>,<var>y</var>) control point in this grid, the points are transformed by a call |
| * to {@code mt.inverse().transform(…)} and the result is compared with the pixel indices of that point. |
| * Those statistics are identified by labels like “i ← P′” and “j ← P′” where <var>P′</var> stands for |
| * the control point.</li> |
| * </ol> |
| * |
| * @param mt the transform to test. |
| * @return statistics of difference between computed values and expected values for each target dimension. |
| * @throws NoninvertibleTransformException if an error occurred while inverting a transform. |
| * |
| * @since 1.0 |
| */ |
| public Statistics[] error(final MathTransform mt) throws NoninvertibleTransformException { |
| final int tgtDim = mt.getTargetDimensions(); |
| final double[] point = new double[Math.max(tgtDim, SOURCE_DIMENSION)]; |
| final Statistics[] stats = new Statistics[tgtDim + SOURCE_DIMENSION]; |
| final StringBuilder buffer = new StringBuilder(); |
| for (int i=0; i<stats.length; i++) { |
| buffer.setLength(0); |
| if (i < tgtDim) { |
| buffer.append("P → "); |
| if (i < 3) { |
| buffer.append((char) ('x' + i)); |
| } else { |
| buffer.append('z').append(i - 1); // After (x,y,z) continue with z2, z3, z4, etc. |
| } |
| } else { |
| buffer.append((char) ('i' + (i - tgtDim))).append(" ← P′"); |
| } |
| stats[i] = new Statistics(buffer.toString()); |
| } |
| /* |
| * If a linearizer has been applied, all target coordinates in this builder have been projected using |
| * that transform. We will need to apply the inverse transform in order to get back the original values. |
| */ |
| final Optional<MathTransform> linearizer = linear.linearizer(); |
| final MathTransform complete = linearizer.isPresent() ? linearizer.get().inverse() : null; |
| final MathTransform inverse = mt.inverse(); |
| final int width = linear.gridSize(0); |
| final int height = linear.gridSize(1); |
| for (int y=0; y<height; y++) { |
| for (int x=0; x<width; x++) { |
| point[0] = tmp[0] = x; |
| point[1] = tmp[1] = y; |
| final double[] expected; |
| try { |
| mt.transform(point, 0, point, 0, 1); |
| expected = linear.getControlPoint(tmp); |
| if (complete != null) { |
| complete.transform(expected, 0, expected, 0, 1); |
| } |
| } catch (TransformException e) { |
| continue; // Ignore the points that we fail to transform. |
| } |
| for (int i=0; i<tgtDim; i++) { |
| stats[i].accept(point[i] - expected[i]); |
| } |
| /* |
| * Transform the geographic point back to grid indices and check error. |
| */ |
| try { |
| inverse.transform(expected, 0, expected, 0, 1); |
| } catch (TransformException e) { |
| continue; // Ignore the points that we fail to transform. |
| } |
| for (int i=0; i<SOURCE_DIMENSION; i++) { |
| stats[tgtDim + i].accept(expected[i] - tmp[i]); |
| } |
| } |
| } |
| return stats; |
| } |
| |
| /** |
| * Returns a string representation of this builder in the given locale. |
| * Current implementation shows the following information: |
| * |
| * <ul> |
| * <li>Number of points.</li> |
| * <li>Linearizers and their correlation coefficients (if available).</li> |
| * <li>The linear component of the transform.</li> |
| * <li>Error statistics, as documented in the {@link #error(MathTransform)} method.</li> |
| * </ul> |
| * |
| * The string representation may change in any future version. |
| * |
| * @param locale the locale for formatting messages and some numbers, or {@code null} for the default. |
| * @return a string representation of this builder. |
| * |
| * @since 1.0 |
| */ |
| public String toString(final Locale locale) { |
| final StringBuilder buffer = new StringBuilder(400); |
| String lineSeparator = null; |
| try { |
| lineSeparator = linear.appendTo(buffer, getClass(), locale, Vocabulary.Keys.LinearTransformation); |
| if (transform != null) { |
| buffer.append(Strings.CONTINUATION_ITEM); |
| final Vocabulary vocabulary = Vocabulary.getResources(locale); |
| vocabulary.appendLabel(Vocabulary.Keys.Errors, buffer); |
| buffer.append(lineSeparator); |
| final StatisticsFormat sf; |
| if (locale != null) { |
| sf = StatisticsFormat.getInstance(locale); |
| } else { |
| sf = StatisticsFormat.getInstance(); |
| } |
| sf.format(error(transform), buffer); |
| } |
| } catch (IOException e) { |
| throw new UncheckedIOException(e); |
| } catch (NoninvertibleTransformException e) { |
| // Ignore - we will not report error statistics. |
| } |
| Strings.insertLineInLeftMargin(buffer, lineSeparator); |
| return buffer.toString(); |
| } |
| |
| /** |
| * Returns a string representation of this builder for debugging purpose. |
| * The string representation is for debugging purpose and may change in any future version. |
| * The default implementation delegates to {@link #toString(Locale)} with a null locale. |
| * |
| * @return a string representation of this builder. |
| * |
| * @since 1.0 |
| */ |
| @Override |
| public String toString() { |
| return toString(null); |
| } |
| } |
