Linear interpolation vectorized for faster calculation matrixwise

Question

To calibrate large amounts of picture-Data in Python, I have to apply a correction elementwise to every pixel. Therefore, I've got a matrix with 200x200x2x9 values. The 200x200 representing the pixelposition, the 2x9 are 9 pairs of values, representing a base and a corresponding correction value for this particular pixel which need to be linear interpolated, polynominal interpolation does not work very well. For example this looks like this:

Now I have a large amount of pictures and to save time I want to apply the correction of every pixel as fast as possible. I was thinking of a vectorized function. The input value from the picture should be a correction based on the position within the multiple points of the interpolation.

A for loop costs me a huge amount of time.

Can somebody help me with that?

***Edit

I.e., for a single pixel, how do you wish to combine a scalar pixel value and the 2x9 interpolation matrix. What is the shape of your result, and how would the interpolation work?

Every pixel in the picture has a value between 0 and 16383. I want to find out, where in the series of linear interpolations this value is placed. In the example shown above 9 pairs of grey-value and the corresponding correction. So if the pixel has e.g. a value of ~6000, the calculated correction is -15. After the calculation, I'd like to have a 200x200 matrix with values (in this example) between ±15.

Answer 1

We start by defining some test code so that we have a minimum reproducible example. For simplicity I use the same base points for the interpolation per pixel. However, to test the robustness, I define the range so that we will have pixels below the lowest base point and ones above the highest.

image = np.random.randint(0, 16383, (200, 200), dtype='i2')

correction = np.empty((200, 200, 2, 9), dtype='i2')
correction[..., 0, :] = np.linspace(
    0, 16383, num=11, endpoint=False, dtype='i2')[1:-1]
correction[..., 1, :] = np.random.randint(-16, 17, (200, 200, 9), dtype='i2')

To make indexing simpler, we flatten the image to 1D and separate the components from the correction matrix.

origshape = image.shape
image = image.ravel()
correction = correction.reshape(image.shape + correction.shape[2:])
bases = correction[:, 0].T
basevalues = correction[:, 1].T

Now we want to find the value range per pixel.

inrange = (image >= bases[:-1]) & (image <= bases[1:])
belowrange = image < bases[0]
aboverange = image > bases[-1]

Then we pick the corresponding values from the correction. We can use np.choose or np.select for this. This step could use some tuning so that the approach works with fancy-indexing. The results will be invalid for values that are above or below the base points. We correct this at a later step. Note that the lower and upper bases default to 0 and 1 for these cases. That avoids issues with the interpolation in these cases.

lowerbases = np.select(inrange, bases[:-1])
upperbases = np.select(inrange, bases[1:], default=1)
lowervalues = np.select(inrange, basevalues[:-1])
uppervalues = np.select(inrange, basevalues[1:])

Now we apply a plain old linear interpolation.

upperpercent = (image - lowerbases).astype('f4') / (upperbases - lowerbases)
lowerpercent = 1. - upperpercent
correctionvalues = (lowerpercent * lowervalues
                    + upperpercent * uppervalues).astype('i2')

For the values out of range, we can apply whatever condition we want. I choose to correct with the lowest or highest value. You may choose to extrapolate instead.

correctionvalues[belowrange] = basevalues[0, belowrange]
correctionvalues[aboverange] = basevalues[-1, aboverange]

Now all that remains is applying the values and turning the image back into 2D

image += correctionvalues
image = image.reshape(origshape)