反转实值索引网格

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和和的图像A的像素值y.

我的问题是:给定索引网格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？

FWIW,图像和索引图X和Y是相同的形状.然而,索引图X和Y没有先验结构.例如,它们不一定是仿射或刚性变换.它们甚至可能是不可逆的,例如,如果X, Y将多个像素映射A到B中的相同精确像素坐标.我正在寻找一种方法的想法,如果存在,将找到合理的逆映射.

解决方案不必是基于OpenCV的,因为我没有使用OpenCV,而是另一个具有remap()实现的库.虽然欢迎任何建议,但我特别热衷于"数学上正确"的事情,即如果我的地图M完全可逆,那么该方法应该在机器精度的一小部分范围内找到完美的反转.

Answer 1

迭代解决方案

上述许多解决方案对我来说不起作用，当地图不可逆时失败，或者速度不是很快。

我提出了另一种 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')

Answer 2

好吧，我只需要自己解决此重新映射反转问题，然后我将概述解决方案。

给定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])

首先，我为2D点集建立一个KD树，{(X[i,j],Y[i,j]}这样我就可以有效地找到N给定点的最近邻居，而(x,y).我使用欧几里得距离作为距离度量。我在GitHub上为KD-Trees找到了一个很棒的C ++头文件库。

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

• 如果距离d_k == 0对于一些k则Xinv[x,y] = i_k和Yinv[x,y] = j_k，否则......

• 使用反距离权重（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()用于消除第一张图像中的镜头畸变。第三幅图像是逆过程的结果。

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://stackoverflow.com/a/18332009/1560876
    """

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