/*
 * 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.Map;
import java.util.Arrays;
import java.util.TreeMap;
import java.util.Objects;
import java.awt.image.RenderedImage;
import org.opengis.util.FactoryException;
import org.opengis.geometry.DirectPosition;
import org.opengis.metadata.extent.GeographicBoundingBox;
import org.opengis.referencing.datum.PixelInCell;
import org.opengis.referencing.operation.Matrix;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.operation.NoninvertibleTransformException;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.apache.sis.coverage.SampleDimension;
import org.apache.sis.feature.internal.Resources;
import org.apache.sis.util.ArraysExt;
import org.apache.sis.util.ArgumentChecks;
import org.apache.sis.util.internal.CollectionsExt;
import org.apache.sis.coverage.grid.j2d.ImageUtilities;
import org.apache.sis.referencing.CRS;
import org.apache.sis.referencing.util.DirectPositionView;
import org.apache.sis.referencing.util.WraparoundAxesFinder;
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.TransformSeparator;
import org.apache.sis.util.logging.Logging;

// Specific to the main branch:
import org.apache.sis.coverage.CannotEvaluateException;
import org.apache.sis.coverage.PointOutsideCoverageException;


/**
 * Default implementation of {@link GridCoverage.Evaluator} for interpolating values at given positions.
 * Values are computed by calls to {@link #apply(DirectPosition)} and are returned as {@code double[]}.
 *
 * <h2>Multi-threading</h2>
 * Evaluators are not thread-safe. An instance of {@code DefaultEvaluator} should be created
 * for each thread that need to compute sample values.
 *
 * <h2>Limitations</h2>
 * Current implementation performs nearest-neighbor sampling only. A future version will provide interpolations.
 *
 * @author  Johann Sorel (Geomatys)
 * @author  Martin Desruisseaux (Geomatys)
 *
 * @see GridCoverage#evaluator()
 */
class DefaultEvaluator implements GridCoverage.Evaluator {
    /**
     * The coverage in which to evaluate sample values.
     */
    private final GridCoverage coverage;

    /**
     * The coordinate reference system of input points given to this converter,
     * or {@code null} if assumed the same than the coverage CRS.
     * This is used by {@link #toGridPosition(DirectPosition)} for checking if {@link #inputToGrid} needs
     * to be recomputed. As long at the evaluated points have the same CRS, the same transform is reused.
     */
    private CoordinateReferenceSystem inputCRS;

    /**
     * The transform from {@link #inputCRS} to grid coordinates.
     * This is cached for avoiding the costly process of fetching a coordinate operation
     * in the common case where the coordinate reference systems did not changed.
     */
    private MathTransform inputToGrid;

    /**
     * Grid coordinates after {@link #inputToGrid} conversion.
     *
     * @see #toGridPosition(DirectPosition)
     */
    private FractionalGridCoordinates.Position position;

    /**
     * Array where to store sample values computed by {@link #apply(DirectPosition)}.
     * For performance reasons, the same array may be recycled on every method call.
     */
    double[] values;

    /**
     * Whether to return {@code null} instead of throwing an exception if given point
     * is outside coverage bounds.
     *
     * @see #isNullIfOutside()
     */
    private boolean nullIfOutside;


    // ―――――――― Following fields are for the handling of wraparound axes. ―――――――――――――――――――――

    /**
     * A bitmask of grid dimensions that need to be verified for wraparound axes.
     */
    private long wraparoundAxes;

    /**
     * Coverage extent converted to floating point numbers, only for the grid dimensions having
     * a bit set to 1 in {@link #wraparoundAxes} bitmask. The length of this array is the number
     * of bits set in {@link #wraparoundAxes} multiplied by 2. Elements are (lower, upper) tuples.
     */
    private double[] wraparoundExtent;

    /**
     * Transform from grid coordinates to the CRS where wraparound axes may exist.
     * It is sometimes the same transform than {@code gridToCRS} but not always.
     * It may differ for example if a projected CRS has been replaced by a geographic CRS.
     */
    private MathTransform gridToWraparound;

    /**
     * The span (maximum - minimum) of wraparound axes, with 0 value for axes that are not wraparound.
     * The length of this array may be shorter than the CRS number of dimensions if all remaining axes
     * are not wraparound axes.
     */
    private double[] periods;

    /**
     * The slice where to perform evaluation, or {@code null} if not yet computed.
     * This information allows to specify for example two-dimensional points for
     * evaluating in a three-dimensional data cube. This is used for completing
     * the missing coordinate values.
     *
     * @see #getDefaultSlice()
     * @see #setDefaultSlice(Map)
     */
    private Map<Integer,Long> slice;

    /**
     * Creates a new evaluator for the given coverage. This constructor is protected for allowing
     * {@link GridCoverage} subclasses to provide their own {@code DefaultEvaluator} implementations.
     * For using an evaluator, invoke {@link GridCoverage#evaluator()} instead.
     *
     * @param  coverage  the coverage for which to create an evaluator.
     *
     * @see GridCoverage#evaluator()
     */
    protected DefaultEvaluator(final GridCoverage coverage) {
        ArgumentChecks.ensureNonNull("coverage", coverage);
        this.coverage = coverage;
    }

    /**
     * Returns the coverage from which this evaluator is fetching sample values.
     * This is the coverage on which the {@link GridCoverage#evaluator()} method has been invoked.
     *
     * @return the source of sample values for this evaluator.
     */
    @Override
    public final GridCoverage getCoverage() {
        return coverage;
    }

    /**
     * Returns the default slice where to perform evaluation, or an empty map if unspecified.
     * Keys are dimensions from 0 inclusive to {@link GridGeometry#getDimension()} exclusive,
     * and values are the grid coordinates of the slice in the dimension specified by the key.
     *
     * <p>This information allows to invoke {@link #apply(DirectPosition)} with for example two-dimensional points
     * even if the underlying coverage is three-dimensional. The missing coordinate values are replaced by the
     * values provided in the map.</p>
     *
     * @return the default slice where to perform evaluation, or an empty map if unspecified.
     */
    @Override
    @SuppressWarnings("ReturnOfCollectionOrArrayField")     // Because the map is unmodifiable.
    public final Map<Integer,Long> getDefaultSlice() {
        if (slice == null) {
            final GridExtent extent = coverage.getGridGeometry().getExtent();
            slice = CollectionsExt.unmodifiableOrCopy(extent.getSliceCoordinates());
        }
        return slice;
    }

    /**
     * Sets the default slice where to perform evaluation when the points do not have enough dimensions.
     * A {@code null} argument restores the default value, which is to infer the slice from the coverage
     * grid geometry.
     *
     * @param  slice  the default slice where to perform evaluation, or an empty map if none.
     * @throws IllegalArgumentException if the map contains an illegal dimension or grid coordinate value.
     *
     * @see GridExtent#getSliceCoordinates()
     */
    @Override
    @SuppressWarnings("AssignmentToCollectionOrArrayFieldFromParameter")
    public void setDefaultSlice(Map<Integer,Long> slice) {
        if (!Objects.equals(this.slice, slice)) {
            if (slice != null) {
                slice = CollectionsExt.unmodifiableOrCopy(new TreeMap<>(slice));
                final GridExtent extent = coverage.getGridGeometry().getExtent();
                final int max = extent.getDimension() - 1;
                for (final Map.Entry<Integer,Long> entry : slice.entrySet()) {
                    final int dim = entry.getKey();
                    ArgumentChecks.ensureBetween("slice.key", 0, max, dim);
                    ArgumentChecks.ensureBetween(extent.getAxisIdentification(dim, dim).toString(),
                                        extent.getLow(dim), extent.getHigh(dim), entry.getValue());
                }
            }
            this.slice  = slice;
            inputCRS    = null;
            inputToGrid = null;
        }
    }

    /**
     * Returns {@code true} if this evaluator is allowed to wraparound coordinates that are outside the grid.
     * The initial value is {@code false}. This method may continue to return {@code false} even after a call
     * to {@code setWraparoundEnabled(true)} if no wraparound axis has been found in the coverage CRS.
     *
     * @return {@code true} if this evaluator may wraparound coordinates that are outside the grid.
     */
    @Override
    public boolean isWraparoundEnabled() {
        return (wraparoundAxes != 0);
    }

    /**
     * Specifies whether this evaluator is allowed to wraparound coordinates that are outside the grid.
     * If {@code true} and if a given coordinate is outside the grid, then this evaluator may translate
     * the point along a wraparound axis in an attempt to get the point inside the grid. For example, if
     * the coverage CRS has a longitude axis, then the evaluator may translate the longitude value by a
     * multiple of 360°.
     *
     * @param  allow  whether to allow wraparound of coordinates that are outside the grid.
     */
    @Override
    public void setWraparoundEnabled(final boolean allow) {
        wraparoundAxes = 0;
        if (allow) try {
            final WraparoundAxesFinder f = new WraparoundAxesFinder(coverage.getCoordinateReferenceSystem());
            if ((periods = f.periods()) != null) {
                final GridGeometry gridGeometry = coverage.getGridGeometry();
                final GridExtent extent = gridGeometry.getExtent();
                MathTransform gridToCRS = gridGeometry.getGridToCRS(PixelInCell.CELL_CENTER);
                gridToWraparound = MathTransforms.concatenate(gridToCRS, f.preferredToSpecified.inverse());
                final Matrix m = gridToWraparound.derivative(new DirectPositionView.Double(extent.getPointOfInterest(PixelInCell.CELL_CENTER)));
                /*
                 * `gridToWraparound` is the transform from grid coordinates to a CRS where wraparound axes exist.
                 * It may be the coverage CRS or its base CRS. The wraparound axes are identified by `periods`.
                 * The Jacobian matrix tells us which grid axes are dependencies of the wraparound axes.
                 *
                 * Note: the use of this matrix is not fully reliable with non-linear `gridToCRS` transforms
                 * because it is theoretically possible that a coefficient is zero at the point of interest
                 * by pure coincidence but would be non-zero at another location,
                 * But we think that it would require a transform with unlikely properties
                 * (e.g. transforming parallels to vertical straight lines at some places).
                 */
                for (int j = periods.length; --j >= 0;) {
                    if (periods[j] > 0) {                       // Find target dimensions (CRS) having wraparound axis.
                        for (int i = Math.min(m.getNumCol(), Long.SIZE); --i >= 0;) {
                            if (m.getElement(j, i) != 0) {
                                wraparoundAxes |= (1L << i);    // Mark sources (grid dimensions) dependent of target (CRS dimensions).
                            }
                        }
                    }
                }
                /*
                 * Get the grid extent only for the grid axes that are connected to wraparound CRS axes.
                 * There is at least one such axis, otherwise `periods` would have been null.
                 */
                wraparoundExtent = new double[Long.bitCount(wraparoundAxes) << 1];
                long axes = wraparoundAxes;
                int j = 0;
                do {
                    final int i = Long.numberOfTrailingZeros(axes);
                    wraparoundExtent[j++] = extent.getLow(i);
                    wraparoundExtent[j++] = extent.getHigh(i);
                    axes &= ~(1L << i);
                } while (axes != 0);
                assert wraparoundExtent.length == j : j;
            }
        } catch (TransformException e) {
            recoverableException("setWraparoundEnabled", e);
        }
    }

    /**
     * Returns whether to return {@code null} instead of throwing an exception if a point is outside coverage bounds.
     * The default value is {@code false}, which means that the default {@link #apply(DirectPosition)} behavior is to
     * throw {@link PointOutsideCoverageException} for points outside bounds.
     *
     * @return whether {@link #apply(DirectPosition)} return {@code null} for points outside coverage bounds.
     */
    @Override
    public boolean isNullIfOutside() {
        return nullIfOutside;
    }

    /**
     * Sets whether to return {@code null} instead of throwing an exception if a point is outside coverage bounds.
     * The default value is {@code false}. Setting this flag to {@code true} may improve performances if the caller
     * expects that many points will be outside coverage bounds, since it reduces the amount of exceptions to be
     * created.
     *
     * @param  flag  whether {@link #apply(DirectPosition)} should use {@code null} return value instead of
     *               {@link PointOutsideCoverageException} for signaling that a point is outside coverage bounds.
     */
    @Override
    public void setNullIfOutside(final boolean flag) {
        nullIfOutside = flag;
    }

    /**
     * Returns a sequence of double values for a given point in the coverage.
     * The CRS of the given point may be any coordinate reference system;
     * coordinate conversions will be applied as needed.
     * If the CRS of the point is undefined, then it is assumed to be the {@linkplain #getCoverage() coverage} CRS.
     * The returned sequence includes a value for each {@linkplain SampleDimension sample dimension}.
     *
     * <p>The default interpolation type used when accessing grid values for points which fall between
     * grid cells is nearest neighbor. This default interpolation method may change in future version.</p>
     *
     * <p>The default implementation invokes {@link GridCoverage#render(GridExtent)} for a small region
     * around the point. Subclasses should override with more efficient implementation.</p>
     *
     * @param  point   the position where to evaluate.
     * @return the sample values at the specified point, or {@code null} if the point is outside the coverage.
     *         For performance reason, this method may return the same array
     *         on every method call by overwriting previous values.
     *         Callers should not assume that the array content stay valid for a long time.
     * @throws PointOutsideCoverageException if the evaluation failed because the input point
     *         has invalid coordinates and the {@link #isNullIfOutside()} flag is {@code false}.
     * @throws CannotEvaluateException if the values cannot be computed at the specified coordinates
     *         for another reason. This exception may be thrown if the coverage data type cannot be
     *         converted to {@code double} by an identity or widening conversion.
     *         Subclasses may relax this constraint if appropriate.
     */
    @Override
    public double[] apply(final DirectPosition point) throws CannotEvaluateException {
        /*
         * TODO: instead of restricting to a single point, keep the automatic size (1 or 2),
         * invoke render for each slice, then interpolate. We would keep a value of 1 in the
         * size array if we want to disable interpolation in some particular axis (e.g. time).
         *
         * See https://issues.apache.org/jira/browse/SIS-576
         */
        final GridGeometry gridGeometry = coverage.gridGeometry;
        final long[] size = new long[gridGeometry.getDimension()];
        Arrays.fill(size, 1);
        try {
            final FractionalGridCoordinates gc = toGridPosition(point);
            try {
                final GridExtent subExtent = gc.toExtent(gridGeometry.extent, size, nullIfOutside);
                if (subExtent != null) {
                    return evaluate(coverage.render(subExtent), 0, 0);
                }
            } catch (ArithmeticException | IndexOutOfBoundsException | DisjointExtentException ex) {
                if (!nullIfOutside) {
                    throw (PointOutsideCoverageException) new PointOutsideCoverageException(
                            gc.pointOutsideCoverage(gridGeometry.extent)).initCause(ex);
                }
            }
        } catch (PointOutsideCoverageException ex) {
            throw ex;
        } catch (RuntimeException | FactoryException | TransformException ex) {
            throw new CannotEvaluateException(ex.getMessage(), ex);
        }
        return null;        // May reach this point only if `nullIfOutside` is true.
    }

    /**
     * Gets sample values from the given image at the given index. This method does not verify
     * explicitly if the coordinates are out of bounds; we rely on the checks performed by the
     * image and sample model implementations.
     *
     * @param  data  the data from which to get the sample values.
     * @param  x     column index of the value to get.
     * @param  y     row index of the value to get.
     * @return the sample values. The same array may be recycled on every method call.
     * @throws ArithmeticException if an integer overflow occurred while computing indices.
     * @throws IndexOutOfBoundsException if a coordinate is out of bounds.
     */
    final double[] evaluate(final RenderedImage data, final int x, final int y) {
        final int tx = ImageUtilities.pixelToTileX(data, x);
        final int ty = ImageUtilities.pixelToTileY(data, y);
        return values = data.getTile(tx, ty).getPixel(x, y, values);
    }

    /**
     * Converts the specified geospatial position to grid coordinates.
     * If the given position is associated to a non-null coordinate reference system (CRS) different than the
     * {@linkplain #getCoverage() coverage} CRS, then this method automatically transforms that position to the
     * {@linkplain GridCoverage#getCoordinateReferenceSystem() coverage CRS} before to compute grid coordinates.
     *
     * <p>This method does not put any restriction on the grid coordinates result.
     * The result may be outside the {@linkplain GridGeometry#getExtent() grid extent}
     * if the {@linkplain GridGeometry#getGridToCRS(PixelInCell) grid to CRS} transform allows it.</p>
     *
     * @param  point  geospatial coordinates (in arbitrary CRS) to transform to grid coordinates.
     * @return the grid coordinates for the given geospatial coordinates.
     * @throws IncompleteGridGeometryException if the {@linkplain GridCoverage#getGridGeometry() grid geometry}
     *         does not define a "grid to CRS" transform, or if the given point has a non-null CRS but the
     *         coverage does not {@linkplain GridCoverage#getCoordinateReferenceSystem() have a CRS}.
     * @throws TransformException if the given coordinates cannot be transformed.
     *
     * @see FractionalGridCoordinates#toPosition(MathTransform)
     */
    @Override
    public FractionalGridCoordinates toGridCoordinates(final DirectPosition point) throws TransformException {
        ArgumentChecks.ensureNonNull("point", point);
        try {
            return new FractionalGridCoordinates(toGridPosition(point));
        } catch (FactoryException e) {
            throw new TransformException(e.getMessage(), e);
        }
    }

    /**
     * Updates the grid {@linkplain #position} with the given geospatial position.
     * This is the implementation of {@link #toGridCoordinates(DirectPosition)} except
     * that it avoid creating a new {@link FractionalGridCoordinates} on each method call.
     *
     * @param  point  the geospatial position.
     * @return the given position converted to grid coordinates (possibly out of grid bounds).
     * @throws FactoryException if no operation is found form given point CRS to coverage CRS.
     * @throws TransformException if the given position cannot be converted.
     */
    final FractionalGridCoordinates.Position toGridPosition(final DirectPosition point)
            throws FactoryException, TransformException
    {
        /*
         * If the `inputToGrid` transform has not yet been computed or is outdated, compute now.
         * The result will be cached and reused as long as the `inputCRS` is the same.
         */
        final CoordinateReferenceSystem crs = point.getCoordinateReferenceSystem();
        if (crs != inputCRS || inputToGrid == null) {
            setInputCRS(crs);
        }
        /*
         * Transform geospatial coordinates to grid coordinates. Result is unconditionally stored
         * in the `position` object, which will be copied by the caller if needed for public API.
         */
        final DirectPosition result = inputToGrid.transform(point, position);
        if (result != position) {
            // Should not happen, but be paranoiac.
            final double[] coordinates = position.coordinates;
            System.arraycopy(result.getCoordinate(), 0, coordinates, 0, coordinates.length);
        }
        /*
         * In most cases, the work of this method ends here. The remaining code in this method
         * is for handling wraparound axes. If a coordinate is outside the coverage extent,
         * check if a wraparound on some axes would bring the coordinates inside the extent.
         */
        long axes = wraparoundAxes;
        if (axes != 0) {
            double[] coordinates = position.coordinates;
            long outsideAxes = 0;
            int j = 0;
            do {
                final int i = Long.numberOfTrailingZeros(axes);
                final double c = coordinates[i];
                double border;
                if (c < (border = wraparoundExtent[j++]) || c > (border = wraparoundExtent[j])) {
                    /*
                     * Detected that the point is outside the grid extent along an axis where wraparound is possible.
                     * The first time that we find such axis, expand the coordinates array for storing two points.
                     * The two points will have the same coordinates, except on all axes where the point is outside.
                     * On those axes, the coordinate of the first point is set to the closest border of the grid.
                     */
                    if (outsideAxes == 0) {
                        final int n = coordinates.length;
                        coordinates = Arrays.copyOf(coordinates, 2*Math.max(n, gridToWraparound.getTargetDimensions()));
                        System.arraycopy(coordinates, 0, coordinates, n, n);
                    }
                    coordinates[i] = border;
                    outsideAxes |= (1L << i);
                }
                j++;    // Outside above `if (…)` statement because increment needs to be unconditional.
                axes &= ~(1L << i);
            } while (axes != 0);
            assert wraparoundExtent.length == j : j;
            /*
             * If a coordinate was found outside the grid, transform to a CRS where we can apply shift.
             * It may be the same CRS than the coverage CRS or the source CRS, but not necessarily.
             * For example if the CRS is projected, then we need to use a geographic intermediate CRS.
             * In the common case where the source CRS is already geographic, the second point in the
             * `coordinates` array after `transform(…)` will contain the same coordinates as `point`,
             * but potentially with more dimensions.
             */
            if (outsideAxes != 0) {
                gridToWraparound.transform(coordinates, 0, coordinates, 0, 2);
                final int s = gridToWraparound.getTargetDimensions();
                for (int i = periods.length; --i >= 0;) {
                    final double period = periods[i];
                    if (period > 0) {
                        /*
                         * Compute the shift that was necessary for moving the point inside the grid,
                         * then round that shift to an integer amount of periods. Modify the original
                         * coordinate by applying that modified translation.
                         */
                        final int oi = i + s;                       // Index of original coordinates.
                        double shift = coordinates[i] - coordinates[oi];
                        shift = Math.copySign(Math.ceil(Math.abs(shift) / period), shift) * period;
                        coordinates[oi] += shift;
                    }
                }
                /*
                 * Convert back the shifted point to grid coordinates, then check again if the new point
                 * is inside the grid extent. If this is not the case, we will return the old position
                 * on the assumption that it will be less confusing to the user.
                 */
                gridToWraparound.inverse().transform(coordinates, s, coordinates, 0, 1);
                j = 0;
                do {
                    final int i = Long.numberOfTrailingZeros(outsideAxes);
                    final double c = coordinates[i];
                    if (c < wraparoundExtent[j++] || c > wraparoundExtent[j++]) {
                        return position;
                    }
                    outsideAxes &= ~(1L << i);
                } while (outsideAxes != 0);
                /*
                 * Copy shifted coordinate values to the final `FractionalGridCoordinates`, except the NaN values.
                 * NaN values may exist if the given `point` has less dimensions than the grid geometry, in which
                 * case missing values have been replaced by `slice` values in the `target` array but not in the
                 * `coordinates` array. We want to keep the `slice` values in the `target` array.
                 *
                 * TODO: to be strict, we should skip the copy only if `slice.containsKey(i)` is true, because it
                 * could happen that a transform resulted in NaN values in other dimensions. But that check would
                 * be costly, so we avoid it for now.
                 */
                final double[] target = position.coordinates;
                for (int i = target.length; --i >= 0;) {
                    final double value = coordinates[i];
                    if (!Double.isNaN(value)) {
                        target[i] = value;
                    }
                }
            }
        }
        return position;
    }

    /**
     * Recomputes the {@link #inputToGrid} field. This method should be invoked when the transform
     * has not yet been computed or became outdated because {@link #inputCRS} needs to be changed.
     *
     * @param  crs  the new value to assign to {@link #inputCRS}.
     */
    private void setInputCRS(final CoordinateReferenceSystem crs)
            throws FactoryException, NoninvertibleTransformException
    {
        final GridGeometry gridGeometry = coverage.getGridGeometry();
        MathTransform gridToCRS = gridGeometry.getGridToCRS(PixelInCell.CELL_CENTER);
        MathTransform crsToGrid = gridToCRS.inverse();
        if (crs != null) {
            final CoordinateReferenceSystem stepCRS = coverage.getCoordinateReferenceSystem();
            final GeographicBoundingBox areaOfInterest = gridGeometry.geographicBBox();
            try {
                CoordinateOperation op = CRS.findOperation(crs, stepCRS, areaOfInterest);
                crsToGrid = MathTransforms.concatenate(op.getMathTransform(), crsToGrid);
            } catch (FactoryException main) {
                /*
                 * Above block tried to compute a "CRS to grid" transform in the most direct way.
                 * It covers the usual case where the point has the required number of dimensions,
                 * and works better if the point has more dimensions (extra dimensions are ignored).
                 * The following block covers the opposite case, where the point does not have enough
                 * dimensions. We try to fill missing dimensions with the help of the `slice` map.
                 */
                final Map<Integer,Long> slice = getDefaultSlice();
                try {
                    CoordinateOperation op = CRS.findOperation(stepCRS, crs, areaOfInterest);
                    gridToCRS = MathTransforms.concatenate(gridToCRS, op.getMathTransform());
                    final TransformSeparator ts = new TransformSeparator(gridToCRS);
                    final int  crsDim = gridToCRS.getTargetDimensions();
                    final int gridDim = gridToCRS.getSourceDimensions();
                    int[] mandatory = new int[gridDim];
                    int n = 0;
                    for (int i=0; i<gridDim; i++) {
                        if (!slice.containsKey(i)) {
                            mandatory[n++] = i;
                        }
                    }
                    mandatory = ArraysExt.resize(mandatory, n);
                    ts.addSourceDimensions(mandatory);          // Retain grid dimensions having no default value.
                    ts.setSourceExpandable(true);               // Retain more grid dimensions if they are required.
                    ts.addTargetDimensionRange(0, crsDim);      // Force retention of all CRS dimensions.
                    gridToCRS = ts.separate();
                    crsToGrid = gridToCRS.inverse();            // With less source dimensions, may be invertible now.
                    mandatory = ts.getSourceDimensions();       // Output grid dimensions computed by `crsToGrid`.
                    final int valueColumn = mandatory.length;   // Matrix column where to write default values.
                    final MatrixSIS m = Matrices.createZero(gridDim+1, valueColumn+1);
                    m.setElement(gridDim, valueColumn, 1);
                    n = 0;
                    for (int j=0; j<gridDim; j++) {
                        if (Arrays.binarySearch(mandatory, j) >= 0) {
                            m.setElement(j, n++, 1);            // Computed value to pass through.
                        } else {
                            final Long value = slice.get(j);
                            if (value == null) {
                                final GridExtent extent = gridGeometry.extent;
                                throw new FactoryException(Resources.format(Resources.Keys.NoNDimensionalSlice_3,
                                                crsDim, extent.getAxisIdentification(j, j), extent.getSize(j)));
                            }
                            m.setElement(j, valueColumn, value);
                        }
                    }
                    crsToGrid = MathTransforms.concatenate(crsToGrid, MathTransforms.linear(m));
                } catch (RuntimeException | FactoryException | NoninvertibleTransformException ex) {
                    main.addSuppressed(ex);
                    throw main;
                }
            }
        }
        // Modify fields only after everything else succeeded.
        position    = new FractionalGridCoordinates.Position(crsToGrid.getTargetDimensions());
        inputCRS    = crs;
        inputToGrid = crsToGrid;
    }

    /**
     * Invoked when a recoverable exception occurred.
     * Those exceptions must be minor enough that they can be silently ignored in most cases.
     *
     * @param  caller     the method where exception occurred.
     * @param  exception  the exception that occurred.
     */
    private static void recoverableException(final String caller, final TransformException exception) {
        Logging.recoverableException(GridExtent.LOGGER, DefaultEvaluator.class, caller, exception);
    }
}
