pythonopencvinterpolationremapspatial-interpolation

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.


Solution

  • 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()
    

    result: pyramid with jagged edges

    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()
    

    result: clearly piecewise-constant function

    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()
    

    pretty smooth pyramid

    Much better.