5

I've an image of about 8000x9000 size as a numpy matrix. I also have a list of indices in a numpy 2xn matrix. These indices are fractional as well as may be out of image size. I need to interpolate the image and find the values for the given indices. If the indices fall outside, I need to return numpy.nan for them. Currently I'm doing it in for loop as below

def interpolate_image(image: numpy.ndarray, indices: numpy.ndarray) -> numpy.ndarray:
    """

    :param image:
    :param indices: 2xN matrix. 1st row is dim1 (rows) indices, 2nd row is dim2 (cols) indices
    :return:
    """
    # Todo: Vectorize this
    M, N = image.shape
    num_indices = indices.shape[1]
    interpolated_image = numpy.zeros((1, num_indices))
    for i in range(num_indices):
        x, y = indices[:, i]
        if (x < 0 or x > M - 1) or (y < 0 or y > N - 1):
            interpolated_image[0, i] = numpy.nan
        else:
            # Todo: Do Bilinear Interpolation. For now nearest neighbor is implemented
            interpolated_image[0, i] = image[int(round(x)), int(round(y))]
    return interpolated_image

But the for loop is taking huge amount of time (as expected). How can I vectorize this? I found scipy.interpolate.interp2d, but I'm not able to use it. Can someone explain how to use this or any other method is also fine. I also found this, but again it is not according to my requirements. Given x and y indices, these generated interpolated matrices. I don't want that. For the given indices, I just want the interpolated values i.e. I need a vector output. Not a matrix.

I tried like this, but as said above, it gives a matrix output

f = interpolate.interp2d(numpy.arange(image.shape[0]), numpy.arange(image.shape[1]), image, kind='linear')
interp_image_vect = f(indices[:,0], indices[:,1])
RuntimeError: Cannot produce output of size 73156608x73156608 (size too large)

For now, I've implemented nearest-neighbor interpolation. scipy interp2d doesn't have nearest neighbor. It would be good if the library function as nearest neighbor (so I can compare). If not, then also fine.

Nagabhushan S N
  • 6,407
  • 8
  • 44
  • 87

1 Answers1

2

It looks like scipy.interpolate.RectBivariateSpline will do the trick:

from scipy.interpolate import RectBivariateSpline
image = # as given
indices = # as given

spline = RectBivariateSpline(numpy.arange(M), numpy.arange(N), image)

interpolated = spline(indices[0], indices[1], grid=False)

This gets you the interpolated values, but it doesn't give you nan where you need it. You can get that with where:

nans = numpy.zeros(interpolated.shape) + numpy.nan
x_in_bounds = (0 <= indices[0]) & (indices[0] < M)
y_in_bounds = (0 <= indices[1]) & (indices[1] < N)
bounded = numpy.where(x_in_bounds & y_in_bounds, interpolated, nans)

I tested this with a 2624x2624 image and 100,000 points in indices and all told it took under a second.

Nathan Vērzemnieks
  • 5,495
  • 1
  • 11
  • 23
  • 1
    Thank you so much. In the meantime, my code changed a bit. I'm doing shifting of origin and stuff. Right now, directly applying your code, gave all NaN's. I'll try to set it properly and check. I'll accept the answer after that. Thanks a ton! – Nagabhushan S N Mar 30 '19 at 13:18
  • Somehow the `nan`s in the original image are causing problem. `interpolated` is all `nan`s. Instead, if I replace all the `nan`s in image with `0` or `-1`, then it works. `image1 = numpy.where(numpy.isnan(image), -1, image)` Any idea why? Or how to use `nan`s in this interpolation? – Nagabhushan S N Mar 30 '19 at 18:10
  • It's not clear what it would mean to interpolate with `nan` - it is, after all, not a number. I think it makes sense to replace the `nan` values in the original image with zero, although it may depend on why they're there. – Nathan Vērzemnieks Mar 30 '19 at 18:37
  • I'm implementing Panorama. So, I need to superimpose images. At every pixel, if 2 or more images are present, then I need to take their mean. Now each image size can vary. So, what I did was to put all images in a sufficiently large matrices. At a pixel, if image is not present, I'll represent it with `nan` so that I can do `numpy.nanmean` later. – Nagabhushan S N Mar 30 '19 at 20:22