/*
 * Licensed to the Apache Software Foundation (ASF) under one or more
 * contributor license agreements.  See the NOTICE file distributed with
 * this work for additional information regarding copyright ownership.
 * The ASF licenses this file to You under the Apache License, Version 2.0
 * (the "License"); you may not use this file except in compliance with
 * the License.  You may obtain a copy of the License at
 *
 *     http://www.apache.org/licenses/LICENSE-2.0
 *
 * Unless required by applicable law or agreed to in writing, software
 * distributed under the License is distributed on an "AS IS" BASIS,
 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 * See the License for the specific language governing permissions and
 * limitations under the License.
 */
package org.apache.sis.coverage.grid;

import java.util.Arrays;
import java.util.Objects;
import java.util.function.Supplier;
import javax.measure.Quantity;
import javax.measure.quantity.Length;
import org.opengis.geometry.Envelope;
import org.opengis.geometry.DirectPosition;
import org.opengis.util.FactoryException;
import org.opengis.referencing.datum.PixelInCell;
import org.opengis.referencing.crs.CoordinateReferenceSystem;
import org.opengis.referencing.operation.CoordinateOperation;
import org.opengis.referencing.operation.MathTransform;
import org.opengis.referencing.operation.TransformException;
import org.apache.sis.util.ArraysExt;
import org.apache.sis.util.collection.BackingStoreException;
import org.apache.sis.referencing.CRS;
import org.apache.sis.referencing.util.CoordinateOperations;
import org.apache.sis.referencing.util.WraparoundApplicator;
import org.apache.sis.metadata.iso.extent.DefaultGeographicBoundingBox;
import org.apache.sis.referencing.operation.transform.MathTransforms;
import org.apache.sis.referencing.operation.transform.WraparoundTransform;
import org.apache.sis.geometry.AbstractDirectPosition;
import org.apache.sis.geometry.Envelopes;
import org.apache.sis.image.ImageProcessor;
import org.apache.sis.measure.Quantities;
import org.apache.sis.measure.Units;
import org.apache.sis.util.logging.Logging;
import org.apache.sis.util.internal.Numerics;


/**
 * Finds a transform from grid cells in a source coverage to geospatial positions in the CRS of a target coverage.
 * This class differs from {@link CRS#findOperation CRS#findOperation(…)} because of the gridded aspect of inputs.
 * {@linkplain GridGeometry grid geometries} give more information about how referencing is applied on datasets.
 * With them, this class provides three additional benefits:
 *
 * <ul class="verbose">
 *   <li>Detect dimensions where target coordinates are constrained to constant values because
 *     the {@linkplain GridExtent#getSize(int) grid size} is only one cell in those dimensions.
 *     This is an important because it makes possible to find operations normally impossible:
 *     we can still produce an operation to a target CRS even if some dimensions have no corresponding
 *     source CRS.</li>
 *
 *   <li>Detect how to handle wraparound axes. For example if a raster spans 150° to 200° of longitude,
 *     this class understands that -170° of longitude should be translated to 190° for thar particular
 *     raster. This will work even if the minimum and maximum values declared in the longitude axis do
 *     not match that range.</li>
 *
 *   <li>Use the area of interest and grid resolution for refining the coordinate operation between two CRS.</li>
 * </ul>
 *
 * <b>Note:</b> except for the {@link #gridToGrid()} convenience method,
 * this class does not provide the complete chain of operations from grid to grid.
 * It provides only the operation from <em>cell indices</em> in source grid to <em>coordinates in the CRS</em>
 * of destination grid. Callers must add the last step (conversion from target CRS to cell indices) themselves.
 *
 * @author  Alexis Manin (Geomatys)
 * @author  Martin Desruisseaux (Geomatys)
 */
final class CoordinateOperationFinder implements Supplier<double[]> {
    /**
     * Whether the operation is between cell centers or cell corners.
     *
     * @see #setAnchor(PixelInCell)
     */
    private PixelInCell anchor;

    /**
     * The grid geometry which is the source/target of the coordinate operation to find.
     */
    private final GridGeometry source, target;

    /**
     * The target coordinate values, computed only if needed. This is computed by {@link #get()}, which is
     * itself invoked indirectly by {@link org.apache.sis.referencing.operation.CoordinateOperationFinder}.
     * The value is cached in case {@code get()} is invoked many times during the same finder execution.
     *
     * @see #get()
     */
    private double[] coordinates;

    /**
     * The coordinate operation from source to target CRS, computed when first needed.
     * May be {@code null} if there is no information about the source and target CRS.
     * Note that identity operation is not equivalent to {@code null} because identity operation will still
     * be checked for wraparound axes (because the CRS is known), while null operation will have no check.
     *
     * @see #knowChangeOfCRS
     * @see #changeOfCRS()
     */
    private CoordinateOperation changeOfCRS;

    /**
     * Whether the {@link #changeOfCRS} operation has been determined.
     * Note that result of determining that operation may be {@code null}, which is why we need this flag.
     */
    private boolean knowChangeOfCRS;

    /**
     * The {@link #changeOfCRS} transform together with {@link WraparoundTransform} if needed.
     * The wraparound is used for handling images crossing the anti-meridian.
     *
     * @see #gridToCRS()
     */
    private MathTransform forwardChangeOfCRS;

    /**
     * Inverse of {@link #changeOfCRS} transform together with {@link WraparoundTransform} if needed.
     * The wraparound is used for handling images crossing the anti-meridian.
     *
     * <p>Contrarily to {@link #forwardChangeOfCRS}, the process that determine this {@code inverseChangeOfCRS}
     * transform should check if wraparound is really needed. This is because {@code inverseChangeOfCRS} will be
     * used much more extensively (for every pixels) than other transforms.</p>
     *
     * @see #isWraparoundNeeded
     * @see #isWraparoundApplied
     * @see #inverse()
     */
    private MathTransform inverseChangeOfCRS;

    /**
     * Transform from “grid coordinates of the source” to “geospatial coordinates of the target”.
     * This is the concatenation of {@link #source} "grid to CRS" with {@link #forwardChangeOfCRS},
     * possibly with wraparound handling and cached for reuse by {@link #inverse()}:
     *
     * {@snippet lang="java" :
     *     forwardChangeOfCRS = changeOfCRS.getMathTransform();
     *     // + wraparound handling if applicable.
     *     gridToCRS = source.getGridToCRS(anchor);
     *     gridToCRS = MathTransforms.concatenate(gridToCRS, forwardChangeOfCRS);
     *     }
     *
     * @see #gridToCRS()
     */
    private MathTransform gridToCRS;

    /**
     * Transform from the target CRS to the source grid, with {@link WraparoundTransform} applied if needed.
     * This is the concatenation of {@link #inverseChangeOfCRS} with inverse of {@link #source} "grid to CRS",
     * possibly with wraparound handling:
     *
     * {@snippet lang="java" :
     *     inverseChangeOfCRS = forwardChangeOfCRS.inverse();
     *     // + wraparound handling if applicable.
     *     crsToGrid = gridToCRS.inverse();
     *     crsToGrid = MathTransforms.concatenate(inverseChangeOfCRS, crsToGrid);
     *     }
     *
     * @see #inverse()
     * @see #applyWraparound(MathTransform)
     */
    private MathTransform crsToGrid;

    /**
     * Whether the {@link #isWraparoundNeeded} value has been determined. This flag controls whether to perform
     * a more extensive check of wraparound occurrence. This flag should be {@code false} the first time that
     * {@link #inverse()} is invoked and {@code true} the next time.
     */
    private boolean isWraparoundNeedVerified;

    /**
     * Whether {@link #inverseChangeOfCRS} needs to include a {@link WraparoundTransform} step. We do this check
     * only for {@link #inverseChangeOfCRS} because it is the transform which will be executed for every pixels.
     * By contrast, {@link #forwardChangeOfCRS} will systematically contain a {@link WraparoundTransform} step
     * because we use it only for transforming envelopes and for the test that determines the value of
     * this {@code isWraparoundNeeded} flag.
     */
    private boolean isWraparoundNeeded;

    /**
     * Whether {@link WraparoundTransform} has been applied on {@link #inverseChangeOfCRS}. This field complements
     * {@link #isWraparoundNeeded} because a delay may exist between the time we detected that wraparound is needed
     * and the time we applied the necessary operation steps.
     *
     * <p>Note that despite this field name, a {@code true} value does not imply that {@link #inverseChangeOfCRS}
     * and {@link #crsToGrid} transforms really contain some {@link WraparoundTransform} steps. It only means that
     * the {@link WraparoundApplicator#forDomainOfUse WraparoundApplicator.forDomainOfUse(…)} method has been invoked.
     * That method may have decided to not insert any wraparound steps.</p>
     *
     * @see #applyWraparound(MathTransform)
     */
    private boolean isWraparoundApplied;

    /**
     * Whether to disable completely all wraparounds checks.
     * If {@code true}, then calculation done in this class should be equivalent to following code:
     *
     * {@snippet lang="java" :
     *     forwardChangeOfCRS = changeOfCRS.getMathTransform();
     *     inverseChangeOfCRS = forwardChangeOfCRS.inverse();
     *     gridToCRS          = source.getGridToCRS(anchor);
     *     crsToGrid          = gridToCRS.inverse();
     *     gridToCRS          = MathTransforms.concatenate(gridToCRS, forwardChangeOfCRS);
     *     crsToGrid          = MathTransforms.concatenate(inverseChangeOfCRS, crsToGrid);
     *     }
     *
     * <b>Tip:</b> searching usage of this field should help to identify code doing wraparound handling.
     *
     * @see #nowraparound()
     */
    private boolean isWraparoundDisabled;

    /**
     * Creates a new finder initialized to {@link PixelInCell#CELL_CORNER} anchor.
     *
     * @param  source  the grid geometry which is the source of the coordinate operation to find.
     * @param  target  the grid geometry which is the target of the coordinate operation to find.
     */
    CoordinateOperationFinder(final GridGeometry source, final GridGeometry target) {
        this.source = source;
        this.target = target;
        this.anchor = PixelInCell.CELL_CORNER;
    }

    /**
     * Verifies the presence of a CRS considered mandatory,
     * unless the CRS of the opposite grid is also missing.
     *
     * @param  rs  {@code true} is source CRS is mandatory, {@code false} if target CRS is mandatory.
     * @throws IncompleteGridGeometryException if a mandatory CRS is missing.
     */
    final void verifyPresenceOfCRS(final boolean rs) {
        if ((rs ? target : source).isDefined(GridGeometry.CRS)) {
            Objects.requireNonNull((rs ? source : target).getCoordinateReferenceSystem());
        }
    }

    /**
     * Sets whether operations will be between cell centers or cell corners.
     * This method must be invoked before any other method in this class.
     * The {@link PixelInCell#CELL_CORNER} value should be used first
     * in order to cache values computed relative to pixel corners.
     *
     * @param  newValue  whether operations will be between cell centers or cell corners.
     */
    final void setAnchor(final PixelInCell newValue) {
        /*
         * If `coordinates` is non-null, it means that `CoordinateOperationContext.getConstantCoordinates()`
         * has been invoked during previous `gridToCRS()` execution, which implies that `changeOfCRS` depends
         * on the `anchor` value. In such case we will need to recompute all fields that depend on `changeOfCRS`.
         */
        anchor = newValue;
        gridToCRS = null;
        crsToGrid = null;
        if (coordinates != null) {
            coordinates        = null;
            changeOfCRS        = null;
            forwardChangeOfCRS = null;
            inverseChangeOfCRS = null;
            knowChangeOfCRS    = false;
            // Do not clear `isWraparoundNeeded`; its value is still valid.
        }
    }

    /**
     * Disables completely all wraparounds operation.
     *
     * @see #isWraparoundDisabled
     */
    final void nowraparound() {
        gridToCRS                = null;        // For forcing recomputation.
        crsToGrid                = null;
        forwardChangeOfCRS       = null;
        inverseChangeOfCRS       = null;
        isWraparoundNeeded       = false;
        isWraparoundApplied      = true;
        isWraparoundNeedVerified = true;
        isWraparoundDisabled     = true;
    }

    /**
     * Returns the target of the "corner to CRS" transform.
     * May be {@code null} if the neither the source and target grid geometry define a CRS.
     */
    final CoordinateReferenceSystem getTargetCRS() {
        return (changeOfCRS != null) ? changeOfCRS.getTargetCRS() :
                source.isDefined(GridGeometry.CRS) ? source.getCoordinateReferenceSystem() : null;
    }

    /**
     * Returns the coordinate operation from source CRS to target CRS. It may be the identity operation.
     * We try to take envelopes in account because the operation choice may depend on the geographic area.
     *
     * @todo Specify also the desired resolution, computed from target grid geometry. It will require
     *       more direct use of {@link org.apache.sis.referencing.operation.CoordinateOperationContext}.
     *       As a side effect, we could remove {@link CoordinateOperations#CONSTANT_COORDINATES} hack.
     *
     * @return operation from source CRS to target CRS, or {@code null} if a CRS is not specified.
     * @throws FactoryException if no operation can be found between the source and target CRS.
     * @throws TransformException if some coordinates cannot be transformed to the specified target.
     */
    private CoordinateOperation changeOfCRS() throws FactoryException, TransformException {
        if (!knowChangeOfCRS) {
            final Envelope sourceEnvelope = source.envelope;
            final Envelope targetEnvelope = target.envelope;
            try {
                CoordinateOperations.CONSTANT_COORDINATES.set(this);
                if (sourceEnvelope != null && targetEnvelope != null) {
                    changeOfCRS = Envelopes.findOperation(sourceEnvelope, targetEnvelope);
                }
                if (changeOfCRS == null && source.isDefined(GridGeometry.CRS) && target.isDefined(GridGeometry.CRS)) {
                    final CoordinateReferenceSystem sourceCRS = source.getCoordinateReferenceSystem();
                    /*
                     * Unconditionally create operation even if CRS are the same. A non-null operation trig
                     * the check for wraparound axes, which is necessary even if the transform is identity.
                     */
                    DefaultGeographicBoundingBox areaOfInterest = null;
                    if (sourceEnvelope != null || targetEnvelope != null) try {
                        areaOfInterest = new DefaultGeographicBoundingBox();
                        areaOfInterest.setBounds(targetEnvelope != null ? targetEnvelope : sourceEnvelope);
                    } catch (TransformException e) {
                        areaOfInterest = null;
                        recoverableException("changeOfCRS", e);
                    }
                    changeOfCRS = CRS.findOperation(sourceCRS, target.getCoordinateReferenceSystem(), areaOfInterest);
                }
            } catch (BackingStoreException e) {                         // May be thrown by getConstantCoordinates().
                throw e.unwrapOrRethrow(TransformException.class);
            } finally {
                CoordinateOperations.CONSTANT_COORDINATES.remove();
            }
            knowChangeOfCRS = true;
        }
        return changeOfCRS;
    }

    /**
     * Computes the transform from “grid coordinates of the source” to “grid coordinates of the target”.
     * This is a concatenation of {@link #gridToCRS()} with target "CRS to grid" transform.
     *
     * @return operation from source grid indices to target grid indices.
     * @throws FactoryException if no operation can be found between the source and target CRS.
     * @throws TransformException if some coordinates cannot be transformed to the specified target.
     * @throws IncompleteGridGeometryException if required CRS or a "grid to CRS" information is missing.
     */
    final MathTransform gridToGrid() throws FactoryException, TransformException {
        final MathTransform step1 = gridToCRS();
        final MathTransform step2 = target.getGridToCRS(anchor);
        if (step1.equals(step2)) {                                          // Optimization for a common case.
            return MathTransforms.identity(step1.getSourceDimensions());
        } else {
            return MathTransforms.concatenate(step1, step2.inverse());
        }
    }

    /**
     * Computes the transform from “grid coordinates of the source” to “geospatial coordinates of the target”.
     * It may be the identity operation. We try to take envelopes in account because the operation choice may
     * depend on the geographic area.
     *
     * <p>The transform returned by this method applies wraparound checks systematically on every axes having
     * wraparound range. This method does not verify whether those checks are needed (i.e. whether wraparound
     * can possibly happen). This is okay because this transform is used only for transforming envelopes;
     * it is not used for transforming pixel coordinates.</p>
     *
     * <h4>Implementation note</h4>
     * After invocation of this method, the following fields are valid:
     * <ul>
     *   <li>{@link #changeOfCRS} — cached for {@link #inverse()} usage.</li>
     *   <li>{@link #forwardChangeOfCRS} — cached for next invocation of this {@code gridToCRS()} method.</li>
     * </ul>
     *
     * @return operation from source grid indices to target geospatial coordinates.
     * @throws FactoryException if no operation can be found between the source and target CRS.
     * @throws TransformException if some coordinates cannot be transformed to the specified target.
     * @throws IncompleteGridGeometryException if required CRS or a "grid to CRS" information is missing.
     */
    final MathTransform gridToCRS() throws FactoryException, TransformException {
        if (gridToCRS == null) {
            /*
             * The following line may throw IncompleteGridGeometryException, which is desired because
             * if that transform is missing, we cannot continue (we have no way to guess it).
             */
            gridToCRS = source.getGridToCRS(anchor);
            final CoordinateOperation changeOfCRS = changeOfCRS();
            if (changeOfCRS != null) {
                /*
                 * At this point we have a "grid → source CRS" transform. Append a "source CRS → target CRS" transform,
                 * which may be identity. A wraparound may be applied for keeping target coordinates inside the expected
                 * target domain.
                 */
apply:          if (forwardChangeOfCRS == null) {
                    forwardChangeOfCRS = changeOfCRS.getMathTransform();
                    if (!isWraparoundDisabled) {
                        DirectPosition sourceMedian = median(source, forwardChangeOfCRS);
                        DirectPosition targetMedian = median(target, null);
                        if (targetMedian == null) {
                            if (sourceMedian == null) {
                                break apply;
                            }
                            targetMedian = sourceMedian;
                            sourceMedian = null;
                        }
                        final WraparoundApplicator ap = new WraparoundApplicator(sourceMedian,
                                targetMedian, changeOfCRS.getTargetCRS().getCoordinateSystem());
                        forwardChangeOfCRS = ap.forDomainOfUse(forwardChangeOfCRS);
                    }
                }
                gridToCRS = MathTransforms.concatenate(gridToCRS, forwardChangeOfCRS);
            }
        }
        return gridToCRS;
    }

    /**
     * Computes the transform from “geospatial coordinates of the target” to “grid coordinates of the source”.
     * This is similar to invoking {@link MathTransform#inverse()} on {@link #gridToCRS()}, except in the way
     * wraparounds are handled.
     *
     * @return operation from target geospatial coordinates to source grid indices.
     * @throws FactoryException if no operation can be found between the source and target CRS.
     * @throws TransformException if some coordinates cannot be transformed.
     */
    final MathTransform inverse() throws FactoryException, TransformException {
        final MathTransform sourceCrsToGrid = source.getGridToCRS(anchor).inverse();
        final CoordinateOperation changeOfCRS = changeOfCRS();
        if (changeOfCRS == null) {
            return sourceCrsToGrid;
        }
        if (inverseChangeOfCRS == null) {
            inverseChangeOfCRS = changeOfCRS.getMathTransform().inverse();
            if (!isWraparoundDisabled) {
                isWraparoundApplied = false;
                if (!isWraparoundNeedVerified) {
                    isWraparoundNeedVerified = true;
                    /*
                     * Need to compute transform with wraparound checks, but contrarily to `gridToCRS()` we do not want
                     * `WraparoundTransform` to be systematically inserted. This is for performance reasons, because the
                     * transform returned by this method will be applied on every pixels of destination image. We create
                     * both transforms with and without wraparound, and check if their results differ.
                     *
                     * We give precedence to corners specified by target extent. The corners specified by source extent
                     * are used only as a fallback if the target extent has not been specified, in which case we assume
                     * that caller will fallback on source extent transformed to target coordinates.  The target extent
                     * is preferred because it may cover only a sub-region of the source, or conversely it may be world.
                     * If smaller, wraparound may become useless (i.e. sub-region may not cross anti-meridian anymore).
                     * If larger with [-180 … +180]° longitude range, the use of source extent may fail to detect that
                     * a part of the raster need to be rendered on each side of the [-180 … +180]° range.
                     */
                    final MathTransform inverseNoWrap = inverseChangeOfCRS;
                    final MathTransform crsToGridNoWrap = MathTransforms.concatenate(inverseNoWrap, sourceCrsToGrid);
                    if (target.isDefined(GridGeometry.EXTENT | GridGeometry.CRS)) {
                        if (applyWraparound(sourceCrsToGrid)) {
                            isWraparoundNeeded = isWraparoundNeeded(target.getExtent(),
                                    target.getGridToCRS(anchor), crsToGridNoWrap, null);
                        }
                    } else if (source.isDefined(GridGeometry.EXTENT)) {
                        isWraparoundNeeded = isWraparoundNeeded(source.getExtent(),
                                gridToCRS(), crsToGridNoWrap, sourceCrsToGrid);
                    }
                    if (!isWraparoundNeeded) {
                        inverseChangeOfCRS = inverseNoWrap;     // Discard the transform that was applying wraparound.
                        crsToGrid = crsToGridNoWrap;
                    }
                }
                if (isWraparoundNeeded) {
                    applyWraparound(sourceCrsToGrid);           // Update `inverseChangeOfCRS` if possible.
                }
            }
        }
        /*
         * Here, `inverseChangeOfCRS` already contains the wraparound step if needed.
         */
        if (crsToGrid == null) {
            crsToGrid = MathTransforms.concatenate(inverseChangeOfCRS, sourceCrsToGrid);
        }
        return crsToGrid;
    }

    /**
     * Verifies whether wraparound is needed for a "CRS to grid" transform.
     * This method converts coordinates of all corners of a grid (source or target) to the target CRS,
     * then (potentially back) to the source grid. This method uses one transform applying wraparounds
     * and another transform without wraparounds. By checking whether grid coordinates are equal with
     * both transforms, we determine if wraparound is necessary or not.
     *
     * @param  extent           the grid extent which is providing all corners to project.
     * @param  extentToCRS      transform from {@code extent} to target CRS.
     * @param  crsToGridNoWrap  inverse of {@code gridToCRS} but without handling of wraparound axes.
     * @param  sourceCrsToGrid  if {@code extent} is {@link #target} extent, shall be {@code null}.
     *                          If {@code extent} is {@link #source} extent, shall be the transform to
     *                          concatenate with {@link #inverseChangeOfCRS} for creating {@link #crsToGrid}.
     * @return whether wraparound transform seems needed.
     * @throws TransformException if an error occurred while transforming coordinates.
     */
    private boolean isWraparoundNeeded(final GridExtent extent, final MathTransform extentToCRS,
            final MathTransform crsToGridNoWrap, final MathTransform sourceCrsToGrid)
            throws FactoryException, TransformException
    {
        final boolean  mapCorner    = (anchor == PixelInCell.CELL_CORNER);
        final int      extentDim    = extent.getDimension();
        final int      gridDim      = crsToGridNoWrap.getTargetDimensions();
        final double[] buffer       = new double[Math.max(extentToCRS.getTargetDimensions(), gridDim)];
        final double[] reference    = new double[Math.max(extentDim, gridDim)];
        final double[] withoutWrap  = new double[gridDim];
        long maskOfUppers = Numerics.bitmask(extentDim);
        while (--maskOfUppers != 0) {
            for (int i=0; i<extentDim; i++) {
                final long bit = 1L << i;
                long cc;                                                // Grid coordinate of a corner.
                if ((maskOfUppers & bit) == 0) {
                    cc = extent.getLow(i);
                } else {
                    cc = extent.getHigh(i);                             // Inclusive.
                    if (mapCorner && cc != Long.MAX_VALUE) cc++;        // Make exclusive.
                }
                reference[i] = cc;
            }
            /*
             * Transform corner from the extent to target CRS, then (potentially back) to source grid coordinates.
             * We do not concatenate the forward and inverse transforms because we do not want MathTransformFactory
             * to simplify the transformation chain (e.g. replacing "Mercator → Inverse of Mercator" by an identity
             * transform), because such simplification would erase wraparound effects.
             */
            extentToCRS.transform(reference, 0, buffer, 0, 1);              // To coordinates in target CRS.
            crsToGridNoWrap.transform(buffer, 0, withoutWrap, 0, 1);        // To coordinates in source grid.
            /*
             * The reference must be a corner in the `source` grid. If the given extent was from `target` grid,
             * convert to source grid coordinates by completing the "target → CRS → source" chain of transforms.
             * The `crsToGrid` transform includes the wraparound, contrarily to `crsToGridNoWrap` used above.
             */
            if (sourceCrsToGrid == null) {
                // `applyWraparound()` already invoked by caller.
                crsToGrid.transform(buffer, 0, reference, 0, 1);
            }
            /*
             * Compare coordinates without wraparound with the reference. If they differ, we may consider that
             * wraparounds are needed.
             */
            boolean isBufferTransformed = false;
            for (int i=0; i<gridDim; i++) {
                final double error = Math.abs(withoutWrap[i] - reference[i]);
                if (!(error <= 1)) {                                            // Use `!` for catching NaN.
                    if (sourceCrsToGrid == null) {
                        /*
                         * If the `extent` specified in argument was the target extent, then the `reference`
                         * corner has already been computed using a transform that include wraparound checks.
                         * We do not consider NaN in `reference` as a need for wraparound because NaN values
                         * may occur when an operation between two CRS reduces the number of dimensions.
                         */
                        if (!Double.isNaN(reference[i])) {
                            return true;
                        }
                    } else {
                        /*
                         * If the `extent` specified in argument was the source extent, then the `reference`
                         * corner has been computed with a transform that does not use `inverseChangeOfCRS`.
                         * So it may be worth to double-check with another calculation of the corner, this
                         * time including wraparound steps.
                         */
                        if (!isBufferTransformed) {
                            isBufferTransformed = true;
                            if (!applyWraparound(sourceCrsToGrid)) {
                                return false;
                            }
                            crsToGrid.transform(buffer, 0, buffer, 0, 1);
                        }
                        if (Math.abs(buffer[i] - reference[i]) < (error <= Double.MAX_VALUE ? error : 1)) {
                            return true;
                        }
                    }
                }
            }
        }
        return false;
    }

    /**
     * Inserts {@link WraparoundTransform} steps into {@link #inverseChangeOfCRS} transform if possible.
     * IF this method returns {@code true}, then the {@link #inverseChangeOfCRS} and {@link #crsToGrid}
     * fields have been updated to transforms applying wraparound.
     *
     * @param  sourceCrsToGrid  value of {@code source.getGridToCRS(anchor).inverse()}.
     * @return whether at least one wraparound step has been added.
     * @throws TransformException if some coordinates cannot be transformed.
     */
    private boolean applyWraparound(final MathTransform sourceCrsToGrid) throws FactoryException, TransformException {
        if (!isWraparoundApplied) {
            isWraparoundApplied = true;
            DirectPosition sourceMedian = median(source, null);
            DirectPosition targetMedian = median(target, inverseChangeOfCRS);
            if (sourceMedian == null) {
                sourceMedian = targetMedian;
                targetMedian = null;
            }
            if (sourceMedian != null) {
                final MathTransform inverseNoWrap = inverseChangeOfCRS;
                final WraparoundApplicator ap = new WraparoundApplicator(targetMedian,
                        sourceMedian, changeOfCRS().getSourceCRS().getCoordinateSystem());
                inverseChangeOfCRS = ap.forDomainOfUse(inverseNoWrap);
                if (inverseChangeOfCRS != inverseNoWrap) {
                    crsToGrid = MathTransforms.concatenate(inverseChangeOfCRS, sourceCrsToGrid);
                    return true;
                }
            }
        }
        return false;
    }

    /**
     * Returns the point of interest converted to the Coordinate Reference System.
     * If the grid does not define a point of interest or does not define a CRS,
     * then this method returns {@code null}.
     *
     * @param  grid         the source or target grid providing the point of interest.
     * @param  changeOfCRS  transform from source CRS to target CRS, or {@code null} if none.
     */
    private static DirectPosition median(final GridGeometry grid, final MathTransform changeOfCRS) throws TransformException {
        if (!grid.isDefined(GridGeometry.EXTENT | GridGeometry.GRID_TO_CRS)) {
            return null;
        }
        return new AbstractDirectPosition() {
            /** The coordinates, computed when first needed. */
            private double[] coordinates;

            /** Returns the number of dimensions. */
            @Override public int getDimension() {
                return coordinates().length;
            }

            /** Returns the coordinate tuple, computed when first needed. */
            @SuppressWarnings("ReturnOfCollectionOrArrayField")
            private double[] coordinates() {
                if (coordinates == null) try {
                    final double[] poi = grid.getExtent().getPointOfInterest(PixelInCell.CELL_CENTER);
                    MathTransform tr = grid.getGridToCRS(PixelInCell.CELL_CENTER);
                    if (changeOfCRS != null) {
                        tr = MathTransforms.concatenate(tr, changeOfCRS);
                    }
                    coordinates = new double[tr.getTargetDimensions()];
                    tr.transform(poi, 0, coordinates, 0, 1);
                } catch (TransformException e) {
                    throw new BackingStoreException(e);
                }
                return coordinates;
            }

            /**
             * Returns the median rounded to a value having an exact representation in base 2 using about 10 bits.
             * The intent is to reduce the risk of rounding errors with add/subtract operations.
             */
            @Override public double getOrdinate(final int i) {
                final double m = coordinates()[i];
                final int power = 10 - Math.getExponent(m);
                return Math.scalb(Math.rint(Math.scalb(m, power)), -power);
            }
        };
    }

    /**
     * Invoked when the target CRS has some dimensions that the source CRS does not have.
     * For example, this is invoked during the conversion from (<var>x</var>, <var>y</var>)
     * coordinates to (<var>x</var>, <var>y</var>, <var>t</var>). If constant values can
     * be given to the missing dimensions, than those values are returned. Otherwise this
     * method returns {@code null}.
     *
     * <p>The returned array has a length equals to the number of dimensions in the target CRS.
     * Only coordinates in dimensions without source (<var>t</var> in above example) will be used.
     * All other coordinate values will be ignored.</p>
     *
     * @see org.apache.sis.referencing.operation.CoordinateOperationContext#getConstantCoordinates()
     */
    @Override
    @SuppressWarnings("ReturnOfCollectionOrArrayField")
    public double[] get() {
        if (coordinates == null && target.isDefined(GridGeometry.EXTENT | GridGeometry.GRID_TO_CRS)) {
            final MathTransform tr = target.getGridToCRS(anchor);
            coordinates = new double[tr.getTargetDimensions()];
            double[] gc = new double[tr.getSourceDimensions()];
            Arrays.fill(gc, Double.NaN);
            final GridExtent extent = target.getExtent();
            for (int i=0; i<gc.length; i++) {
                final long low = extent.getLow(i);
                if (low == extent.getHigh(i)) {
                    gc[i] = low;
                }
            }
            /*
             * At this point, the only grid coordinates with finite values are the ones where the
             * grid size is one cell (i.e. conversion to target CRS can produce only one value).
             * After conversion with `getGidToCRS(anchor)`, the corresponding target dimensions
             * will have non-NaN coordinate values only if they do not depend on any dimension
             * other than the one having a grid size of 1.
             */
            try {
                tr.transform(gc, 0, coordinates, 0, 1);
            } catch (TransformException e) {
                throw new BackingStoreException(e);
            }
        }
        return coordinates;
    }

    /**
     * Configures the accuracy hints on the given processor.
     *
     * <h4>Prerequisite</h4>
     * This method assumes that {@link #gridToCRS()} or {@link #inverse()}
     * has already been invoked before this method.
     */
    final void setAccuracyOf(final ImageProcessor processor) {
        final double accuracy = CRS.getLinearAccuracy(changeOfCRS);
        if (accuracy > 0) {
            Length qm = Quantities.create(accuracy, Units.METRE);
            Quantity<?>[] hints = processor.getPositionalAccuracyHints();       // Array is already a copy.
            for (int i=0; i<hints.length; i++) {
                if (Units.isLinear(hints[i].getUnit())) {
                    hints[i] = qm;
                    qm = null;
                }
            }
            if (qm != null) {
                hints = ArraysExt.append(hints, qm);
            }
            processor.setPositionalAccuracyHints(hints);                        // Null elements will be ignored.
        }
    }

    /**
     * Invoked when an ignorable exception occurred.
     *
     * @param  caller  the method where the exception occurred.
     * @param  e       the ignorable exception.
     */
    private static void recoverableException(final String caller, final Exception e) {
        Logging.recoverableException(GridExtent.LOGGER, CoordinateOperationFinder.class, caller, e);
    }
}
