Ошибка интерполяции переназначения OpenCV?

Question

Я использую функцию переназначения opencv, чтобы отобразить изображение в другую систему координат. Тем не менее, мои первоначальные тесты показывают, что есть некоторые проблемы с интерполяцией. Здесь я приведу простой пример постоянного сдвига в 0,1 пикселя для изображения, которое равно 0 везде, но в положении [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])

дает

0.90625
0.09375

Тем не менее, я бы ожидал 0,9 и 0,1 вместо этого, учитывая метод линейной интерполяции. Я делаю что-то не так или это какая-то цифровая проблема? Существуют ли более точные алгоритмы переназначения?

Благодарю.

Answer 1

Хорошо поймал. На мой взгляд, ваши ожидания верны, о чем свидетельствует np.interp дающий 0.1 а также 0.9 ценности.

Давайте нарисуем пирамиду (интерполирующую в диапазон 49:51 квадратных пикселей):

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

Заметили что-нибудь? При построении сетки 200x200 на пирамиде видны очень заметные ступени. Давайте посмотрим на сечение нашего результата:

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

Как видите, результатом является кусочно-постоянная функция, то есть огромная ошибка дискретизации в выходных данных. Чтобы быть точным, мы видим шаги 0.03125 == 1/32 в результате.

Я подозреваю, что cv2.remap не предназначен для использования для субпиксельных манипуляций, но для более масштабного отображения из одной сетки в другую. Другой вариант заключается в том, что внутренняя точность была принесена в жертву улучшениям производительности. В любом случае, вы не сходите с ума: вы должны видеть 0.1 а также 0.9 в результате точной (би) линейной интерполяции.

Если вы не привязаны к openCV из-за других задач, это отображение, т.е. 2d-интерполяция, может выполняться с различными битами scipy.interpolateа именно его части сделаны для 2d интерполяции. Для вашего особого случая линейной интерполяции на регулярной сетке, scipy.interpolate.RegularGridInterpolator или что-то подобное может быть уместным.

Или даже лучше (но я еще не использовал этот подмодуль): scipy.ndimage.map_coordinates кажется, что именно то, что вы ищете:

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

Применительно к примеру пирамиды:

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

Намного лучше.