/*
 * 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.0
 *
 * @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;

    /**
     * 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 {@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);
        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;
                        }
                    }
                } catch (TransformException e) {
                    throw new FactoryException(e);                                          // Should never happen.
                }
                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));
                }
            }
            /*
             * 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;
        }
        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);
    }
}