Python计算机视觉编程第五章

最新推荐文章于 2024-08-21 10:28:58 发布

asdaasddsa

最新推荐文章于 2024-08-21 10:28:58 发布

阅读量669

点赞数 18

文章标签：计算机视觉人工智能

本文链接：https://blog.csdn.net/zxsdss/article/details/140946011

版权

一、外极几何

多视图几何是利用在不同视点所拍摄图像间的关系，来研究照相机之间或者特征之间关系的一门科学。图像的特征通常是兴趣点。多视图几何中最重要的内容是双视图几何。

如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。通过外极几何来描述这些几何关系。

当分析双视图几何关系时，可以将照相机间的相对位置关系用单应性矩阵加以简化。这里的单应性矩阵通常只做刚体变换，即只通过单应矩阵变换了坐标系。同一个图像点经过不同的投影矩阵产生的不同投影点必须满足：

$\mathbf{x}_2^T\boldsymbol{F}\mathbf{x}_1=0$

其为外极约束条件。矩阵 F 为基础矩阵。基础矩阵可以由两照相机的参数矩阵表示。其表明可以借助 F 来恢复出照相机参数，而 F 可以从对应的投影图像点计算出来。在知道照相机内参数的情况下，重建是基于度量的。所谓度量重建，即能够在三维重建中正确地表示距离和角度。

外极几何的示意图如下，

1.一个简单的数据集

首先在网上下载指定的资源之后，把这些资源加载下来：

img1 = np.array(Image.open(r'D:\test\merton\001.jpg'))
img2 = np.array(Image.open(r'D:\test\merton\002.jpg'))
points2D = [np.loadtxt(r'D:\test\merton\2D\00'+str(i+1)+'.corners').T for i in range(3)]
points3D = np.loadtxt(r'D:\test\merton/p3d').T
corr = np.genfromtxt(r'D:\test\merton\2D\nview-corners',dtype='int',missing='*')
P = [camera.Camera(np.loadtxt(r'D:\test\merton\2D\00'+str(i+1)+'.P')) for i in range(3)]

接着对其进行可视化操作：

X = np.vstack( (points3D,np.ones(points3D.shape[1])) )
x = P[0].project(X)
# 在视图 1 中绘制点
plt.figure()
plt.subplot(121)
plt.imshow(img1)
plt.plot(points2D[0][0],points2D[0][1],'*')
plt.axis('off')
plt.subplot(122)
plt.imshow(img1)
plt.plot(x[0],x[1],'r.')
plt.axis('off')
plt.show()

结果为：

左图上的点为图像点，右图上的点为投影的三维点。右图比左图多一些点，这些多出的点是从视图 2 和视图 3 重建出来的。

2.用Matplotlib绘制三维数据

为了可视化三维重建结果，我们需要绘制出三维图像。可以使用 Matplotlib 中的 mplot3d 工具包来进行绘制。通过在 axes 对象中加上 projection="3d" 关键字实现三维绘图，即：

from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import axes3d

fig = plt.figure()
ax = fig.gca(projection="3d")
X,Y,Z = axes3d.get_test_data(0.25)
ax.plot(X.flatten(),Y.flatten(),Z.flatten(),'o')
plt.show()

get_test_data() 函数在 x, y 空间按照设定的空间间隔参数来产生均匀的采样点。之后可以将产生的三列数据点输入到plot()函数，下面为一个对前面的例子的数据三维点的效果：

X = np.vstack( (points3D,np.ones(points3D.shape[1])) )
x = P[0].project(X)

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(points3D[0],points3D[1],points3D[2],'k.')
plt.show()

可以得到如下的结果：

其次还可以使用旋转工具观察不同视图的三维图像点。

3.计算F：八点法

八点法是通过对应点来计算基础矩阵的算法。外极约束可以写成线性系统的形式：

$\begin{bmatrix}x_2^1x_1^1&x_2^1y_1^1&x_2^1w_1^1&\cdots&w_2^1w_1^1\\x_2^2x_1^2&x_2^2y_1^2&x_2^2w_1^2&\cdots&w_2^2w_1^2\\\vdots&\vdots&\vdots&\ddots&\vdots\\x_2^nx_1^n&x_2^ny_1^n&x_2^nw_1^n&\cdots&w_2^nw_1^n\end{bmatrix}\begin{bmatrix}F_{11}\\F_{12}\\F_{13}\\\vdots\\F_{33}\end{bmatrix}=A\boldsymbol{f}=0$

因为算法中需要 8 个对应点来计算基础矩阵 F，所以该算法叫做八点法。其实现的函数为：

def compute_fundamental(x1,x2):
    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")
    A = np.zeros((n,9))
    for i in range(n):
        A[i] = [x1[0,i]*x2[0,i], x1[0,i]*x2[1,i], x1[0,i]*x2[2,i],
            x1[1,i]*x2[0,i], x1[1,i]*x2[1,i], x1[1,i]*x2[2,i],
            x1[2,i]*x2[0,i], x1[2,i]*x2[1,i], x1[2,i]*x2[2,i] ]
    U,S,V = linalg.svd(A)
    F = V[-1].reshape(3,3)
    U,S,V = linalg.svd(F)
    S[2] = 0
    F = np.dot(U,np.dot(np.diag(S),V))
    return F

通常使用 SVD 算法来计算最小二乘解。但是上面的函数忽略了对图像坐标进行归一化这一步骤，可能会带来数值问题。

4.外极点和外极线

外极点满可以通过计算 F 的零空间来得到，如果想获得另一幅图像的外极点（，只需将 F 转置后输入下面的函数即可：

def compute_epipole(F):
    U,S,V = linalg.svd(F)
    e = V[-1]
    return e/e[2]

def plot_epipolar_line(im,F,x,epipole=None,show_epipole=True):
    m,n = im.shape[:2]
    line = np.dot(F,x)
    # 外极线参数和值
    t = np.linspace(0,n,100)
    lt = np.array([(line[2]+line[0]*tt)/(-line[1]) for tt in t])
    # 仅仅处理位于图像内部的点和线
    ndx = (lt>=0) & (lt<m)
    plt.plot(t[ndx],lt[ndx],linewidth=2)
    if show_epipole:
        if epipole is None:
            epipole = compute_epipole(F)
        plt.plot(epipole[0]/epipole[2],epipole[1]/epipole[2],'r*')

ndx = (corr[:,0]>=0) & (corr[:,1]>=0)
x1 = points2D[0][:,corr[ndx,0]]
x1 = np.vstack( (x1,np.ones(x1.shape[1])) )
x2 = points2D[1][:,corr[ndx,1]]
x2 = np.vstack( (x2,np.ones(x2.shape[1])) )
F = compute_fundamental(x1,x2)
e = compute_epipole(F)

plt.figure()
plt.subplot(221)
plt.imshow(img1)
plt.axis('off')

plt.subplot(222)
plt.imshow(img1)
for i in range(5):
 plot_epipolar_line(img1,F,x2[:,i],e,False)
plt.axis('off')

plt.subplot(223)
plt.imshow(img2)
plt.axis('off')

plt.subplot(224)
plt.imshow(img2)
for i in range(5):
 plt.plot(x2[0,i],x2[1,i],'o')
plt.axis('off')
plt.show()

结果为：

二、照相机和三维结构的计算

1.三角剖分

给定照相机参数模型，图像点可以通过三角剖分来恢复出这些点的三维位置。其基本思想为对于两个照相机 $P_1$ 和 $P_2$ 的视图，三维实物点 X 的投影点为 $x_1$ 和 $x_2$ ，照相机方程定义了如下的关系：

$\begin{bmatrix}P_1&-\mathbf{x}_1&0\\P_2&0&-\mathbf{x}_2\end{bmatrix}\begin{bmatrix}\mathbf{X}\\\lambda_1\\\lambda_2\end{bmatrix}=0$

可以通过 SVD 算法来得到三维点的最小二乘估值。计算一个点对的最小二乘三角剖分的函数为：

def triangulate_point(x1,x2,P1,P2):
    M = np.zeros((6,6))
    M[:3,:4] = P1
    M[3:,:4] = P2
    M[:3,4] = -x1
    M[3:,5] = -x2
    U,S,V = linalg.svd(M)
    X = V[-1,:4]
    return X / X[3]
def triangulate(x1,x2,P1,P2):
    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")
    X = [ triangulate_point(x1[:,i],x2[:,i],P1,P2) for i in range(n)]
    return np.array(X).T

第一个函数为一个点的，其最后一个特征向量的前 4 个值就是齐次坐标系下的对应三维坐标。第二个函数是实现多个点的三角刨分。其输入是两个图像点数组，输出为一个三维坐标数组。下面为一个对前面所做实验的数据集进行三角刨分的实验。

ndx = (corr[:,0]>=0) & (corr[:,1]>=0)
# 获取坐标，并用齐次坐标表示
x1 = points2D[0][:,corr[ndx,0]]
x1 = np.vstack( (x1,np.ones(x1.shape[1])) )
x2 = points2D[1][:,corr[ndx,1]]
x2 = np.vstack( (x2,np.ones(x2.shape[1])) )
Xtrue = points3D[:,ndx]
Xtrue = np.vstack( (Xtrue,np.ones(Xtrue.shape[1])) )
# 检查前三个点
Xest = triangulate(x1,x2,P[0].P,P[1].P)
print (Xest[:,:3])
print (Xtrue[:,:3])
# 绘制图像
from mpl_toolkits.mplot3d import axes3d
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(Xest[0],Xest[1],Xest[2],'ko')
ax.plot(Xtrue[0],Xtrue[1],Xtrue[2],'r.')
plt.axis('auto')
plt.show()

其结果为：

2.由三维点计算照相机矩阵

我们可以使用直接线性变换的方法来计算一些三维点及其图像投影的照相机矩阵P。其本质上为是三角剖分方法的逆问题，有时我们将其称为照相机反切法。可以使用 SVD 分解估计出照相机矩阵。其实现代码为：

def compute_P(x,X): 
    n = x.shape[1]
    if X.shape[1] != n:
        raise ValueError("Number of points don't match.")
    M = np.zeros((3*n,12+n))
    for i in range(n):
        M[3*i,0:4] = X[:,i]
        M[3*i+1,4:8] = X[:,i]
        M[3*i+2,8:12] = X[:,i]
        M[3*i:3*i+3,i+12] = -x[:,i]
    U,S,V = linalg.svd(M)
    return V[-1,:12].reshape((3,4))

下面为一个对前面所做实验的数据集进行估计照相机矩阵的实验：

corr = corr[:,0] # 视图 1
ndx3D = np.where(corr>=0)[0] # 丢失的数值为 -1
ndx2D = corr[ndx3D]
# 选取可见点，并用齐次坐标表示
x = points2D[0][:,ndx2D] # 视图 1
x = np.vstack( (x,np.ones(x.shape[1])) )
X = points3D[:,ndx3D]
X = np.vstack( (X,np.ones(X.shape[1])) )
# 估计 P
Pest = camera.Camera(compute_P(x,X))
# 比较！
print (Pest.P / Pest.P[2,3])
print (P[0].P / P[0].P[2,3])
xest = Pest.project(X)
# 绘制图像
plt.figure()
plt.imshow(img1)
plt.plot(x[0],x[1],'bo')
plt.plot(xest[0],xest[1],'r.')
plt.axis('off')
plt.show()

其结果为：

上面的第一张图是照相机矩阵归一化的格式，其中的上半部分是估计出的照相机矩阵，下面为数据集的创建者计算出的照相机矩阵。第二张图是使用估计出的照相机矩阵投影这些三维点，然后绘制出投影后的结果。

3.由基础矩阵计算照相机矩阵

3.1 未标定的情况——投影重建

在没有任何照相机内参数知识的情况下，照相机矩阵只能通过射影变换恢复出来。在无标定的情况下，第二个照相机矩阵可以使用一个（3×3）的射影变换得出。一个简单的方法为：

$P_2=[S_{\mathrm{e}}F|\mathrm{e}]$

但是使用该矩阵做出的三角形剖分很有可能会发生畸变，其实现代码为：

def compute_P_from_fundamental(F):
    e = compute_epipole(F.T) # 左极点
    Te = skew(e)
    return np.vstack((np.dot(Te,F.T).T,e)).T

def skew(a):
    return np.array([[0,-a[2],a[1]],[a[2],0,-a[0]],[-a[1],a[0],0]])

3.2 已标定的情况——度量重建

在已经标定的情况下，重建会保持欧式空间中的一些度量特性。在标定归一化的坐标系下，基础矩阵称为本质矩阵。为了区别为标定后的情况，以及归一化了的图像坐标，我们通常将其记为 E，而非 F。从本质矩阵恢复出的照相机矩阵中存在度量关系，但有四个可能解。其实现代码为：

def compute_P_from_essential(E):
    U,S,V = np.svd(E)
    if np.det(np.dot(U,V))<0:
        V = -V
    E = np.dot(U,np.dot(np.diag([1,1,0]),V))
    Z = skew([0,0,-1])
    W = np.array([[0,-1,0],[1,0,0],[0,0,1]])
    P2 = [np.vstack((np.dot(U,np.dot(W,V)).T,U[:,2])).T,
            np.vstack((np.dot(U,np.dot(W,V)).T,-U[:,2])).T,
            np.vstack((np.dot(U,np.dot(W.T,V)).T,U[:,2])).T,
            np.vstack((np.dot(U,np.dot(W.T,V)).T,-U[:,2])).T]
    return P2

三、多视图重建

假设照相机已经标定，计算重建可以分为如下的步骤：

检测特征点，然后在两幅图像间匹配
由匹配计算基础矩阵
由基础矩阵计算照相机矩阵
三角剖分这些三维点

1.稳健估计基础矩阵

类似于稳健计算单应性矩阵，当存在噪声和不正确的匹配时，我们需要估计基础矩阵。其实现代码为：

class RansacModel(object):
    def __init__(self,debug=False):
        self.debug = debug
    def fit(self,data):
        data = data.T
        x1 = data[:3,:8]
        x2 = data[3:,:8]
        F = compute_fundamental_normalized(x1,x2)
        return F
    def get_error(self,data,F):
        data = data.T
        x1 = data[:3]
        x2 = data[3:]
        Fx1 = np.dot(F,x1)
        Fx2 = np.dot(F,x2)
        denom = Fx1[0]**2 + Fx1[1]**2 + Fx2[0]**2 + Fx2[1]**2
        err = ( np.diag(np.dot(x1.T,np.dot(F,x2))) )**2 / denom
        return err
def compute_fundamental_normalized(x1,x2):
    n = x1.shape[1]
    if x2.shape[1] != n:
        raise ValueError("Number of points don't match.")
    x1 = x1 / x1[2]
    mean_1 = np.mean(x1[:2],axis=1)
    S1 = np.sqrt(2) / np.std(x1[:2])
    T1 = np.array([[S1,0,-S1*mean_1[0]],[0,S1,-S1*mean_1[1]],[0,0,1]])
    x1 = np.dot(T1,x1)
    x2 = x2 / x2[2]
    mean_2 = np.mean(x2[:2],axis=1)
    S2 = np.sqrt(2) / np.std(x2[:2])
    T2 = np.array([[S2,0,-S2*mean_2[0]],[0,S2,-S2*mean_2[1]],[0,0,1]])
    x2 = np.dot(T2,x2)
    # 使用归一化的坐标计算F
    F = compute_fundamental(x1,x2)
    # 反归一化
    F = np.dot(T1.T,np.dot(F,T2))
    return F/F[2,2]
def F_from_ransac(x1,x2,model,maxiter=5000,match_theshold=1e-6):
    import ransac
    data = np.vstack((x1,x2))
    # 计算F，并返回正确点索引
    F,ransac_data = ransac.ransac(data.T,model,8,maxiter,match_theshold,20,
    return_all=True)
    return F, ransac_data['inliers']

第二个函数的作用是将图像点归一化为零均值固定方差；第三个函数的作用是知道哪些匹配和F矩阵是一致的。

2.示例

代码实现为：


K = np.array([[2394,0,932],[0,2398,628],[0,0,1]])
# 载入图像，并计算特征
im1 = np.array(Image.open(r"D:\test\gbg\DSC_0132.JPG"))
sift.process_image(r"D:\test\gbg\DSC_0132.JPG",'im1.sift')
l1,d1 = sift.read_features_from_file('im1.sift')
im2 = np.array(Image.open(r"D:\test\gbg\DSC_0133.JPG"))
sift.process_image(r"D:\test\gbg\DSC_0133.JPG",'im2.sift')
l2,d2 = sift.read_features_from_file('im2.sift')
# 匹配特征
matches = sift.match_twosided(d1,d2)
ndx = matches.nonzero()[0]
# 使用齐次坐标表示，并使用 inv(K) 归一化
x1 = homography.make_homog(l1[ndx,:2].T)
ndx2 = [int(matches[i]) for i in ndx]
x2 = homography.make_homog(l2[ndx2,:2].T)
x1n = np.dot(inv(K),x1)
x2n = np.dot(inv(K),x2)
# 使用 RANSAC 方法估计 E
model = RansacModel()
E,inliers = F_from_ransac(x1n,x2n,model)
# 计算照相机矩阵（P2 是 4 个解的列表）
P1 = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0]])
P2 = compute_P_from_essential(E)
ind = 0
maxres = 0
for i in range(4):
    X = triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[i])
    d1 = np.dot(P1,X)[2]
    d2 = np.dot(P2[i],X)[2]
    if sum(d1>0)+sum(d2>0) > maxres:
        maxres = sum(d1>0)+sum(d2>0)
        ind = i
        infront = (d1>0) & (d2>0)
    # 三角剖分正确点，并移除不在所有照相机前面的点
    X = triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[ind])
    X = X[:,infront]

fig = plt.figure()
plt.subplot(223)
ax = fig.gca(projection='3d')
ax.plot(-X[0],X[1],X[2],'k.')
plt.axis('off')

plt.subplot(224)
ax = fig.gca(projection='3d')
ax.plot(-X[0],X[1],X[2],'k.')
plt.axis('off')

cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2[ind])
x1p = cam1.project(X)
x2p = cam2.project(X)
x1p = np.dot(K,x1p)
x2p = np.dot(K,x2p)

plt.subplot(221)
plt.imshow(im1)
plt.gray()
plt.plot(x1p[0],x1p[1],'o')
plt.plot(x1[0],x1[1],'r.')
plt.axis('off')

plt.subplot(222)
plt.imshow(im2)
plt.gray()
plt.plot(x2p[0],x2p[1],'o')
plt.plot(x2[0],x2[1],'r.')
plt.axis('off')
plt.show()

其结果为：

很明显的存在问题。

3.扩展

3.1 多视图

当我们有同一场景的多个视图时，三维重建会变得更准确，包括更多的细节信息。对于视频序列，我们可以使用时序关系，在连续的帧对中匹配特征。但是存在误差会在加入的视图中积累的问题。对于静止的图像，一个办法是找到一个中央参考视图，然后计算与之有关的所有其他照相机矩阵。另一个办法是计算一个图像对的照相机矩阵和三维重建，然后增量地加入新的图像和三维点。

3.2光束法平差

多视图重建的最后一步，通常是通过优化三维点的位置和照相机参数来减少二次投影误差。该过程称为光束法平差。

3.3 自标定

在未标定照相机的情形中，有时可以从图像特征来计算标定矩阵。该过程称为自标定。

四、立体图形

一个多视图成像的特殊例子是立体视觉，即使用两台只有水平偏移的照相机观测同一场景。当照相机的位置如上设置，两幅图像具有相同的图像平面，图像的行是垂直对齐的，那么称图像对是经过矫正的。该设置在机器人学中很常见，常被称为立体平台。

立体重建（有时也称为致密深度重建）就是恢复深度图（或者视差图），图像中每个像素的深度（或者视差）都需要计算出来。

计算视差图

在该立体重建算法中，我们将对于每个像素尝试不同的偏移，并按照局部图像周围归一化的互相关值，选择具有最好分数的偏移，然后记录下该最佳偏移。因为每个偏移在某种程度上对应于一个平面，所以该过程有时称为扫平面法。当密集地应用在图像中时，归一化的互相关值可以很快地计算出来。扫平面法的实现代码为：

def plane_sweep_ncc(im_l,im_r,start,steps,wid):
    m,n = im_l.shape
    mean_l = np.zeros((m,n))
    mean_r = np.zeros((m,n))
    s = np.zeros((m,n))
    s_l = np.zeros((m,n))
    s_r = np.zeros((m,n))
    dmaps = np.zeros((m,n,steps))
    filters.uniform_filter(im_l,wid,mean_l)
    filters.uniform_filter(im_r,wid,mean_r)
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    for displ in range(steps):
        filters.uniform_filter(np.roll(norm_l,-displ-start)*norm_r,wid,s)
        filters.uniform_filter(np.roll(norm_l,-displ-start)*np.roll(norm_l,-displ-start),wid,s_l)
        filters.uniform_filter(norm_r*norm_r,wid,s_r)
        dmaps[:,:,displ] = s/np.sqrt(s_l*s_r)
    return np.argmax(dmaps,axis=2)
def plane_sweep_gauss(im_l,im_r,start,steps,wid):
    m,n = im_l.shape
    mean_l = np.zeros((m,n))
    mean_r = np.zeros((m,n))
    s = np.zeros((m,n))
    s_l = np.zeros((m,n))
    s_r = np.zeros((m,n))
    dmaps = np.zeros((m,n,steps))
    filters.gaussian_filter(im_l,wid,0,mean_l)
    filters.gaussian_filter(im_r,wid,0,mean_r)
    norm_l = im_l - mean_l
    norm_r = im_r - mean_r
    for displ in range(steps):
        filters.gaussian_filter(np.roll(norm_l,-displ-start)*norm_r,wid,0,s)
        filters.gaussian_filter(np.roll(norm_l,-displ-start)*np.roll(norm_l,-displ-start),wid, 0,s_l)
        filters.gaussian_filter(norm_r*norm_r,wid,0,s_r) # 和反归一化
        dmaps[:,:,displ] = s/np.sqrt(s_l*s_r)
    return np.argmax(dmaps,axis=2)

im_l = np.array(Image.open(r'D:\test\scene1.row3.col3.ppm').convert('L'),'f')
im_r = np.array(Image.open(r'D:\test\scene1.row3.col4.ppm').convert('L'),'f')
# 开始偏移，并设置步长
steps = 12
start = 4
wid = 9
res = plane_sweep_ncc(im_l,im_r,start,steps,wid)
res2 = plane_sweep_gauss(im_l,im_r,start,steps,wid)
plt.figure()
plt.subplot(221)
plt.axis('off')
plt.imshow(im_l)
plt.subplot(222)
plt.axis('off')
plt.imshow(im_r)
plt.subplot(223)
plt.axis('off')
plt.imshow(res)
plt.subplot(224)
plt.axis('off')
plt.imshow(res2)
plt.show()

左下图为是标准的 NCC 扫平面法重建的视差图，右下图为高斯版本的视差图。与标准版本相比，高斯版本具有较少的噪声，但缺少很多细节信息。

asdaasddsa

关注

18
点赞
踩
14

收藏

觉得还不错? 一键收藏
0
评论
Python计算机视觉编程第五章

在该立体重建算法中，我们将对于每个像素尝试不同的偏移，并按照局部图像周围归一化的互相关值，选择具有最好分数的偏移，然后记录下该最佳偏移。当照相机的位置如上设置，两幅图像具有相同的图像平面，图像的行是垂直对齐的，那么称图像对是经过矫正的。如果有一个场景的两个视图以及视图中的对应图像点，那么根据照相机间的空间相对位置关系、照相机的性质以及三维场景点的位置，可以得到对这些图像点的一些几何关系约束。上面的第一张图是照相机矩阵归一化的格式，其中的上半部分是估计出的照相机矩阵，下面为数据集的创建者计算出的照相机矩阵。
复制链接

扫一扫