Numpy实现透视变换

2019-12-30

其实就是用numpy.linalg.solve解方程啦

透视变换就是三维空间上的一个点乘上变换矩阵转移到别的坐标（因为图片本身是二维的，所以新旧坐标的z值都为1），即：

$$\begin{bmatrix} XW \\ YW \\ W \end{bmatrix} = \begin{bmatrix} a & b & c \\ d & e & f \\ g & h & 1 \end{bmatrix} \begin{bmatrix} x \\ y \\ 1 \end{bmatrix}$$

解得新坐标：

$$X = \frac {ax+by+c} W = \frac {ax+by+c} {gx+hy+1}$$

$$Y = \frac {dx+ey+f} W = \frac {dx+ey+f} {gx+hy+1}$$

但其实a-h是未知的，所以得先通过几对给定的新旧点来解出它们，才能得到新旧图片上所有点的映射。

解a-h的步骤如下：首先把分母与左边的X/Y相乘，再整理等式把单独的X/Y拿到等号右边可以得到两个方程：

$$ax+by+\phantom 1 c-0d+0e+0f-Xxg-Xyh=X$$

$$0a+0b+0c+xd+ye+\phantom 1 f-Yxg-Yyh=Y$$

转换成矩阵形式：

$$\begin{bmatrix} x & y & 1 & 0 & 0 & -1 & -Xx & -Xy \\ 0 & 0 & 0 & x & y & 1 & -Yx & -Yy \\ & & & & & \vdots \end{bmatrix} \begin{bmatrix} a \\ b \\ c \\ d \\ e \\ f \\ g \\ h \end{bmatrix} = \begin{bmatrix} X \\ Y \\ \vdots \end{bmatrix}$$

8个未知数需要8个方程来解，每对新旧点提供两个方程，所以需要4对点来填满矩阵，求解a-h相当于解形如Ax=b的方程。numpy.linalg.solve就是用来解Ax=b的函数，用A与b当参数调用就可以得到x（a-h），这样我们就得到了变换矩阵。

构造A、b与解变换矩阵的代码如下：

A_list = []
b_list = []
for xy, XY in zip(old_points, new_points):
    A_list.append([xy[0], xy[1], 1, 0, 0, 0, -XY[0] * xy[0], -XY[0] * xy[1]])
    A_list.append([0, 0, 0, xy[0], xy[1], 1, -XY[1] * xy[0], -XY[1] * xy[1]])
    b_list.extend(XY)

a, b, c, d, e, f, g, h = np.linalg.solve(A_list, b_list)
H = np.array([[a, b, c], [d, e, f], [g, h, 1]])

拿到变换矩阵遍历新图的每一处(X, Y)坐标，再通过numpy.linalg.solve就可以算出变换前的(x, y)坐标：

$$\begin{bmatrix} a & b & c \\ d & e & f \\ g & h & 1 \end{bmatrix} \begin{bmatrix} x/W \\ y/W \\ 1/W \end{bmatrix} = \begin{bmatrix} X \\ Y \\ 1 \end{bmatrix}$$

xy1 = np.linalg.solve(H, [X, Y, 1])
x, y, _ = (1 / xy1[2] * xy1).astype(int)

然后把原图(x, y)处的灰度值拿来放到新图的(X, Y)处即可（是的，这里没考虑插值啥的）：

if 0 <= x < im.shape[0] and 0 <= y < im.shape[1]:
    newim[X, Y] = im[x, y]

效果图：

完整代码在这里

参考了这篇文章