反转实值索引网格

Question

OpenCV的remap()使用实值索引网格对图像进行双线性插值采样，返回采样网格作为新图像。

具体而言，请设：

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 的图像A的像素值。

我的问题是：给定索引网格 X，Y，如何生成“反向网格” 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？

顺便说一下，图像和索引映射X和Y具有相同的形状。但是，索引映射X和Y没有先验结构。例如，它们不一定是仿射或刚性变换。它们甚至可能无法反转，例如，如果X，Y将A中的多个像素映射到B中的完全相同的像素坐标，则无法反转。我正在寻找一种方法，如果存在合理的反向映射，则会找到它。

解决方案不必基于OpenCV，因为我没有使用OpenCV，而是另一个具有remap()实现的库。虽然欢迎任何建议，但我特别希望得到一些“数学上正确”的东西，即如果我的映射M是完美可逆的，则该方法应在某些机器精度小的范围内找到完美的反向映射。

- SuperElectric

Answer 1

我刚刚不得不自己解决了这个重新映射反转问题，现在我将概述我的解决方案。

remap()函数的输入参数为X和Y，其功能如下：

B[i, j] = A(X[i, j], Y[i, j])   

我计算了Xinv和Yinv，可以用于remap()函数中，以反转该过程：

A[x, y] = B(Xinv[x,y],Yinv[x,y])

首先，我为二维点集{(X[i,j],Y[i,j])}构建了一个KD-Tree，以便可以高效地找到给定点(x,y)的N个最近邻居。我的距离度量使用欧几里得距离。我在GitHub上找到了一个很好的C++头文件库用于KD-Tree。

然后，我循环遍历A的网格中的所有(x,y)值，并在我的点集中找到N = 5个最近邻居{(X[i_k,j_k],Y[i_k,j_k]) | k = 0 .. N-1}。

• 如果对于某些 k，距离d_k == 0，那么Xinv[x,y] = i_k并且Yinv[x,y] = j_k，否则...

• 使用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()函数来消除第一张图片的镜头畸变。第三张图片是反转这个过程的结果。

- wcochran

---

你好，能分享一下你的代码吗？我尝试实现了一下，但结果图像没有反转镜头校正。 - SCaffrey

@SCaffrey 很高兴它起作用了。我本来可能会发布我的代码链接，但它是我无法分享的一个更大项目的一部分。我可以提取有趣的部分并发布它，但听起来你已经得到了它。 - wcochran

您可以在此处找到这个聪明想法的Python实现。 - basfest

我可能漏掉了什么，但这个解决方案似乎既武断又浪费。武断：为什么是5个邻居？为什么是反距离加权？浪费：为什么不使用网格结构来获取邻居？唯一有意义的应用是那些映射既平滑又双射的情况，在这种情况下不需要昂贵的邻居搜索。 - pthibault

@pthibault 这些点 {(X[i,j],Y[i,j])} 不是 在一个规则的网格上。在径向镜头畸变的情况下，它们位于被拉伸的像牛肉干一样的网格上（例如，在角落处看到第二张图像）。换句话说，这些点代表了非线性函数的样本 -- 当你有重度鱼眼畸变并添加切向畸变时，你可以忘记使用网格来查找邻居。N = 5个邻居是一个很好的金发姑娘值 -- 不太少也不太多。IDW是插值非规则网格上采样的2D表面的绝佳方法。 - wcochran

Answer 2

迭代解法

许多以上的解决方案对我来说并不适用，当地图不可逆时会失败，或者速度不够快。

我提供一种替代方法，使用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) 是一个映射，那么可以通过 (x, y) = -F(x', y') 来近似求解其反函数，只要 F 的梯度很小。

我们可以继续改进我们的映射，上面的式子得到了我们的第一次预测（I 是“恒等映射”）：

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

- Hannesh

---

我喜欢这个，因为它简单易懂且有坚实的数学基础。 - Christoph Rackwitz

重新映射有限精度，因为标准插值模式在内部使用固定点数学。据我所知，它仅限于5个小数位。 - Christoph Rackwitz

如果解决方案在收敛时出现问题，或者在图像边缘附近存在大的偏差，这可能是由于“P += I - cv.remap（...）”步骤中的振荡引起的，那么只需计算“P +=（I - cv.remap（...））* 0.5”即可。通过这种调整，迭代被阻尼，并且在仅5-10次迭代内收敛（下降到remap的固定点数值限制，这给您提供了5个小数位）。 - Christoph Rackwitz

很棒的答案。非常简单。在实践中，前向映射可能需要丢失源数据。您能修改您的答案以应对这种情况吗？测试代码的这个调整说明了这一点：将xmap = (xx-dx).astype(np.float32)替换为xmap = (0.1*N+xx*0.8-dx).astype(np.float32)。我们现在看到图像左侧出现了损坏。有趣的是右侧没有！请注意，在cv.remap()行中添加borderMode=cv.BORDER_CONSTANT, borderValue=0.5有助于说明预期的无效像素和任何不正确像素之间的差异。 - ianinini

很棒的回答。非常简单易懂。在实践中，前向映射可能需要丢失源数据。你能修改你的回答来应对这个情况吗？对测试代码进行的这个调整阐明了这一点：用xmap = (0.1*N+xx*0.8-dx).astype(np.float32)替换xmap = (xx-dx).astype(np.float32)。现在我们看到图像左侧出现了损坏。有趣的是右侧没有出现！注意，向cv.remap()行添加borderMode=cv.BORDER_CONSTANT, borderValue=0.5有助于说明预期的无效像素和任何不正确像素之间的差异。 - undefined

Answer 3

这是一个重要的问题，而且我很惊讶任何标准库都没有更好地解决它（至少据我所知）。

我对已接受的解决方案并不满意，因为它没有利用变换的隐式平滑性。我可能会错过一些重要情况，但我无法想象既具有可逆性又在像素尺度上非平滑的映射。

平滑性意味着无需计算最近邻：最近的点是原始网格上已经靠近的那些点。

我的解决方案利用了这样一个事实，即在原始映射中，一个正方形[(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 之外没有任何其他依赖项。逻辑是遍历所有四边形，并逐步构建反向映射。我在这里复制代码，希望有足够的注释使想法变得足够清晰。

关于不太明显的内容有一些评论：

• 反双线性函数通常只返回范围为[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 https://dev59.com/5HRA5IYBdhLWcg3w1BmW#18332009
    """

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

- pthibault

---

sh_dest = ... 中存在错误，顺序应为（高度，宽度）-- 代码执行需要一段时间 -- 并且似乎生成的地图包含杂乱的0索引（当需要拉伸输入时）... 我在最后一张图片（结果）中看到了这方面的提示。 - Christoph Rackwitz

Answer 4

您可以在已知点处反转地图并将其插值到新网格中。 只要扭曲不是非常巨大，它就能正常工作。

这是使用scipy.interpolate.griddata的Python非常简单的实现:

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)

- FelisRattus

Answer 5

如果您的地图是从单应性矩阵 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()的说明如下：

该函数实际上为remap()使用的反向映射算法构建映射。也就是说，对于目标图像中的每个像素(u, v)，该函数会计算出源图像中相应的坐标。

在你只有地图的情况下，你需要自己完成这个操作。然而，新地图坐标的插值并不简单，因为一个像素的支持区域可能非常大。

下面是一个简单的Python解决方案，它通过点对点映射来反转地图。这可能会导致一些坐标未分配，而其他坐标会被更新多次。因此地图上可能会有空洞。

下面是一个小的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()

- Tobias

---

谢谢，但正如问题所述，地图没有先验结构。它不一定是一个单应性变换。 - SuperElectric

您可以使用中值模糊函数cv2::medianBlur()填充从点到点映射中的空洞。 - Tobias

当然可以，但这会引入一个任意的参数（核大小），并且对于那些没有很多空洞的图像来说，可能会对准确性造成更多的伤害，因为会模糊不需要模糊的区域。在这里，准确性指的是M[M_inv] - I的Frobenius范数，其中M是我们试图反转的索引映射，a[b]表示将a重新映射为b，M_inv是M的逆映射，而I是恒等映射，即对于所有i，j，I[i，j，：] == [i，j]的映射。（在我的用例中，M_inv必须通过这些措施“在数学上正确”，而不仅仅是产生合理的图像）。 - SuperElectric

我认为解决方案将涉及点对点和单应性映射的混合。设M是从图像A到B的映射，N是从图像B到A的映射（我们要找的东西）。 M中的每个2x2像素块都定义了一个单应性H，因为它们的像素值（A的坐标）构成了一个四边形，该四边形被映射到由它们的像素坐标（B的坐标）形成的正方形。为了定义N，我们需要遍历A的像素坐标[i，j]，并为每个像素坐标确定它所在的M中的四边形，并将[i，j]乘以该四边形的H^-1。 - SuperElectric

这个解决方案对我有效，但是我花了一些时间才意识到结果是正确的，因为图像的某些部分被边缘裁剪了。 - Georgy

对我来说，即使我将两个参数作为cv2.fisheye.initUndistortRectifyMap的结果得到，invert_maps仍会崩溃并出现assert错误。 - Stepan Yakovenko

Answer 6

这是@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)

- SCaffrey

Answer 7

一个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

- Roulbacha

Answer 8

使用OpenCV没有任何标准方法来完成此任务。

如果您正在寻找完整的即用型解决方案，我不确定能否提供帮助，但我至少可以描述一种我几年前用于执行此任务的方法。

首先，您应该创建与源图像具有相同尺寸的重新映射映射。我创建了具有较大尺寸的映射以进行更简单的插值，并在最后一步将它们裁剪到适当的大小。然后，您应该使用先前的重新映射映射中存在的值填充它们（不太困难：只需迭代它们，如果映射坐标x和y位于图像的限制范围内，则将它们的行和列作为新的y和x，并将其放置在旧的x和y列和行的新地图中）。这是一个相当简单的解决方案，但它可以得到相当好的结果。对于完美的解决方案，您应该使用您的插值方法和邻近像素将旧的x和y插值为整数值。

之后，您应该手动实际重新映射像素颜色，或者完全填充您的重新映射映射与像素坐标并使用OpenCV版本。

你将面临相当具有挑战性的任务：你需要在空白区域内插值像素。换句话说，你应该根据最近的非零像素坐标距离混合颜色（如果重新映射颜色）或坐标（如果进行完整的地图计算）。实际上，对于线性插值来说也并不是很困难，你甚至可以查看OpenCV github page中remap()的实现。对于NN插值来说，它会更简单 - 只需取最近邻的颜色/坐标。

最后一个任务是超出重新映射像素区域边界的区域外推。同样可以使用OpenCV中的算法作为参考。

- avtomaton

Answer 9

一种方法是获取原始地图，遍历其条目并取x和y值的floor和ceil。这会在原始源图像的坐标系中给出(x,y)周围的四个最近整数：(x_f,y_f)、(x_c,y_f)、(x_f,y_c)和(x_c,y_c)。然后，您可以使用每个整数作为索引填充一个结构，其中包含像素值和权重，并使用所述数据使用首选的不规则网格插值。

使用反距离插值很容易实现，因为该结构可以是图像数组累积，而权重是标量。F是原始源，G是扭曲图像，F'是恢复的图像。映射是M。

将F'初始化为0。创建一个与F'大小相同的浮点数的0初始化权重数组W。

遍历M。对于M中的每个元素，找到4个整数对及其与(x,y)的距离。从G中取出相应的像素值，通过其倒数距离加权，并将其累积到F'中，如下所示：

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'。

此时，图像通常已经接近完成，但是在高下采样比率下，F'中的某些像素可能没有填充。因此，您需要通过W进行几次前后传递，以查找0权重条目，并从其非空邻居插值这些像素。由于通常没有太多这样的像素，因此也可以使用KNN搜索和插值来完成此部分。

它比KNN方法更易于实现和扩展（虽然我认为对于小图像来说这很棒）。缺点是逆距离不是最好的插值方案，但如果映射不太杂乱且原始图像没有被大量降采样，则似乎效果还不错。当然，如果降采样比例很高，您将不得不推断出许多丢失的信息，因此结果本质上会变得粗糙。

如果您想尽可能地从地图反演中挤取出更多内容，可以尝试解决由原始插值方案定义的（潜在欠定的）方程组；这并非不可能，但具有挑战性。

- Zgyz McFitzgyzson

Answer 10

楼上的，我想我找到了答案。我还没有实现它，如果有人提供了一个更简单的解决方案(或者发现这个解决方案有问题)，我会选择他们的答案。

问题陈述

设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.

...其中方括号表示使用整数索引进行数组查找，圆括号表示使用浮点索引进行双线性插值查找。我们可以使用更经济的符号重新表述上述内容：

B = A(M)

我们希望找到一个反向映射 N，它可以尽可能好地将 B 映射回 A：

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具有相同维度的恒等映射。

存在两个问题：

1. M不是同构，因此上述方法将在N[i，j，：]处留下“空洞”，对于任何不在M值中的[i，j]。
2. M的值是浮点坐标[i，j]，而不是整数坐标。我们不能简单地为浮点值i，j分配双线性插值数量N（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，:)等于M中所有像素映射[i，j] -> [k，l]的所需值[k，l]。

解决方案

将空N构造为浮点数的3D张量：

N = zeros(size=(A.shape[0], A.shape[1], 2))

对于 A 坐标空间中的每个坐标 [i, j]，执行以下操作：

1. 在 M 中找到包含 [i, j] 的 A 坐标的 2x2 网格。计算将这些 A 坐标映射到其相应的 B 坐标（由 2x2 网格的像素索引给出）的单应矩阵 H。
2. 设置 N[i, j, :] = matmul(H, [i, j])

这里可能会产生昂贵的步骤，即在第 1 步中搜索包围 [i, j] 的 A 坐标的 2x2 网格。暴力搜索会使整个算法的时间复杂度为 O(n*m)，其中 n 是 A 中的像素数，m 是 B 中的像素数。

为了将其降低到 O(n)，可以在每个 A 坐标四边形内运行扫描线算法，以识别它包含的所有整数值坐标 [i, j]。这可以预先计算为哈希映射，将整数值 A 坐标 [i, j] 映射到其包围四边形的 B 坐标 [k, l] 的左上角。

- SuperElectric