反转实值索引网格
Inverting a real-valued index grid
OpenCV 的 remap() 使用实值索引网格从使用双线性插值的图像中采样值网格，并且 returns 样本网格作为新图像。
准确地说，让：
A = an image
X = a grid of real-valued X coords into the image.
Y = a grid of real-valued Y coords into the image.
B = remap(A, X, Y)
那么对于所有的像素坐标i,j,
B[i, j] = A(X[i, j], Y[i, j])
其中圆括号符号A(x, y)表示使用双线性插值法使用浮点值坐标x和y.[=27=求解图像A的像素值]
我的问题是：给定索引网格 X、Y，我如何生成 "inverse grid" X^-1、Y^-1 这样：
X(X^-1[i, j], Y^-1[i, j]) = i
Y(X^-1[i, j], Y^-1[i, j]) = j
和
X^-1(X[i, j], Y[i, j]) = i
Y^-1(X[i, j], Y[i, j]) = j
对于所有整数像素坐标i, j?
FWIW，图像和索引图 X 和 Y 的形状相同。但是，索引映射 X 和 Y 没有先验结构。例如，它们不一定是仿射或刚性变换。它们甚至可能是不可逆的，例如如果 X, Y 将 A 中的多个像素映射到 B 中相同的精确像素坐标。我正在寻找一种方法的想法，如果存在的话，它将找到一个合理的逆映射。
解决方案不需要基于 OpenCV，因为我没有使用 OpenCV，而是另一个具有 remap() 实现的库。虽然欢迎任何建议，但我特别热衷于 "mathematically correct"，即如果我的地图 M 是完全可逆的，该方法应该在机器精度的一些小范围内找到完美的逆。
OpenCV没有任何标准的方法。
如果您正在寻找完整的 ready-to-use 解决方案，我不确定我是否可以提供帮助，但我至少可以描述我几年前用来完成此任务的方法。
首先，您应该创建与源图像尺寸相同的重映射图。我创建了更大尺寸的地图以简化插值，并在最后一步将它们裁剪到合适的尺寸。然后你应该用以前的重新映射地图中存在的值填充它们（不是那么困难：只需迭代它们并且如果地图坐标 x 和 y 位于图像的限制内，将它们的行和列作为新的 y 和 x，并放入旧的新地图的 x 和 y 列和行）。这是相当简单的解决方案，但它给出了相当好的结果。对于完美的一个，您应该使用插值方法和相邻像素将旧的 x 和 y 插值到整数值。
在此之后，您应该手动重新映射像素颜色，或者使用像素坐标完全填充您的重新映射地图并使用来自 OpenCV 的版本。
您将遇到相当具有挑战性的任务：您应该在空白区域插入像素。换句话说，您应该取最近的 non-zero 像素坐标的距离，并根据这些距离混合颜色（如果您重新映射颜色）或坐标（如果您继续进行完整地图计算）分数。其实线性插值也没有那么难，你甚至可以看看remap()在OpenCV github page中的实现。对于 NN 插值，它会简单得多 - 只需取最近邻居的 color/coordinate。
最后一项任务是将区域外推到重映射像素区域的边界之外。也可以参考OpenCV的算法
据我了解，您有一个原始图像和一个转换后的图像，并且您希望在不知情的情况下恢复已应用的转换的性质，但假设它是合理的，例如旋转或fish-eye扭曲。
我会尝试对图像进行阈值处理以将其转换为二进制图像，包括索引图像和普通图像。然后尝试识别对象。大多数映射将至少保留连通性和欧拉数，索引中最大的对象大多仍将是平原中最大的对象。
然后花点时间查看匹配的图像/索引对，看看是否可以删除平移、旋转和缩放。这会为您提供多个反向地图，然后您可以尝试将它们拼接在一起。 （如果变换不简单则很难，但是无法解决重构任何变换的一般问题）。
如果您的地图是从单应性 H 派生的，您可以反转 H 并直接使用 cv::initUndistortRectifyMap() 创建反转地图。
例如在 Python:
import numpy as np.
map_size = () # fill in your map size
H_inv = np.linalg.inv(H)
map1, map2 = cv2.initUndistortRectifyMap(cameraMatrix=np.eye(3), distCoeffs=np.zeros(5), R=H_inv, newCameraMatrix=np.eye(3), size=map_size, m1type=cv2.CV_32FC1)
OpenCV 文档说明 initUndistortRectifyMap()：
The function actually builds the maps for the inverse mapping
algorithm that is used by remap(). That is, for each pixel (u, v) in
the destination image, the function computes the corresponding
coordinates in the source image.
如果你刚刚给了地图，你必须自己做。
然而，新地图坐标的插值并不是微不足道的，因为一个像素的支持区域可能非常大。
这是一个简单的 Python 解决方案，它通过 point-to-point 映射来反转地图。这可能会留下一些未分配的坐标，而其他坐标将被更新几次。所以地图可能有漏洞
这是一个演示这两种方法的 Python 小程序：
import cv2
import numpy as np
def invert_maps(map_x, map_y):
assert(map_x.shape == map_y.shape)
rows = map_x.shape[0]
cols = map_x.shape[1]
m_x = np.ones(map_x.shape, dtype=map_x.dtype) * -1
m_y = np.ones(map_y.shape, dtype=map_y.dtype) * -1
for i in range(rows):
for j in range(cols):
i_ = round(map_y[i, j])
j_ = round(map_x[i, j])
if 0 <= i_ < rows and 0 <= j_ < cols:
m_x[i_, j_] = j
m_y[i_, j_] = i
return m_x, m_y
def main():
img = cv2.imread("pigeon.png", cv2.IMREAD_GRAYSCALE)
# a simply rotation by 45 degrees
H = np.array([np.sin(np.pi/4), -np.cos(np.pi/4), 0, np.cos(np.pi/4), np.sin(np.pi/4), 0, 0, 0, 1]).reshape((3,3))
H_inv = np.linalg.inv(H)
map_size = (img.shape[1], img.shape[0])
map1, map2 = cv2.initUndistortRectifyMap(cameraMatrix=np.eye(3), distCoeffs=np.zeros(5), R=H, newCameraMatrix=np.eye(3), size=map_size, m1type=cv2.CV_32FC1)
map1_inv, map2_inv = cv2.initUndistortRectifyMap(cameraMatrix=np.eye(3), distCoeffs=np.zeros(5), R=H_inv, newCameraMatrix=np.eye(3), size=map_size, m1type=cv2.CV_32FC1)
map1_simple_inv, map2_simple_inv = invert_maps(map1, map2)
img1 = cv2.remap(src=img, map1=map1, map2=map2, interpolation=cv2.INTER_LINEAR)
img2 = cv2.remap(src=img1, map1=map1_inv, map2=map2_inv, interpolation=cv2.INTER_LINEAR)
img3 = cv2.remap(src=img1, map1=map1_simple_inv, map2=map2_simple_inv,
interpolation=cv2.INTER_LINEAR)
cv2.imshow("Original image", img)
cv2.imshow("Mapped image", img1)
cv2.imshow("Mapping forth and back with H_inv", img2)
cv2.imshow("Mapping forth and back with invert_maps()", img3)
cv2.waitKey(0)
if __name__ == '__main__':
main()
OP在这里。我想我找到了答案。我还没有实现它，如果有人想出了一个不那么繁琐的解决方案（或者发现这个有问题），我会选择他们的答案。
问题陈述
设A为源图像，B为目标图像，M为A坐标到B坐标的映射，即：
B[k, l, :] == A(M[k, l, 0], M[k, l, 1], :)
for all k, l in B's coords.
...其中方括号表示使用整数索引进行数组查找，圆括号表示使用 floating-point 索引进行双线性插值查找。我们使用更经济的符号重申以上内容：
B = A(M)
我们希望找到一个尽可能最好地将 B 映射回 A 的逆映射 N:
Find N s.t. A \approx B(N)
问题可以不参考A或B来陈述：
Find N = argmin_N || M(N) - I_n ||
...其中||*||表示Frobenius范数，I_n是与N维度相同的恒等映射，即映射其中：
I_n[i, j, :] == [i, j]
for all i, j
天真的解决方案
如果M的值都是整数，并且M是同构的，那么可以直接构造N为：
N[M[k, l, 0], M[k, l, 1], :] = [k, l]
for all k, l
或者在我们的简化符号中：
N[M] = I_m
...其中 I_m 是与 M 具有相同维度的恒等映射。
有两个问题：
- M 不是同构，因此对于不在 M.
- M 的值是 floating-point 坐标 [i, j]，而不是整数坐标。我们不能简单地为 bilinearly-interpolated 个数量 N(i, j, :) 赋值，因为 float-valued i, j。为了达到相同的效果，我们必须改为设置 [i, j] 的四个周围角的值 N[floor(i), floor(j), :], N[floor(i), ceil(j), :], N[ceil(i), floor(j), :], N[ceil(i), ceil(j), :] 这样插值 N(i, j, :) 等于期望值 [ k, l], 对于所有像素映射 [i, j] --> [k, l] in M.
解决方案
将空 N 构造为浮点数的 3D 张量：
N = zeros(size=(A.shape[0], A.shape[1], 2))
对于A的坐标space中的每个坐标[i,j]，做：
- 在 M 中找到 [i, j] 所在的 A-coordinates 的 2x2 网格。
计算将那些 A-coordinates 映射到它们对应的 B-coordinates 的单应矩阵 H（由 2x2 网格的像素索引给出）。
- 设 N[i, j, :] = matmul(H, [i, j])
此处可能代价高昂的步骤是在步骤 1 中搜索 M 中 A-coordinates 的 2x2 网格，该网格环绕 [i, j]。 brute-force 搜索会使整个算法复杂度为 O(n*m)，其中 n 是 A 中的像素数，m 是 B 中的像素数。
为了将其减少到 O(n)，可以改为 运行 每个 A-coordinate 四边形内的扫描线算法来识别它包含的所有 integer-valued 坐标 [i, j] .这可以预先计算为一个 hashmap，它将 integer-valued A 坐标 [i, j] 映射到其环绕四边形的 B 坐标 [k, l] 的 upper-left 角。
好吧，我必须自己解决这个 重映射反转问题，我将概述我的解决方案。
给定 X，Y 用于执行以下操作的 remap() 函数：
B[i, j] = A(X[i, j], Y[i, j])
我计算了 Xinv，Yinv 可以被 remap() 函数用来 反转 过程：
A[x, y] = B(Xinv[x,y],Yinv[x,y])
首先，我在 GitHub 上构建了一个 KD-Tree for the 2D point set {(X[i,j],Y[i,j]} so I can efficiently find the N nearest neighbors to a given point (x,y). I use Euclidian distance for my distance metric. I found a great C++ header lib for KD-Trees。
然后我遍历 A 网格中的所有 (x,y) 值，并在我的点集中找到 N = 5 最近的邻居 {(X[i_k,j_k],Y[i_k,j_k]) | k = 0 .. N-1}。
If distance d_k == 0 for some k then Xinv[x,y] = i_k and Yinv[x,y] = j_k, otherwise...
使用Inverse Distance Weighting (IDW)计算内插值：
- 让体重
w_k = 1 / pow(d_k, p)（我用p = 2）
Xinv[x,y] = (sum_k w_k * i_k)/(sum_k w_k) Yinv[x,y] = (sum_k w_k * j_k)/(sum_k w_k)
请注意，如果 B 是 W x H 图像，则 X 和 Y 是 W x H 浮点数组。如果 A 是 w x h 图像，则 Xinv 和 Yinv 是 w x h 浮点数组。与图像和地图大小保持一致很重要。
很有魅力！我的第一个版本尝试了暴力搜索，我什至从未等待它完成。我切换到 KD-Tree 然后我开始获得合理的 运行 次。如果我有时间，我想将其添加到 OpenCV。
下面的第二张图片使用 remap() 消除了第一张图片的镜头畸变。第三张图片是反转过程的结果。
您可以在已知点反转地图并将其插入到新网格中。
它会很好地工作，而失真不是很大。
这是在 Python 中使用 scipy.interpolate.griddata 的非常简单的实现：
map_x, map_y = cv2.initUndistortRectifyMap(K, D, None, new_K, image_size, cv2.CV_32FC1)
points = np.stack([map_x.flatten(), map_y.flatten()], axis=1)
grid = np.mgrid[:map_x.shape[0], :map_y.shape[1]]
values = grid.reshape(2, -1).T[..., ::-1]
from scipy.interpolate import griddata
grid_y, grid_x = grid
map_back = griddata(points, values, (grid_x, grid_y), method='cubic').astype(map_undistort.dtype)
如果对地图使用CV_32FC2，可以简化点的构造：
map_undistort, _ = cv2.initUndistortRectifyMap(K, D, None, new_K, image_size, cv2.CV_32FC2)
points = map_undistort.reshape(-1, 2)
这是@wcochran 的回答的一个实现。我试图恢复由 lensfunpy 导致的镜头校正。
mod = lensfunpy.Modifier(lens, cam.crop_factor, width, height)
mod.initialize(focal_length, aperture, distance)
undist_coords = mod.apply_geometry_distortion()
## the lens correction part
# im_undistorted = cv2.remap(im, undist_coords, None, cv2.INTER_CUBIC)
# im_undistorted = cv2.remap(im, undist_coords, None, cv2.INTER_LANCZOS4)
# cv2.imwrite(undistorted_image_path, im_undistorted)
undist_coords_f = undist_coords.reshape((-1, 2))
tree = KDTree(undist_coords_f)
def calc_val(point_pos):
nearest_dist, nearest_ind = tree.query([point_pos], k=5)
if nearest_dist[0][0] == 0:
return undist_coords_f[nearest_ind[0][0]]
# starts inverse distance weighting
w = np.array([1.0 / pow(d, 2) for d in nearest_dist])
sw = np.sum(w)
# embed()
x_arr = np.floor(nearest_ind[0] / 1080)
y_arr = (nearest_ind[0] % 1080)
xx = np.sum(w * x_arr) / sw
yy = np.sum(w * y_arr) / sw
return (xx, yy)
un_correction_x = np.zeros((720, 1080))
un_correction_y = np.zeros((720, 1080))
## reverse the lens correction
for i in range(720):
print("row %d operating" % i)
for j in range(1080):
un_correction_x[i][j], un_correction_y[i][j] = calc_val((i, j))
# print((i, j), calc_val((j, i)))
dstMap1, dstMap2 = cv2.convertMaps(un_correction_x.astype(np.float32), un_correction_y.astype(np.float32), cv2.CV_32FC2)
im_un_undistorted = cv2.remap(im_undistorted, dstMap1, dstMap2, cv2.INTER_LANCZOS4)
这是一个重要的问题，我很惊讶没有在任何标准库中更好地解决这个问题（至少据我所知）。
我对公认的解决方案不满意，因为它没有使用转换的隐式平滑度。我可能会错过重要的案例，但我无法想象映射在任何有用的意义上都是可逆的，但在像素尺度上是不平滑的。
平滑意味着不需要计算最近邻：最近的点是那些在原始网格上已经很近的点。
我的解决方案使用的事实是，在原始映射中，正方形 [(i,j), (i+1, j), (i+1, j+1), (i, j+1 )] 映射到内部没有其他点的四边形 [(X[i,j], Y[i,j], X[i+1,j], Y[i+1,j], ...] .然后逆映射只需要在四边形内插值。为此我使用逆双线性插值，这将在顶点和任何其他仿射变换给出精确的结果。
除了 numpy 之外，该实现没有其他依赖项。逻辑是 运行 遍历所有四边形并逐步构建反向映射。我把代码复制到这里，希望有足够的评论让这个想法足够清晰。
关于不太明显的东西的一些评论：
- 反双线性函数通常 return 坐标仅在 [0,1] 范围内。我删除了裁剪操作，因此超出范围的值意味着坐标在四边形之外（这是解决多边形点问题的一种扭曲方式！）。为了避免丢失边缘上的点，我实际上允许 [0,1] 范围之外的点，这通常意味着两个相邻的四边形可能会拾取一个索引。在这些罕见的情况下，我只是让结果成为两个结果的平均值，相信超出范围的点是以合理的方式“外推”的。
- 一般来说，所有的四边形都有不同的形状，它们与规则网格的重叠可以从无到有变化很多点。该例程一次求解所有四边形（利用
bilinear_inverse 的矢量化性质，但在每次迭代中仅选择坐标（到其边界框的偏移量）有效的四边形。
import numpy as np
def bilinear_inverse(p, vertices, numiter=4):
"""
Compute the inverse of the bilinear map from the unit square
[(0,0), (1,0), (1,1), (0,1)]
to the quadrilateral vertices = [p0, p1, p2, p4]
Parameters:
----------
p: array of shape (2, ...)
Points on which the inverse transforms are applied.
vertices: array of shape (4, 2, ...)
Coordinates of the vertices mapped to the unit square corners
numiter:
Number of Newton interations
Returns:
--------
s: array of shape (2, ...)
Mapped points.
This is a (more general) python implementation of the matlab implementation
suggested in
"""
p = np.asarray(p)
v = np.asarray(vertices)
sh = p.shape[1:]
if v.ndim == 2:
v = np.expand_dims(v, axis=tuple(range(2, 2 + len(sh))))
# Start in the center
s = .5 * np.ones((2,) + sh)
s0, s1 = s
for k in range(numiter):
# Residual
r = v[0] * (1 - s0) * (1 - s1) + v[1] * s0 * (1 - s1) + v[2] * s0 * s1 + v[3] * (1 - s0) * s1 - p
# Jacobian
J11 = -v[0, 0] * (1 - s1) + v[1, 0] * (1 - s1) + v[2, 0] * s1 - v[3, 0] * s1
J21 = -v[0, 1] * (1 - s1) + v[1, 1] * (1 - s1) + v[2, 1] * s1 - v[3, 1] * s1
J12 = -v[0, 0] * (1 - s0) - v[1, 0] * s0 + v[2, 0] * s0 + v[3, 0] * (1 - s0)
J22 = -v[0, 1] * (1 - s0) - v[1, 1] * s0 + v[2, 1] * s0 + v[3, 1] * (1 - s0)
inv_detJ = 1. / (J11 * J22 - J12 * J21)
s0 -= inv_detJ * (J22 * r[0] - J12 * r[1])
s1 -= inv_detJ * (-J21 * r[0] + J11 * r[1])
return s
def invert_map(xmap, ymap, diagnostics=False):
"""
Generate the inverse of deformation map defined by (xmap, ymap) using inverse bilinear interpolation.
"""
# Generate quadrilaterals from mapped grid points.
quads = np.array([[ymap[:-1, :-1], xmap[:-1, :-1]],
[ymap[1:, :-1], xmap[1:, :-1]],
[ymap[1:, 1:], xmap[1:, 1:]],
[ymap[:-1, 1:], xmap[:-1, 1:]]])
# Range of indices possibly within each quadrilateral
x0 = np.floor(quads[:, 1, ...].min(axis=0)).astype(int)
x1 = np.ceil(quads[:, 1, ...].max(axis=0)).astype(int)
y0 = np.floor(quads[:, 0, ...].min(axis=0)).astype(int)
y1 = np.ceil(quads[:, 0, ...].max(axis=0)).astype(int)
# Quad indices
i0, j0 = np.indices(x0.shape)
# Offset of destination map
x0_offset = x0.min()
y0_offset = y0.min()
# Index range in x and y (per quad)
xN = x1 - x0 + 1
yN = y1 - y0 + 1
# Shape of destination array
sh_dest = (1 + x1.max() - x0_offset, 1 + y1.max() - y0_offset)
# Coordinates of destination array
yy_dest, xx_dest = np.indices(sh_dest)
xmap1 = np.zeros(sh_dest)
ymap1 = np.zeros(sh_dest)
TN = np.zeros(sh_dest, dtype=int)
# Smallish number to avoid missing point lying on edges
epsilon = .01
# Loop through indices possibly within quads
for ix in range(xN.max()):
for iy in range(yN.max()):
# Work only with quads whose bounding box contain indices
valid = (xN > ix) * (yN > iy)
# Local points to check
p = np.array([y0[valid] + ix, x0[valid] + iy])
# Map the position of the point in the quad
s = bilinear_inverse(p, quads[:, :, valid])
# s out of unit square means p out of quad
# Keep some epsilon around to avoid missing edges
in_quad = np.all((s > -epsilon) * (s < (1 + epsilon)), axis=0)
# Add found indices
ii = p[0, in_quad] - y0_offset
jj = p[1, in_quad] - x0_offset
ymap1[ii, jj] += i0[valid][in_quad] + s[0][in_quad]
xmap1[ii, jj] += j0[valid][in_quad] + s[1][in_quad]
# Increment count
TN[ii, jj] += 1
ymap1 /= TN + (TN == 0)
xmap1 /= TN + (TN == 0)
if diagnostics:
diag = {'x_offset': x0_offset,
'y_offset': y0_offset,
'mask': TN > 0}
return xmap1, ymap1, diag
else:
return xmap1, ymap1
这是一个测试例子
import cv2 as cv
from scipy import ndimage as ndi
# Simulate deformation field
N = 500
sh = (N, N)
t = np.random.normal(size=sh)
dx = ndi.gaussian_filter(t, 40, order=(0,1))
dy = ndi.gaussian_filter(t, 40, order=(1,0))
dx *= 30/dx.max()
dy *= 30/dy.max()
# Test image
img = np.zeros(sh)
img[::10, :] = 1
img[:, ::10] = 1
img = ndi.gaussian_filter(img, 0.5)
# Apply forward mapping
yy, xx = np.indices(sh)
xmap = (xx-dx).astype(np.float32)
ymap = (yy-dy).astype(np.float32)
warped = cv.remap(img, xmap, ymap ,cv.INTER_LINEAR)
plt.imshow(warped, cmap='gray')
# Now invert the mapping
xmap1, ymap1 = invert_map(xmap, ymap)
unwarped = cv.remap(warped, xmap1.astype(np.float32), ymap1.astype(np.float32) ,cv.INTER_LINEAR)
plt.imshow(unwarped, cmap='gray')
迭代求解
上面的许多解决方案对我都不起作用，当贴图不可逆时失败，或者速度不是很快。
我提出了另一种 6 行迭代解决方案。
def invert_map(F):
I = np.zeros_like(F)
I[:,:,1], I[:,:,0] = np.indices(sh)
P = np.copy(I)
for i in range(10):
P += I - cv.remap(F, P, None, interpolation=cv.INTER_LINEAR)
return P
效果如何？
对于我为航空摄影反转地形校正图的用例，这种方法在 10 个步骤中轻松收敛到 1/10 像素。它也非常快，因为所有繁重的计算都隐藏在 OpenCV
中
它是如何工作的？
该方法使用的思想是，如果 (x', y') = F(x, y) 是一个映射，那么只要 F 的梯度很小，就可以用 (x, y) = -F(x', y') 来近似逆。
我们可以继续完善我们的映射，上面得到了我们的第一个预测（我是一个“身份映射”）：
G_1 = I - F
我们的第二个预测可以改编自：
G_2 = G_1 + I - F(G_1)
等等：
G_n+1 = G_n + I - F(G_n)
证明 G_n 收敛到逆 F^-1 很难，但我们可以轻松证明的是，如果 G 已经收敛，它将保持收敛。
假设G_n = F^-1，那么我们可以代入：
G_n+1 = G_n + I - F(G_n)
然后得到：
G_n+1 = F^-1 + I - F(F^-1)
G_n+1 = F^-1 + I - I
G_n+1 = F^-1
Q.E.D.
测试脚本
import cv2 as cv
from scipy import ndimage as ndi
import numpy as np
from matplotlib import pyplot as plt
# Simulate deformation field
N = 500
sh = (N, N)
t = np.random.normal(size=sh)
dx = ndi.gaussian_filter(t, 40, order=(0,1))
dy = ndi.gaussian_filter(t, 40, order=(1,0))
dx *= 10/dx.max()
dy *= 10/dy.max()
# Test image
img = np.zeros(sh)
img[::10, :] = 1
img[:, ::10] = 1
img = ndi.gaussian_filter(img, 0.5)
# Apply forward mapping
yy, xx = np.indices(sh)
xmap = (xx-dx).astype(np.float32)
ymap = (yy-dy).astype(np.float32)
warped = cv.remap(img, xmap, ymap ,cv.INTER_LINEAR)
plt.imshow(warped, cmap='gray')
def invert_map(F: np.ndarray):
I = np.zeros_like(F)
I[:,:,1], I[:,:,0] = np.indices(sh)
P = np.copy(I)
for i in range(10):
P += I - cv.remap(F, P, None, interpolation=cv.INTER_LINEAR)
return P
# F: The function to invert
F = np.zeros((sh[0], sh[1], 2), dtype=np.float32)
F[:,:,0], F[:,:,1] = (xmap, ymap)
# Test the prediction
unwarped = cv.remap(warped, invert_map(F), None, cv.INTER_LINEAR)
plt.imshow(unwarped, cmap='gray')
一种方法是获取原始地图，遍历其条目并获取 x 和 y 值的下限和上限。这给出了 (x,y), (xf,yf), (xc,yf), (xf,yc), (x c,yc) 在原始源图像的坐标中。然后，您可以填充一个结构，其中每一个都作为包含像素值和权重的索引，并对这些数据使用您喜欢的不规则网格插值。
这很容易通过反距离插值来实现，因为结构可以是图像数组累加，而权重是标量。 F 是原始源，G 是变形后的图像，F' 是恢复后的图像。地图是M.
将 F' 初始化为 0。创建一个与 F' 大小相同的浮点数的 0 初始化权重数组 W。
遍历 M。对于 M 中的每个整数，找到 4 个整数对及其与 (x,y) 的距离。从G中取出对应的像素值，按其距离的倒数加权，累加到F' like
F'(xf|c,yf|c)+=G(i,j)/sqrt((x-xf|c)^2+(y-yf|c)^2)
然后把那个权重累加到
W(xf|c,yf|c)+=1./sqrt((x-xf|c)^2+(y-yf|c)^2).
完成后，通过迭代对 F' 进行归一化，并将每个像素除以其在 W 中的对应条目（如果它不为零）。
在这一点上，图像通常接近完成，但是由于高下采样率，F' 中的一些像素可能没有被填充。所以你在 W 中来回传递几次以找到 0 权重条目，并从它们的非空邻居中插入这些像素。这部分也可以使用 KNN 搜索和插值来完成，因为它们通常不多。
与 KNN 方法相比，它易于实现并且缩放性更好（尽管我认为这对于小图像非常有用）。缺点是反距离不是最好的插值方案，但如果映射不是太块并且原始没有被大量下采样，它似乎工作得很好。当然，如果下采样率很高，你就不得不推断出很多丢失的信息，所以它本质上会给出粗略的结果。
如果你想尽可能多地从地图反演中挤出来，你可以尝试求解由原始插值方案定义的（可能欠定的）方程组；并非不可能，但具有挑战性。
KNNRegressor 具有反转网格映射的所有必要组件！
给你：
from sklearn.neighbors import KNeighborsRegressor
def get_inverse_maps(map1, map2):
regressor = KNeighborsRegressor(3)
X = np.concatenate((map2[..., None], map1[..., None]), axis=-1).reshape(-1, 2)
y = np.indices(map1.shape).transpose((1, 2, 0)).reshape(-1, 2)
regressor.fit(X, y)
map_inv = regressor.predict(y).reshape(map1.shape + (2,)).astype(np.float32)
map_inv2, map_inv1 = map_inv[..., 0], map_inv[..., 1]
return map_inv1, map_inv2
OpenCV 的 remap() 使用实值索引网格从使用双线性插值的图像中采样值网格，并且 returns 样本网格作为新图像。
准确地说，让：
A = an image
X = a grid of real-valued X coords into the image.
Y = a grid of real-valued Y coords into the image.
B = remap(A, X, Y)
那么对于所有的像素坐标i,j,
B[i, j] = A(X[i, j], Y[i, j])
其中圆括号符号A(x, y)表示使用双线性插值法使用浮点值坐标x和y.[=27=求解图像A的像素值]
我的问题是：给定索引网格 X、Y，我如何生成 "inverse grid" X^-1、Y^-1 这样：
X(X^-1[i, j], Y^-1[i, j]) = i
Y(X^-1[i, j], Y^-1[i, j]) = j
和
X^-1(X[i, j], Y[i, j]) = i
Y^-1(X[i, j], Y[i, j]) = j
对于所有整数像素坐标i, j?
FWIW，图像和索引图 X 和 Y 的形状相同。但是，索引映射 X 和 Y 没有先验结构。例如，它们不一定是仿射或刚性变换。它们甚至可能是不可逆的，例如如果 X, Y 将 A 中的多个像素映射到 B 中相同的精确像素坐标。我正在寻找一种方法的想法，如果存在的话，它将找到一个合理的逆映射。
解决方案不需要基于 OpenCV，因为我没有使用 OpenCV，而是另一个具有 remap() 实现的库。虽然欢迎任何建议，但我特别热衷于 "mathematically correct"，即如果我的地图 M 是完全可逆的，该方法应该在机器精度的一些小范围内找到完美的逆。
OpenCV没有任何标准的方法。
如果您正在寻找完整的 ready-to-use 解决方案，我不确定我是否可以提供帮助，但我至少可以描述我几年前用来完成此任务的方法。
首先，您应该创建与源图像尺寸相同的重映射图。我创建了更大尺寸的地图以简化插值，并在最后一步将它们裁剪到合适的尺寸。然后你应该用以前的重新映射地图中存在的值填充它们（不是那么困难：只需迭代它们并且如果地图坐标 x 和 y 位于图像的限制内，将它们的行和列作为新的 y 和 x，并放入旧的新地图的 x 和 y 列和行）。这是相当简单的解决方案，但它给出了相当好的结果。对于完美的一个，您应该使用插值方法和相邻像素将旧的 x 和 y 插值到整数值。
在此之后，您应该手动重新映射像素颜色，或者使用像素坐标完全填充您的重新映射地图并使用来自 OpenCV 的版本。
您将遇到相当具有挑战性的任务：您应该在空白区域插入像素。换句话说，您应该取最近的 non-zero 像素坐标的距离，并根据这些距离混合颜色（如果您重新映射颜色）或坐标（如果您继续进行完整地图计算）分数。其实线性插值也没有那么难，你甚至可以看看remap()在OpenCV github page中的实现。对于 NN 插值，它会简单得多 - 只需取最近邻居的 color/coordinate。
最后一项任务是将区域外推到重映射像素区域的边界之外。也可以参考OpenCV的算法
据我了解，您有一个原始图像和一个转换后的图像，并且您希望在不知情的情况下恢复已应用的转换的性质，但假设它是合理的，例如旋转或fish-eye扭曲。
我会尝试对图像进行阈值处理以将其转换为二进制图像，包括索引图像和普通图像。然后尝试识别对象。大多数映射将至少保留连通性和欧拉数，索引中最大的对象大多仍将是平原中最大的对象。
然后花点时间查看匹配的图像/索引对，看看是否可以删除平移、旋转和缩放。这会为您提供多个反向地图，然后您可以尝试将它们拼接在一起。 （如果变换不简单则很难，但是无法解决重构任何变换的一般问题）。
如果您的地图是从单应性 H 派生的，您可以反转 H 并直接使用 cv::initUndistortRectifyMap() 创建反转地图。
例如在 Python:
import numpy as np.
map_size = () # fill in your map size
H_inv = np.linalg.inv(H)
map1, map2 = cv2.initUndistortRectifyMap(cameraMatrix=np.eye(3), distCoeffs=np.zeros(5), R=H_inv, newCameraMatrix=np.eye(3), size=map_size, m1type=cv2.CV_32FC1)
OpenCV 文档说明 initUndistortRectifyMap()：
The function actually builds the maps for the inverse mapping algorithm that is used by
remap(). That is, for each pixel (u, v) in the destination image, the function computes the corresponding coordinates in the source image.
如果你刚刚给了地图，你必须自己做。 然而，新地图坐标的插值并不是微不足道的，因为一个像素的支持区域可能非常大。
这是一个简单的 Python 解决方案，它通过 point-to-point 映射来反转地图。这可能会留下一些未分配的坐标，而其他坐标将被更新几次。所以地图可能有漏洞
这是一个演示这两种方法的 Python 小程序：
import cv2
import numpy as np
def invert_maps(map_x, map_y):
assert(map_x.shape == map_y.shape)
rows = map_x.shape[0]
cols = map_x.shape[1]
m_x = np.ones(map_x.shape, dtype=map_x.dtype) * -1
m_y = np.ones(map_y.shape, dtype=map_y.dtype) * -1
for i in range(rows):
for j in range(cols):
i_ = round(map_y[i, j])
j_ = round(map_x[i, j])
if 0 <= i_ < rows and 0 <= j_ < cols:
m_x[i_, j_] = j
m_y[i_, j_] = i
return m_x, m_y
def main():
img = cv2.imread("pigeon.png", cv2.IMREAD_GRAYSCALE)
# a simply rotation by 45 degrees
H = np.array([np.sin(np.pi/4), -np.cos(np.pi/4), 0, np.cos(np.pi/4), np.sin(np.pi/4), 0, 0, 0, 1]).reshape((3,3))
H_inv = np.linalg.inv(H)
map_size = (img.shape[1], img.shape[0])
map1, map2 = cv2.initUndistortRectifyMap(cameraMatrix=np.eye(3), distCoeffs=np.zeros(5), R=H, newCameraMatrix=np.eye(3), size=map_size, m1type=cv2.CV_32FC1)
map1_inv, map2_inv = cv2.initUndistortRectifyMap(cameraMatrix=np.eye(3), distCoeffs=np.zeros(5), R=H_inv, newCameraMatrix=np.eye(3), size=map_size, m1type=cv2.CV_32FC1)
map1_simple_inv, map2_simple_inv = invert_maps(map1, map2)
img1 = cv2.remap(src=img, map1=map1, map2=map2, interpolation=cv2.INTER_LINEAR)
img2 = cv2.remap(src=img1, map1=map1_inv, map2=map2_inv, interpolation=cv2.INTER_LINEAR)
img3 = cv2.remap(src=img1, map1=map1_simple_inv, map2=map2_simple_inv,
interpolation=cv2.INTER_LINEAR)
cv2.imshow("Original image", img)
cv2.imshow("Mapped image", img1)
cv2.imshow("Mapping forth and back with H_inv", img2)
cv2.imshow("Mapping forth and back with invert_maps()", img3)
cv2.waitKey(0)
if __name__ == '__main__':
main()
OP在这里。我想我找到了答案。我还没有实现它，如果有人想出了一个不那么繁琐的解决方案（或者发现这个有问题），我会选择他们的答案。
问题陈述
设A为源图像，B为目标图像，M为A坐标到B坐标的映射，即：
B[k, l, :] == A(M[k, l, 0], M[k, l, 1], :)
for all k, l in B's coords.
...其中方括号表示使用整数索引进行数组查找，圆括号表示使用 floating-point 索引进行双线性插值查找。我们使用更经济的符号重申以上内容：
B = A(M)
我们希望找到一个尽可能最好地将 B 映射回 A 的逆映射 N:
Find N s.t. A \approx B(N)
问题可以不参考A或B来陈述：
Find N = argmin_N || M(N) - I_n ||
...其中||*||表示Frobenius范数，I_n是与N维度相同的恒等映射，即映射其中：
I_n[i, j, :] == [i, j]
for all i, j
天真的解决方案
如果M的值都是整数，并且M是同构的，那么可以直接构造N为：
N[M[k, l, 0], M[k, l, 1], :] = [k, l]
for all k, l
或者在我们的简化符号中：
N[M] = I_m
...其中 I_m 是与 M 具有相同维度的恒等映射。
有两个问题：
- M 不是同构，因此对于不在 M.
- M 的值是 floating-point 坐标 [i, j]，而不是整数坐标。我们不能简单地为 bilinearly-interpolated 个数量 N(i, j, :) 赋值，因为 float-valued i, j。为了达到相同的效果，我们必须改为设置 [i, j] 的四个周围角的值 N[floor(i), floor(j), :], N[floor(i), ceil(j), :], N[ceil(i), floor(j), :], N[ceil(i), ceil(j), :] 这样插值 N(i, j, :) 等于期望值 [ k, l], 对于所有像素映射 [i, j] --> [k, l] in M.
解决方案
将空 N 构造为浮点数的 3D 张量：
N = zeros(size=(A.shape[0], A.shape[1], 2))
对于A的坐标space中的每个坐标[i,j]，做：
- 在 M 中找到 [i, j] 所在的 A-coordinates 的 2x2 网格。 计算将那些 A-coordinates 映射到它们对应的 B-coordinates 的单应矩阵 H（由 2x2 网格的像素索引给出）。
- 设 N[i, j, :] = matmul(H, [i, j])
此处可能代价高昂的步骤是在步骤 1 中搜索 M 中 A-coordinates 的 2x2 网格，该网格环绕 [i, j]。 brute-force 搜索会使整个算法复杂度为 O(n*m)，其中 n 是 A 中的像素数，m 是 B 中的像素数。
为了将其减少到 O(n)，可以改为 运行 每个 A-coordinate 四边形内的扫描线算法来识别它包含的所有 integer-valued 坐标 [i, j] .这可以预先计算为一个 hashmap，它将 integer-valued A 坐标 [i, j] 映射到其环绕四边形的 B 坐标 [k, l] 的 upper-left 角。
好吧，我必须自己解决这个 重映射反转问题，我将概述我的解决方案。
给定 X，Y 用于执行以下操作的 remap() 函数：
B[i, j] = A(X[i, j], Y[i, j])
我计算了 Xinv，Yinv 可以被 remap() 函数用来 反转 过程：
A[x, y] = B(Xinv[x,y],Yinv[x,y])
首先，我在 GitHub 上构建了一个 KD-Tree for the 2D point set {(X[i,j],Y[i,j]} so I can efficiently find the N nearest neighbors to a given point (x,y). I use Euclidian distance for my distance metric. I found a great C++ header lib for KD-Trees。
然后我遍历 A 网格中的所有 (x,y) 值，并在我的点集中找到 N = 5 最近的邻居 {(X[i_k,j_k],Y[i_k,j_k]) | k = 0 .. N-1}。
If distance
d_k == 0for somekthenXinv[x,y] = i_kandYinv[x,y] = j_k, otherwise...使用Inverse Distance Weighting (IDW)计算内插值：
- 让体重
w_k = 1 / pow(d_k, p)（我用p = 2） Xinv[x,y] = (sum_k w_k * i_k)/(sum_k w_k)Yinv[x,y] = (sum_k w_k * j_k)/(sum_k w_k)
- 让体重
请注意，如果 B 是 W x H 图像，则 X 和 Y 是 W x H 浮点数组。如果 A 是 w x h 图像，则 Xinv 和 Yinv 是 w x h 浮点数组。与图像和地图大小保持一致很重要。
很有魅力！我的第一个版本尝试了暴力搜索，我什至从未等待它完成。我切换到 KD-Tree 然后我开始获得合理的 运行 次。如果我有时间，我想将其添加到 OpenCV。
下面的第二张图片使用 remap() 消除了第一张图片的镜头畸变。第三张图片是反转过程的结果。
您可以在已知点反转地图并将其插入到新网格中。 它会很好地工作，而失真不是很大。
这是在 Python 中使用 scipy.interpolate.griddata 的非常简单的实现：
map_x, map_y = cv2.initUndistortRectifyMap(K, D, None, new_K, image_size, cv2.CV_32FC1)
points = np.stack([map_x.flatten(), map_y.flatten()], axis=1)
grid = np.mgrid[:map_x.shape[0], :map_y.shape[1]]
values = grid.reshape(2, -1).T[..., ::-1]
from scipy.interpolate import griddata
grid_y, grid_x = grid
map_back = griddata(points, values, (grid_x, grid_y), method='cubic').astype(map_undistort.dtype)
如果对地图使用CV_32FC2，可以简化点的构造：
map_undistort, _ = cv2.initUndistortRectifyMap(K, D, None, new_K, image_size, cv2.CV_32FC2)
points = map_undistort.reshape(-1, 2)
这是@wcochran 的回答的一个实现。我试图恢复由 lensfunpy 导致的镜头校正。
mod = lensfunpy.Modifier(lens, cam.crop_factor, width, height)
mod.initialize(focal_length, aperture, distance)
undist_coords = mod.apply_geometry_distortion()
## the lens correction part
# im_undistorted = cv2.remap(im, undist_coords, None, cv2.INTER_CUBIC)
# im_undistorted = cv2.remap(im, undist_coords, None, cv2.INTER_LANCZOS4)
# cv2.imwrite(undistorted_image_path, im_undistorted)
undist_coords_f = undist_coords.reshape((-1, 2))
tree = KDTree(undist_coords_f)
def calc_val(point_pos):
nearest_dist, nearest_ind = tree.query([point_pos], k=5)
if nearest_dist[0][0] == 0:
return undist_coords_f[nearest_ind[0][0]]
# starts inverse distance weighting
w = np.array([1.0 / pow(d, 2) for d in nearest_dist])
sw = np.sum(w)
# embed()
x_arr = np.floor(nearest_ind[0] / 1080)
y_arr = (nearest_ind[0] % 1080)
xx = np.sum(w * x_arr) / sw
yy = np.sum(w * y_arr) / sw
return (xx, yy)
un_correction_x = np.zeros((720, 1080))
un_correction_y = np.zeros((720, 1080))
## reverse the lens correction
for i in range(720):
print("row %d operating" % i)
for j in range(1080):
un_correction_x[i][j], un_correction_y[i][j] = calc_val((i, j))
# print((i, j), calc_val((j, i)))
dstMap1, dstMap2 = cv2.convertMaps(un_correction_x.astype(np.float32), un_correction_y.astype(np.float32), cv2.CV_32FC2)
im_un_undistorted = cv2.remap(im_undistorted, dstMap1, dstMap2, cv2.INTER_LANCZOS4)
这是一个重要的问题，我很惊讶没有在任何标准库中更好地解决这个问题（至少据我所知）。
我对公认的解决方案不满意，因为它没有使用转换的隐式平滑度。我可能会错过重要的案例，但我无法想象映射在任何有用的意义上都是可逆的，但在像素尺度上是不平滑的。
平滑意味着不需要计算最近邻：最近的点是那些在原始网格上已经很近的点。
我的解决方案使用的事实是，在原始映射中，正方形 [(i,j), (i+1, j), (i+1, j+1), (i, j+1 )] 映射到内部没有其他点的四边形 [(X[i,j], Y[i,j], X[i+1,j], Y[i+1,j], ...] .然后逆映射只需要在四边形内插值。为此我使用逆双线性插值，这将在顶点和任何其他仿射变换给出精确的结果。
除了 numpy 之外，该实现没有其他依赖项。逻辑是 运行 遍历所有四边形并逐步构建反向映射。我把代码复制到这里，希望有足够的评论让这个想法足够清晰。
关于不太明显的东西的一些评论：
- 反双线性函数通常 return 坐标仅在 [0,1] 范围内。我删除了裁剪操作，因此超出范围的值意味着坐标在四边形之外（这是解决多边形点问题的一种扭曲方式！）。为了避免丢失边缘上的点，我实际上允许 [0,1] 范围之外的点，这通常意味着两个相邻的四边形可能会拾取一个索引。在这些罕见的情况下，我只是让结果成为两个结果的平均值，相信超出范围的点是以合理的方式“外推”的。
- 一般来说，所有的四边形都有不同的形状，它们与规则网格的重叠可以从无到有变化很多点。该例程一次求解所有四边形（利用
bilinear_inverse的矢量化性质，但在每次迭代中仅选择坐标（到其边界框的偏移量）有效的四边形。
import numpy as np
def bilinear_inverse(p, vertices, numiter=4):
"""
Compute the inverse of the bilinear map from the unit square
[(0,0), (1,0), (1,1), (0,1)]
to the quadrilateral vertices = [p0, p1, p2, p4]
Parameters:
----------
p: array of shape (2, ...)
Points on which the inverse transforms are applied.
vertices: array of shape (4, 2, ...)
Coordinates of the vertices mapped to the unit square corners
numiter:
Number of Newton interations
Returns:
--------
s: array of shape (2, ...)
Mapped points.
This is a (more general) python implementation of the matlab implementation
suggested in
"""
p = np.asarray(p)
v = np.asarray(vertices)
sh = p.shape[1:]
if v.ndim == 2:
v = np.expand_dims(v, axis=tuple(range(2, 2 + len(sh))))
# Start in the center
s = .5 * np.ones((2,) + sh)
s0, s1 = s
for k in range(numiter):
# Residual
r = v[0] * (1 - s0) * (1 - s1) + v[1] * s0 * (1 - s1) + v[2] * s0 * s1 + v[3] * (1 - s0) * s1 - p
# Jacobian
J11 = -v[0, 0] * (1 - s1) + v[1, 0] * (1 - s1) + v[2, 0] * s1 - v[3, 0] * s1
J21 = -v[0, 1] * (1 - s1) + v[1, 1] * (1 - s1) + v[2, 1] * s1 - v[3, 1] * s1
J12 = -v[0, 0] * (1 - s0) - v[1, 0] * s0 + v[2, 0] * s0 + v[3, 0] * (1 - s0)
J22 = -v[0, 1] * (1 - s0) - v[1, 1] * s0 + v[2, 1] * s0 + v[3, 1] * (1 - s0)
inv_detJ = 1. / (J11 * J22 - J12 * J21)
s0 -= inv_detJ * (J22 * r[0] - J12 * r[1])
s1 -= inv_detJ * (-J21 * r[0] + J11 * r[1])
return s
def invert_map(xmap, ymap, diagnostics=False):
"""
Generate the inverse of deformation map defined by (xmap, ymap) using inverse bilinear interpolation.
"""
# Generate quadrilaterals from mapped grid points.
quads = np.array([[ymap[:-1, :-1], xmap[:-1, :-1]],
[ymap[1:, :-1], xmap[1:, :-1]],
[ymap[1:, 1:], xmap[1:, 1:]],
[ymap[:-1, 1:], xmap[:-1, 1:]]])
# Range of indices possibly within each quadrilateral
x0 = np.floor(quads[:, 1, ...].min(axis=0)).astype(int)
x1 = np.ceil(quads[:, 1, ...].max(axis=0)).astype(int)
y0 = np.floor(quads[:, 0, ...].min(axis=0)).astype(int)
y1 = np.ceil(quads[:, 0, ...].max(axis=0)).astype(int)
# Quad indices
i0, j0 = np.indices(x0.shape)
# Offset of destination map
x0_offset = x0.min()
y0_offset = y0.min()
# Index range in x and y (per quad)
xN = x1 - x0 + 1
yN = y1 - y0 + 1
# Shape of destination array
sh_dest = (1 + x1.max() - x0_offset, 1 + y1.max() - y0_offset)
# Coordinates of destination array
yy_dest, xx_dest = np.indices(sh_dest)
xmap1 = np.zeros(sh_dest)
ymap1 = np.zeros(sh_dest)
TN = np.zeros(sh_dest, dtype=int)
# Smallish number to avoid missing point lying on edges
epsilon = .01
# Loop through indices possibly within quads
for ix in range(xN.max()):
for iy in range(yN.max()):
# Work only with quads whose bounding box contain indices
valid = (xN > ix) * (yN > iy)
# Local points to check
p = np.array([y0[valid] + ix, x0[valid] + iy])
# Map the position of the point in the quad
s = bilinear_inverse(p, quads[:, :, valid])
# s out of unit square means p out of quad
# Keep some epsilon around to avoid missing edges
in_quad = np.all((s > -epsilon) * (s < (1 + epsilon)), axis=0)
# Add found indices
ii = p[0, in_quad] - y0_offset
jj = p[1, in_quad] - x0_offset
ymap1[ii, jj] += i0[valid][in_quad] + s[0][in_quad]
xmap1[ii, jj] += j0[valid][in_quad] + s[1][in_quad]
# Increment count
TN[ii, jj] += 1
ymap1 /= TN + (TN == 0)
xmap1 /= TN + (TN == 0)
if diagnostics:
diag = {'x_offset': x0_offset,
'y_offset': y0_offset,
'mask': TN > 0}
return xmap1, ymap1, diag
else:
return xmap1, ymap1
这是一个测试例子
import cv2 as cv
from scipy import ndimage as ndi
# Simulate deformation field
N = 500
sh = (N, N)
t = np.random.normal(size=sh)
dx = ndi.gaussian_filter(t, 40, order=(0,1))
dy = ndi.gaussian_filter(t, 40, order=(1,0))
dx *= 30/dx.max()
dy *= 30/dy.max()
# Test image
img = np.zeros(sh)
img[::10, :] = 1
img[:, ::10] = 1
img = ndi.gaussian_filter(img, 0.5)
# Apply forward mapping
yy, xx = np.indices(sh)
xmap = (xx-dx).astype(np.float32)
ymap = (yy-dy).astype(np.float32)
warped = cv.remap(img, xmap, ymap ,cv.INTER_LINEAR)
plt.imshow(warped, cmap='gray')
# Now invert the mapping
xmap1, ymap1 = invert_map(xmap, ymap)
unwarped = cv.remap(warped, xmap1.astype(np.float32), ymap1.astype(np.float32) ,cv.INTER_LINEAR)
plt.imshow(unwarped, cmap='gray')
迭代求解
上面的许多解决方案对我都不起作用，当贴图不可逆时失败，或者速度不是很快。
我提出了另一种 6 行迭代解决方案。
def invert_map(F):
I = np.zeros_like(F)
I[:,:,1], I[:,:,0] = np.indices(sh)
P = np.copy(I)
for i in range(10):
P += I - cv.remap(F, P, None, interpolation=cv.INTER_LINEAR)
return P
效果如何？ 对于我为航空摄影反转地形校正图的用例，这种方法在 10 个步骤中轻松收敛到 1/10 像素。它也非常快，因为所有繁重的计算都隐藏在 OpenCV
中它是如何工作的？
该方法使用的思想是，如果 (x', y') = F(x, y) 是一个映射，那么只要 F 的梯度很小，就可以用 (x, y) = -F(x', y') 来近似逆。
我们可以继续完善我们的映射，上面得到了我们的第一个预测（我是一个“身份映射”）：
G_1 = I - F
我们的第二个预测可以改编自：
G_2 = G_1 + I - F(G_1)
等等：
G_n+1 = G_n + I - F(G_n)
证明 G_n 收敛到逆 F^-1 很难，但我们可以轻松证明的是，如果 G 已经收敛，它将保持收敛。
假设G_n = F^-1，那么我们可以代入：
G_n+1 = G_n + I - F(G_n)
然后得到：
G_n+1 = F^-1 + I - F(F^-1)
G_n+1 = F^-1 + I - I
G_n+1 = F^-1
Q.E.D.
测试脚本
import cv2 as cv
from scipy import ndimage as ndi
import numpy as np
from matplotlib import pyplot as plt
# Simulate deformation field
N = 500
sh = (N, N)
t = np.random.normal(size=sh)
dx = ndi.gaussian_filter(t, 40, order=(0,1))
dy = ndi.gaussian_filter(t, 40, order=(1,0))
dx *= 10/dx.max()
dy *= 10/dy.max()
# Test image
img = np.zeros(sh)
img[::10, :] = 1
img[:, ::10] = 1
img = ndi.gaussian_filter(img, 0.5)
# Apply forward mapping
yy, xx = np.indices(sh)
xmap = (xx-dx).astype(np.float32)
ymap = (yy-dy).astype(np.float32)
warped = cv.remap(img, xmap, ymap ,cv.INTER_LINEAR)
plt.imshow(warped, cmap='gray')
def invert_map(F: np.ndarray):
I = np.zeros_like(F)
I[:,:,1], I[:,:,0] = np.indices(sh)
P = np.copy(I)
for i in range(10):
P += I - cv.remap(F, P, None, interpolation=cv.INTER_LINEAR)
return P
# F: The function to invert
F = np.zeros((sh[0], sh[1], 2), dtype=np.float32)
F[:,:,0], F[:,:,1] = (xmap, ymap)
# Test the prediction
unwarped = cv.remap(warped, invert_map(F), None, cv.INTER_LINEAR)
plt.imshow(unwarped, cmap='gray')
一种方法是获取原始地图，遍历其条目并获取 x 和 y 值的下限和上限。这给出了 (x,y), (xf,yf), (xc,yf), (xf,yc), (x c,yc) 在原始源图像的坐标中。然后，您可以填充一个结构，其中每一个都作为包含像素值和权重的索引，并对这些数据使用您喜欢的不规则网格插值。
这很容易通过反距离插值来实现，因为结构可以是图像数组累加，而权重是标量。 F 是原始源，G 是变形后的图像，F' 是恢复后的图像。地图是M.
将 F' 初始化为 0。创建一个与 F' 大小相同的浮点数的 0 初始化权重数组 W。
遍历 M。对于 M 中的每个整数，找到 4 个整数对及其与 (x,y) 的距离。从G中取出对应的像素值，按其距离的倒数加权，累加到F' like
F'(xf|c,yf|c)+=G(i,j)/sqrt((x-xf|c)^2+(y-yf|c)^2)
然后把那个权重累加到
W(xf|c,yf|c)+=1./sqrt((x-xf|c)^2+(y-yf|c)^2).
完成后，通过迭代对 F' 进行归一化，并将每个像素除以其在 W 中的对应条目（如果它不为零）。
在这一点上，图像通常接近完成，但是由于高下采样率，F' 中的一些像素可能没有被填充。所以你在 W 中来回传递几次以找到 0 权重条目，并从它们的非空邻居中插入这些像素。这部分也可以使用 KNN 搜索和插值来完成，因为它们通常不多。
与 KNN 方法相比，它易于实现并且缩放性更好（尽管我认为这对于小图像非常有用）。缺点是反距离不是最好的插值方案，但如果映射不是太块并且原始没有被大量下采样，它似乎工作得很好。当然，如果下采样率很高，你就不得不推断出很多丢失的信息，所以它本质上会给出粗略的结果。
如果你想尽可能多地从地图反演中挤出来，你可以尝试求解由原始插值方案定义的（可能欠定的）方程组；并非不可能，但具有挑战性。
KNNRegressor 具有反转网格映射的所有必要组件！
给你：
from sklearn.neighbors import KNeighborsRegressor
def get_inverse_maps(map1, map2):
regressor = KNeighborsRegressor(3)
X = np.concatenate((map2[..., None], map1[..., None]), axis=-1).reshape(-1, 2)
y = np.indices(map1.shape).transpose((1, 2, 0)).reshape(-1, 2)
regressor.fit(X, y)
map_inv = regressor.predict(y).reshape(map1.shape + (2,)).astype(np.float32)
map_inv2, map_inv1 = map_inv[..., 0], map_inv[..., 1]
return map_inv1, map_inv2