Method grid_apply_internal

src/transformations/gridshift.cpp:606–723

#define TOL 1e-12

PJ_XYZ gridshiftData::grid_apply_internal(
    PJ_CONTEXT *ctx, const std::string &type, bool isVerticalOnly,
    const PJ_XYZ in, PJ_DIRECTION direction, const GenericShiftGrid *grid,
    GenericShiftGridSet *gridset, bool &shouldRetry) {

    shouldRetry = false;
    if (in.x == HUGE_VAL)
        return in;

    /* normalized longitude of input */
    const NS_PROJ::ExtentAndRes *extent;
    PJ_XY normalized_in = normalizeX(grid, in, extent);

    bool biquadraticInterpolationOut = false;
    PJ_XYZ shift = grid_interpolate(ctx, type, normalized_in, grid,
                                    biquadraticInterpolationOut);
    if (grid->hasChanged()) {
        shouldRetry = gridset->reopen(ctx);
        PJ_XYZ out;
        out.x = out.y = out.z = HUGE_VAL;
        return out;
    }
    if (shift.x == HUGE_VAL)
        return shift;

    if (direction == PJ_FWD) {
        PJ_XYZ out = in;
        out.x += shift.x;
        out.y += shift.y;
        out.z += shift.z;
        return out;
    }

    if (isVerticalOnly) {
        PJ_XYZ out = in;
        out.z -= shift.z;
        return out;
    }

    PJ_XY guess;
    guess.x = normalized_in.x - shift.x;
    guess.y = normalized_in.y - shift.y;

    // NOAA NCAT transformer tool doesn't do iteration in the reverse path.
    // Do the same (only for biquadratic, although NCAT applies this logic to
    // bilinear too)
    // Cf
    // https://github.com/noaa-ngs/ncat-lib/blob/77bcff1ce4a78fe06d0312102ada008aefcc2c62/src/gov/noaa/ngs/grid/Transformer.java#L374
    // When trying to do iterative reverse path with biquadratic, we can
    // get convergence failures on points that are close to the boundary of
    // cells or half-cells. For example with
    // echo -122.4250009683  37.8286740788 0 | bin/cct -I +proj=gridshift
    // +grids=tests/us_noaa_nadcon5_nad83_1986_nad83_harn_conus_extract_sanfrancisco.tif
    // +interpolation=biquadratic
    if (!biquadraticInterpolationOut) {
        int i = MAX_ITERATIONS;
        const double toltol = TOL * TOL;
        PJ_XY diff;