OpenCV remap interpolation error?

I'm using opencv remap function to map an image to another coordinate system. However, my initial tests indicate that there are some issues with the interpolation. Here, I give a simple example of a constant 0.1 pixel shift for a image that is 0 everywhere but at position [50,50].

import cv2
import numpy as np

prvs = np.zeros((100,80), dtype=np.float32)
prvs[50:51, 50:51] = 1.

grid_x, grid_y = np.meshgrid(np.arange(prvs.shape[1]), np.arange(prvs.shape[0]))
grid_y = grid_y.astype(np.float32)
grid_x = grid_x.astype(np.float32) + 0.1

prvs_remapped = cv2.remap(prvs, grid_x, grid_y, interpolation=cv2.INTER_LINEAR)

print(prvs_remapped[50,50])
print(prvs_remapped[50,49])

gives

0.90625
0.09375

However, I would expect 0.9 and 0.1 instead, given the linear interpolation method. Am I doing something wrong or is this some numeric issue? Are there any more precise remapping algorithms around?

Thanks.

Nice catch. Your expectations are correct in my opinion, as exemplified by np.interp giving 0.1 and 0.9 values.

Let's plot a pyramid (interpolating into the 49:51 square pixel range):

import numpy as np
import cv2
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

prvs = np.zeros((100,80), dtype=np.float32)
prvs[50:51, 50:51] = 1

lin = np.linspace(49,51,200)
grid_x,grid_y = np.meshgrid(lin,lin)
grid_x = grid_x.astype(np.float32)
grid_y = grid_y.astype(np.float32)
prvs_zoommapped = cv2.remap(prvs, grid_x, grid_y, interpolation=cv2.INTER_LINEAR)

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(grid_x,grid_y,prvs_zoommapped,cmap='viridis')
plt.show()

Notice anything off? With a plotting grid of 200x200, there are very visible steps on the pyramid. Let's take a look at the cross-section of our result:

fig,ax = plt.subplots()
ax.plot(prvs_zoommapped[100,:],'x-')
ax.grid('on')
plt.show()

As you can see, the result is a piece-wise constant function, i.e. there's huge discretization error in the output. To be precise, we see steps of 0.03125 == 1/32 in the result.

My suspicion is that cv2.remap is not meant to be used for sub-pixel manipulations, but for a larger-scale mapping from one grid to another. The other option is that internally precision has been sacrificed for performance improvements. Either way, you're not going crazy: you should be seeing 0.1 and 0.9 as the result of exact (bi)linear interpolation.

If you're not committed to openCV due to other tasks, this mapping i.e. 2d interpolation can be performed with various bits of scipy.interpolate, namely its parts made for 2d interpolation. For your special case of linear interpolation on a regular grid, scipy.interpolate.RegularGridInterpolator or something similar might be appropriate.

Or even better (but I haven't used this submodule yet): scipy.ndimage.map_coordinates seems like exactly what you're looking for:

from scipy import ndimage
ndimage.map_coordinates(prvs, [[50.1, 49.1], [50, 50]], order=1)
# output: array([ 0.89999998,  0.1       ], dtype=float32)

Applied to the pyramid example:

import numpy as np
import cv2
from scipy import ndimage
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

prvs = np.zeros((100,80), dtype=np.float32)
prvs[50:51, 50:51] = 1

lin = np.linspace(49,51,200)
grid_x,grid_y = np.meshgrid(lin,lin)
prvs_zoommapped = ndimage.map_coordinates(prvs, [grid_x, grid_y], order=1)

fig = plt.figure()
ax = fig.add_subplot(111,projection='3d')
ax.plot_surface(grid_x,grid_y,prvs_zoommapped,cmap='viridis')
plt.show()

Much better.