Create function GDALInterpolateAtPoint

[jjimenezshaw]

What does this PR do?

Creates the function GDALInterpolateAtPoint to interpolate a value in a raster band.
Interpolation methods are (so far) nearest, bilinear and cubic.
It can be used via

• C API GDALRasterInterpolateAtPoint
• C++ API Band::InterpolateAtPoint
• python via InterpolateAtPoint in a raster band
• CLI app gdallocationinfo with the new option -r

The C API is

CPLErr GDALRasterInterpolateAtPoint(GDALRasterBandH hBand, double dfPixel,
                                    double dfLine,
                                    GDALRIOResampleAlg eInterpolation,
                                    double *pdfRealValue, double *pdfImagValue)

For complex layers the implementation is not done yet (to be done in a later PR), but the API is ready. In C++ pdfImagValue is by default a null_ptr to simplify the usage.

As sugested by @rouault , the implementation is extracted from the RPC code, that currently has a cache on the raster data. That makes the code more performant if several values are requested in the same area of the raster. The code is removed from alg/gdal_rpc.cpp and pasted with some modifications in alg/gdal_interpolateatpoint.cpp.

Interpolation methods

Only nearest, bilinear and cubic methods are implemented so far. For the borders a minor adjust is done to be consistent between the top and left borders (that were extrapolating) with the bottom and right (that are downgrading to a simpler algorithm). Now both are downgrading to a simpler algorithm.

However this produces discontinuities in the returned values at the borders. In a later PR I would like to change the behaviour in the borders doing something similar as show in the images at https://en.wikipedia.org/wiki/Bicubic_interpolation

Implementation problems

There are some questions/problems described in the comments.

What are related issues/pull requests?

https://lists.osgeo.org/pipermail/gdal-dev/2024-April/058899.html

Tasklist

• Make sure code is correctly formatted (cf pre-commit configuration)
• Add test case(s)
• Add documentation
• Updated Python API documentation (swig/include/python/docs/)
• Review
• Adjust for comments
• All CI builds and checks have passed

[jjimenezshaw]

Problem with cubic interpolation:
As mentioned in the python test test_interpolateatpoint_cubic_several_points

    raster_array_1 = np.array(
        (
            [1.0, 2.0, 1.5, -0.3],
            [1.0, 2.0, 1.5, -0.3],
            [1.0, 2.0, 1.5, -0.3],
            [1.0, 2.0, 1.5, -0.3],
        )
    )
    mem_ds.GetRasterBand(1).WriteArray(raster_array_1)

    # using https://tools.timodenk.com/cubic-spline-interpolation to get some samples
    got_cubic = mem_ds.GetRasterBand(1).InterpolateAtPoint(1.5, 1.5, gdal.GRIORA_Cubic)
    assert got_cubic == pytest.approx(
        1.75, 1e-6
    )  # I would expect 2.0, as it is a input point. To be studied

I would expect a value of 2.0, that is the pixel value at x=1.5. However it returns 1.75.
Debugging the interpolation code (copied from the RPC algorithm) I could see that all points have weights not null. Is that expected?

[jjimenezshaw]

Problem in python return code
As shown in the test

def test_interpolateatpoint_errors():

    ds = gdal.Open("data/byte.tif")

    res = ds.GetRasterBand(1).InterpolateAtPoint(1000, 120, gdal.GRIORA_Bilinear)
    assert res == 0.0  # I would expect None. To be studied

those coordinates are out of bound. I would expect it returning a None, but it return 0.0. I do not know what is happening in the swig code.

[jjimenezshaw]

Please, feel free to suggest new names or locations for functions and files if needed.

[rouault]

2 potential reasons:

• perhaps because the cubic spline in this page is not the same method as the cubic of GDAL. The Cubic in GDAL is a Catmull-Rohm cubic (

  gdal/alg/gdalwarpkernel.cpp

  Line 3391 in b6b5546

  static double GWKCubic(double dfX)

  )
• or the 0.5 pixel shift applied (cf comments "Convert from upper left corner of pixel coordinates to center of" in GDALRPCGetDEMHeight()). The reasoning is that the pixel value is exactly reached at the center of the pixel.

[jjimenezshaw]

I plotted the values with the three interpolations and got this:
(5 pixels in the x axis, the last two equal)
in any case, I will try to fix it in a later PR, because that interpolation affects RPC, and I do not want to add noise there.

def test_plot():
    import matplotlib.pyplot as plt

    mem_ds = gdal.GetDriverByName("MEM").Create(
        "", xsize=5, ysize=4, bands=1, eType=gdal.GDT_Float32
    )
    np = pytest.importorskip("numpy")
    raster_array_1 = np.array(
        (
            [1.0, 2.0, 1.5, 0.3, 0.3],
            [1.0, 2.0, 1.5, 0.3, 0.3],
            [1.0, 2.0, 1.5, 0.3, 0.3],
            [1.0, 2.0, 1.5, 0.3, 0.3],
        )
    )
    band = mem_ds.GetRasterBand(1)
    band.WriteArray(raster_array_1)

    t = np.arange(0.0, 5.0, 0.01)
    y = 1.5
    v_cubic = [band.InterpolateAtPoint(x, y, gdal.GRIORA_Cubic) for x in t]
    v_linear = [band.InterpolateAtPoint(x, y, gdal.GRIORA_Bilinear) for x in t]
    v_near = [band.InterpolateAtPoint(x, y, gdal.GRIORA_NearestNeighbour) for x in t]

    plt.plot(t, v_linear, "r", t, v_cubic, "b", t, v_near, "g")
    plt.grid(True)
    plt.show()

[jjimenezshaw]

• The Cubic in GDAL is a Catmull-Rohm cubic

The code in BiCubicKernel seems to be similar to GWKBSpline

gdal/alg/gdalwarpkernel.cpp

Line 3464 in b6b5546

static double GWKBSpline(double x)

[rouault]

The code in BiCubicKernel seems to be similar to GWKBSpline

oh! you're right. I didn't realize that. This was introduced in 8a7c8dc . Ok, so we should rather map GRIORA_CubicSpline to it, and rename it internally to CubicBSpline
It would be nice to add proper mapping for Catmull-Rohm cubic / GRIORA_Cubic, but that can be done later

[jjimenezshaw]

The code in BiCubicKernel seems to be similar to GWKBSpline

oh! you're right. I didn't realize that. This was introduced in 8a7c8dc . Ok, so we should rather map GRIORA_CubicSpline to it, and rename it internally to CubicBSpline It would be nice to add proper mapping for Catmull-Rohm cubic / GRIORA_Cubic, but that can be done later

That was my impression yesterday night. But it opens questions respect to RPC code.
I see three enums about interpolation methods (are there more?):

• DEMResampleAlg
• GDALRIOResampleAlg
• GDALResampleAlg

The last two have to be equivalent if I understand correctly the comment before their definitions.
The first one is used only in RPC, right?
But based on what we have seen here, DEMResampleAlg::DRA_Cubic is equivalent to GDALResampleAlg::GRA_CubicSpline

At gdal_rpc.cpp is the only place I found that something is set to DRA_Cubic:

    const char *pszDEMInterpolation =
        CSLFetchNameValueDef(papszOptions, "RPC_DEMINTERPOLATION", "bilinear");
    if (EQUAL(pszDEMInterpolation, "near"))
    {
        psTransform->eResampleAlg = DRA_NearestNeighbour;
    }
    else if (EQUAL(pszDEMInterpolation, "bilinear"))
    {
        psTransform->eResampleAlg = DRA_Bilinear;
    }
    else if (EQUAL(pszDEMInterpolation, "cubic"))
    {
        psTransform->eResampleAlg = DRA_Cubic;
    }

that is "conflicting" with other places (where we use the -r argument) about the meaning of cubic.

So, me guess about what to do is:

• Implement an interpolation in InterpolateAtPoint for GDALResampleAlg::GRA_Cubic for cubic and rename the current one for cubicspline.
• When RPC is using DRA_Cubic, use GRA_CubicSpline. That will keep the current (3.9) behavior. Would it be a good idea to change that in GDAL 4.0, so RPC would have cubic and cubicspline, consistent with other places in the code?
• Change the API of the internal function GDALInterpolateAtPoint to use GDALRIOResampleAlg or GDALResampleAlg instead of DEMResampleAlg...
• Question: What is the difference between GDALRIOResampleAlg and GDALResampleAlg? Which one should I use in the internal function GDALInterpolateAtPoint?

[rouault]

The last two have to be equivalent if I understand correctly the comment before their definitions.

The values of the constants are designed to match when they are available in both, but some methods are present in one and not in the other one

The first one is used only in RPC, right?

yes was purely an implementation detail of the RPC code

But based on what we have seen here, DEMResampleAlg::DRA_Cubic is equivalent to GDALResampleAlg::GRA_CubicSpline

yes apparently. so DRA_Cubic should be probably renamed to DRA_CubicSpline for consistence

that is "conflicting" with other places (where we use the -r argument) about the meaning of cubic.

indeed. For existing RPC_DEMINTERPOLATION=cubic, we can't change the name. But for the new -r flag of gdallocationinfo, we should accept "cubicspline" for now, and reserve "cubic" for Catmull-Rohm cubic if we implement it at some later time

Change the API of the internal function GDALInterpolateAtPoint to use GDALRIOResampleAlg

yes we could do that to avoid adding new enumerations.

What is the difference between GDALRIOResampleAlg and GDALResampleAlg?

GDALRIOResampleAlg is for the methods supported by RasterIO() and the overview building code
GDALResampleAlg is for the warping code

Which one should I use in the internal function GDALInterpolateAtPoint?

I would go for GDALRIOResampleAlg

[coveralls]

coverage: 69.298% (-0.003%) from 69.301%
when pulling 3b089b2 on jjimenezshaw:interpolate-at-point
into 15866aa on OSGeo:master.

[rouault]

Problem in python return code

I've pushed a commit in your branch to fix/workaround that.

You didn't do anything wrong. Back in 2006, the IF_ERROR_RETURN_NONE typemap was half-restored in 86d4f1d after having been removed in 1c6e19f . So the end result is that the IF_ERROR_RETURN_NONE typemap currently only eats the CPLErr return code but doesn't do anything with it. I had always that feeling that something didn't work as intended when using GetStatistics() which uses this typemap too, but never investigated in depth. I'm a bit hesitant if the typemap should be fixed now. Probably yes. Anyway this is out of scope for this PR.

[jjimenezshaw]

I hope I addressed all the comments. Including make obvious that the interpolation is splinecubic. Next PR will be to include cubic.
Let me squash the commits.