[VP] 恢复图像镜面扭曲 - 多视觉几何

---

一、算法原理

感谢老哥提供的图！

图像中的四条线在现实世界中应该是平行的，或者可以说现实世界中的平行线经过扭曲变换 H 变成了上图那样。

那么我们拿到一张扭曲后的图，首先应该是寻找四个点，这四个点互相构成的四条线在现实世界应该两两平行，如上图所示。

第二步，图像中的四个点确定的四条线，会有两个交点 (我们称为“消失点”) ，由这两个消失点，我们就可以确定一条 “消失线” 。

由两个点求线的公式：$l=P_1\times P_2$

由两条线求交点公式：$P=l_1\times l_2$

第三步，单应矩阵 $H=H_A{H}_P$，其中 $H_P$为：

$$\left[\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ l_1& l_2& l_3\end{array}\right]$$

其中的$l_1,l_2,l_3$ 就是第二步计算出的消失线， $H_A$可以为任何的仿射变换矩阵。

第四步，知道了 $H_P$就可以用 Skimage.transform.warp 函数来恢复图像，或者DLT算法恢复

恢复后的图如下：

二、代码实现

1、确定四个点

import cv2
import sympy
import numpy as np
from PIL import Image
from skimage.transform import *
import matplotlib.pyplot as plt
from scipy.interpolate import  griddata
from mpl_toolkits.mplot3d import Axes3D

#四个点坐标
Pa = np.array([[270], [300], [1]])
Pb = np.array([[454], [230], [1]])
Pc = np.array([[637], [300], [1]])
Pd = np.array([[455], [400], [1]])

#绘制一下这四个点看看
img = cv2.imread('Ex-2.PNG')
plt.plot(270, 300, 'o', color='r')
plt.plot(454, 230, 'o', color='r')
plt.plot(637, 300, 'o', color='r')
plt.plot(455, 400, 'o', color='r')
plt.imshow(img)

2、由这四个点确定四条线

# Lines
l_ad = np.cross(Pa.T, Pd.T)
l_ad = l_ad / l_ad[0][2]
print(l_ad)

l_bc = np.cross(Pb.T, Pc.T)
l_bc = l_bc / l_bc[0][2]
print(l_bc)

l_ab = np.cross(Pa.T, Pb.T)
l_ab = l_ab / l_ab[0][2]
print(l_ab)

l_cd = np.cross(Pc.T, Pd.T)
l_cd = l_cd / l_cd[0][2]
print(l_cd)

[[ 0.00350877 -0.00649123 1. ]]
[[ 0.00678952 -0.01774976 1. ]]
[[-9.44669366e-04 -2.48313090e-03 1.00000000e+00]]
[[-8.45308538e-04 -1.53846154e-03 1.00000000e+00]]

3、由四条线确定两个消失点

P1 = np.cross(l_ad, l_bc)
P1 = P1 / P1[0][2]
print("P1: ", P1)

P2 = np.cross(l_ab, l_cd)
P2 = P2 / P2[0][2]
print("P2: ", P2)

P1: [[-618.34579439 -180.18691589 1. ]]
P2: [[ 1.46307420e+03 -1.53886926e+02 1.00000000e+00]]

4、由这两个消失点确定消失线

Line = np.cross(P1, P2)
Line = Line / Line[0][2]
Line

array([[-7.33035052e-05, 5.80134750e-03, 1.00000000e+00]])

5、由消失线确定$H_P$

Hp = np.matrix([[1, 0, 0], [0, 1, 0], [Line[0][0], Line[0][1], Line[0][2]]])
Hp

matrix([[ 1.00000000e+00, 0.00000000e+00, 0.00000000e+00],
[ 0.00000000e+00, 1.00000000e+00, 0.00000000e+00],
[-7.33035052e-05, 5.80134750e-03, 1.00000000e+00]])

6、由warp函数和$H_P$复原图像

plt.imshow(warp(img, np.linalg.inv(Hp)))

7、用DLT复原图像

def image_rebound(mm,nn,hh):
    W = np.array([[1, nn, nn, 1 ],[1, 1, mm, mm],[ 1, 1, 1, 1]])
    ws = np.dot(hh,W)
    ### scaling
    xx = np.vstack((ws[2,:],ws[2,:],ws[2,:]))
    wsX =  np.round(ws/xx)
    bounds = [np.min(wsX[1,:]), np.max(wsX[1,:]),np.min(wsX[0,:]), np.max(wsX[0,:])]
    return bounds

def make_transform(imm,hh):   
    mm,nn = imm.shape[0],imm.shape[0]
    bounds = image_rebound(mm,nn,hh)
    nrows = bounds[1] - bounds[0]
    ncols = bounds[3] - bounds[2]
    s = max(nn,mm)/max(nrows,ncols)
    scale = np.array([[s, 0, 0],[0, s, 0], [0, 0, 1]])
    trasf = scale@hh
    trasf_prec =  np.linalg.inv(trasf)
    bounds = image_rebound(mm,nn,trasf)
    nrows = (bounds[1] - bounds[0]).astype(np.int)
    ncols = (bounds[3] - bounds[2]).astype(np.int)
    return bounds, nrows, ncols, trasf, trasf_prec

def get_new_image(nrows,ncols,imm,bounds,trasf_prec,nsamples):
    xx  = np.linspace(1, ncols, ncols)
    yy  = np.linspace(1, nrows, nrows)
    [xi,yi] = np.meshgrid(xx,yy) # 生成网格点坐标矩阵。
    a0 = np.reshape(xi, -1,order ='F')+bounds[2]
    a1 = np.reshape(yi,-1, order ='F')+bounds[0]
    a2 = np.ones((ncols*nrows))
    uv = np.vstack((a0.T,a1.T,a2.T)) 
    new_trasf = np.dot(trasf_prec,uv)
    val_normalization = np.vstack((new_trasf[2,:],new_trasf[2,:],new_trasf[2,:]))
   
    ### The new transformation
    newT = new_trasf/val_normalization
    
    ### 
    xi = np.reshape(newT[0,:],(nrows,ncols),order ='F') 
    yi = np.reshape(newT[1,:],(nrows,ncols),order ='F')
    cols = imm.shape[1]
    rows = imm.shape[0]
    xxq  = np.linspace(1, rows, rows).astype(np.int)
    yyq  = np.linspace(1, cols, cols).astype(np.int)
    [x,y] = np.meshgrid(yyq,xxq) 
    x = (x - 1).astype(np.int) #Offset x and y relative to region origin.
    y = (y - 1).astype(np.int) 
        
    ix = np.random.randint(im.shape[1], size=nsamples)
    iy = np.random.randint(im.shape[0], size=nsamples)
    samples = im[iy,ix]
       
    int_im = griddata((iy,ix), samples, (yi,xi))
    
    #Plotting
    fig = plt.figure(figsize=(8, 8))
    columns = 2
    rows = 1
    fig.add_subplot(rows, columns, 1)
    plt.imshow(im)
    fig.add_subplot(rows, columns, 2)
    
    plt.imshow(int_im.astype(np.uint8))
    plt.show()

im = np.array(Image.open('Ex-2.PNG'))
bounds, nrows, ncols,  trasf, trasf_prec = make_transform(im,Hp)

nn,mm  = im.shape[0],im.shape[0]
if max(nn,mm)>1000:
    kk =6
else: kk =5
nsamples = 10**kk         
get_new_image(nrows,ncols,im,bounds,trasf_prec,nsamples)

8、确定 $H^{-T}{l}^T$在无穷远处的一条线上

即确定$H^{-T}{l}^T=\left[\begin{array}{c}0\\ 0\\ 1\end{array}\right]$，此处使用的是sympy的符号系统，以证明任何数字均满足此公式

from sympy import *

h1 = sympy.Symbol('h1')
h2 = sympy.Symbol('h2')
h3 = sympy.Symbol('h3')
h4 = sympy.Symbol('h4')
h5 = sympy.Symbol('h5')
h6 = sympy.Symbol('h6')

l1 = sympy.Symbol('l1')
l2 = sympy.Symbol('l2')
l3 = sympy.Symbol('l3')

HaMatrix = sympy.Matrix([[h1, h2, h3], [h4, h5, h6], [0, 0, 1]])
HpMatrix = sympy.Matrix([[1, 0, 0], [0, 1, 0], [l1, l2, l3]])
HpMatrix

H = HaMatrix*HpMatrix
H

l = Matrix([[l1], [l2], [l3]])
l

simplify(H.T.inv()*l)