Search code examples
c++rasterosrmbilinear-interpolation

Bilinear Interpolation - OSRM Rastersource

I've got a question about bilinear interpolation in the OSRM-Project. I understand the "normal" bilinear interpolation. Here the picture from Wikipedia, what is insane:

Bilinear Interpolation - Wikipedia

Now I'm trying to understand the bilinear interpolation which is used in the OSRM-Project for raster source data.

// Query raster source using bilinear interpolation
RasterDatum RasterSource::GetRasterInterpolate(const int lon, const int lat) const
{
if (lon < xmin || lon > xmax || lat < ymin || lat > ymax)
{
    return {};
}

const auto xthP = (lon - xmin) / xstep;
const auto ythP =
    (ymax - lat) /
    ystep; // the raster texture uses a different coordinate system with y pointing downwards

const std::size_t top = static_cast<std::size_t>(fmax(floor(ythP), 0));
const std::size_t bottom = static_cast<std::size_t>(fmin(ceil(ythP), height - 1));
const std::size_t left = static_cast<std::size_t>(fmax(floor(xthP), 0));
const std::size_t right = static_cast<std::size_t>(fmin(ceil(xthP), width - 1));

// Calculate distances from corners for bilinear interpolation
const float fromLeft = xthP - left; // this is the fraction part of xthP
const float fromTop = ythP - top;   // this is the fraction part of ythP
const float fromRight = 1 - fromLeft;
const float fromBottom = 1 - fromTop;

return {static_cast<std::int32_t>(raster_data(left, top) * (fromRight * fromBottom) +
                                  raster_data(right, top) * (fromLeft * fromBottom) +
                                  raster_data(left, bottom) * (fromRight * fromTop) +
                                  raster_data(right, bottom) * (fromLeft * fromTop))};
}

Original Code here

Can someone explain me how the code works?

The input format are the SRTM data in ASCII format.

The variables height and width are defined as nrows and ncolumns. The variables xstep and ystep are defined as:

return (max - min) / (static_cast<float>(count) - 1)

Where count is height for ystep and width for xstep, max and min similar.

And another question: Can I use the same code for data in TIF-format and the whole world?

Solution

  • Horizontal pixel coordinates are in the range [0, width - 1]; similarly vertical coordinates are in [0, height - 1]. (Zero-indexing convention used in many many languages including C++)

    The lines

    const auto xthP = (lon - xmin) / xstep; (and for ythP)

    Convert the input image-space coordinates (long, lat) into pixel coordinates. xstep is the width of each pixel in image-space.

    Rounding this down (using floor) gives pixels intersected by the sample area on one side, and rounding up (ceil) gives the pixels on the other side. For the X-coordinate these give left and right.

    The reason for using fmin and fmax are to clamp the coordinates so that they don't exceed the pixel coordinate range.

    EDIT: since you are trying to interpret this picture, I'll list the corresponding parts below:

    • Q11 = (left, top)
    • Q12 - (left, bottom), etc.
    • P = (xthP, ythP)
    • R1 = fromTop, R2 = fromBottom etc.

    A good start point would be http://www.cs.uu.nl/docs/vakken/gr/2011/Slides/06-texturing.pdf, slide 27. In future though, Google is your friend.