5.1 外极几何
多视图几何是利用在不同视点所拍摄图像间的关系,来研究照相机之间或者特征之 间关系的一门科学。多视 图几何中最重要的内容是双视图几何。如果有一个场景的两个视图以及视图中的对应图像点,那么根据照相机间的空间相 对位置关系、照相机的性质以及三维场景点的位置,可以得到对这些图像点的一些 几何关系约束。我们通过外极几何来描述这些几何关系。
三维场景点X经过4×4的单应性矩阵H变换为HX后,则HX在照相机里得到的图像点和X在照相机 P 里得到的图像点相同。利用照相机方程为
,这里的单应性矩阵通常只做刚体变换,即只通过单应矩阵变换了坐标系。将原点和坐标轴与第一个照相机对齐,则:
和
其中和
是标定矩阵,R 是第二个照相机的旋转矩阵,t是第二个照相机的平移向量。利用这些照相机参数矩 阵,我们可以找到点X的投影点
和
(分别对应于投影矩阵
和
)的关系。
同一个图像点经过不同的投影矩阵产生的不同投影点必须满足:,其中F的表达式为
。矩阵
为反对称矩阵:
矩阵F为基础矩阵。基础矩阵可以由两照相机的参数 矩阵(相对旋转R和平移t)表示。由于det(F)=0,所以基础矩阵的秩小于等于2。
第二个视图中的图像点,利用上述公式, 我们找到对应第一幅图像的一条直线:
其中 是第一幅图像中的一条直线。这条线称为对应于点
的外极线,也就是说,点
在第一幅图像中的对应点一定在这条线上。
5.1.1 一个简单的数据集
我们需要一个带有图像点、三维点和照相机参数矩阵的数据集。 我们这里使用一个牛津多视图数据集。
import numpy as np
from PIL import Image
from numpy import genfromtxt
import camera
import matplotlib.pyplot as plt
# 载入一些图像,并确保关闭图像对象
def load_image(file_path):
img = Image.open(file_path)
arr = np.array(img)
img.close() # 关闭图像对象
return arr
im1 = load_image('images001.jpg')
im2 = load_image('images002.jpg')
# 载入每个视图的二维点到列表中
points2D = []
for i in range(3):
points_file = f'00{i + 1}.corners'
try:
points2D.append(np.loadtxt(points_file).T)
except Exception as e:
print(f"Error loading {points_file}: {e}")
# 载入三维点
try:
points3D = np.loadtxt('p3d').T
except Exception as e:
print(f"Error loading 3D points: {e}")
# 载入对应关系
try:
corr = genfromtxt('nview-corners', dtype='int', missing='*')
except Exception as e:
print(f"Error loading correspondence data: {e}")
# 载入照相机矩阵到 Camera 对象列表中
P = []
for i in range(3):
P_file = f'00{i + 1}.P'
try:
P.append(camera.Camera(np.loadtxt(P_file)))
except Exception as e:
print(f"Error loading camera matrix {P_file}: {e}")
将三维的点投影到一个视图,然后和观测到的图像点比较:
# 将三维点转换成齐次坐标表示,并投影
X = np.vstack((points3D, np.ones(points3D.shape[1])))
x = P[0].project(X)
# 在视图1中绘制点
plt.figure()
plt.imshow(im1)
plt.plot(points2D[0][0], points2D[0][1], '*')
plt.axis('off')
plt.figure()
plt.imshow(im1)
plt.plot(x[0], x[1], 'r.')
plt.axis('off')
plt.show()
运行结果如下
5.1.2 用Matplotlib绘制三维数据
Matplotlib中的mplot3d工具 包可以方便地绘制出三维点、线、等轮廓线、表面以及其他基本图形组件,还可以 通过图像窗口控件实现三维旋转和缩放。可以通过在axes对象中加上projection="3d"关键字实现三维绘图,如下:
from mpl_toolkits.mplot3d import axes3d
import matplotlib.pyplot as plt
import numpy as np
# 创建第一个3D图,绘制测试数据
fig1 = plt.figure()
ax1 = fig1.add_subplot(111, projection="3d")
# 生成三维样本
X, Y, Z = axes3d.get_test_data(0.25)
# 在三维中绘制点
ax1.plot(X.flatten(), Y.flatten(), Z.flatten(), 'o')
ax1.set_title("Test Data Points")
plt.show()
# 创建第二个3D图,绘制随机生成的点
points3D = [np.random.rand(10), np.random.rand(10), np.random.rand(10)]
fig2 = plt.figure()
ax2 = fig2.add_subplot(111, projection='3d')
ax2.plot(points3D[0], points3D[1], points3D[2], 'k.')
ax2.set_title("Random 3D Points")
plt.show()
5.1.3 计算F:八点法
八点法是通过对应点来计算基础矩阵的算法。外极约束可以写成线性系统的形式:
其中,f包含F的元素,和
是一对图像点,共有n对对 应点。基础矩阵中有9个元素,由于尺度是任意的,所以只需要8个方程。因为算 法中需要8个对应点来计算基础矩阵F,所以该算法叫做八点法。以下代码是八点法中最小化
的函数:
import numpy as np
def compute_fundamental(x1,x2):
"""使用归一化的八点算法,从对应点(x1,x2 3×n的数组)中计算基础矩阵
每行由如下构成:
[x'*x,x'*y' x', y'*x, y'*y, y', x, y, 1]"""
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 = np.linalg.svd(A)
F = V[-1].reshape(3, 3)
# 受限F
# 通过将最后一个奇异值置0,使秩为2
U, S, V = np.linalg.svd(F)
S[2] = 0
F = np.dot(U, np.dot(np.diag(S), V))
return F
上面的函数忽略了一个重要的步骤:对图像坐标 进行归一化,这可能会带来数值问题。我们会在后面加以解决。
5.1.4 外极点和外极线
外极点满足Fe1 =0,因此可以通过计算F的零空间来得到。计算的函数为
def compute_epipole(F):
""" 从基础矩阵F中计算右极点(可以使用F.T获得左极点)"""
# 返回F的零空间(Fx = 0)
U, S, V = np.linalg.svd(F)
e = V[-1]
return e / e[2]
我们可以在之前样本数据集的前两个视图上运行这两个函数:
import numpy as np
import Merton1
from PIL import Image
import matplotlib.pyplot as plt
def compute_fundamental(x1,x2):
"""使用归一化的八点算法,从对应点(x1,x2 3×n的数组)中计算基础矩阵
每行由如下构成:
[x'*x,x'*y' x', y'*x, y'*y, y', x, y, 1]"""
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 = np.linalg.svd(A)
F = V[-1].reshape(3, 3)
# 受限F
# 通过将最后一个奇异值置0,使秩为2
U, S, V = np.linalg.svd(F)
S[2] = 0
F = np.dot(U, np.dot(np.diag(S), V))
return F
def compute_epipole(F):
""" 从基础矩阵F中计算右极点(可以使用F.T获得左极点)"""
# 返回F的零空间(Fx = 0)
U, S, V = np.linalg.svd(F)
e = V[-1]
return e / e[2]
def plot_epipolar_line(im,F,x,epipole=None,show_epipole=True):
""" 在图像中,绘制外极点和外极线 F×x=0。F是基础矩阵,x是另一幅图像中的点"""
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*')
corr = np.genfromtxt('nview-corners',dtype='int')
#在前面两个视图中点的索引
ndx = (corr[:,0]>=0) & (corr[:,1]>=0)
#获得坐标,并将其用齐次坐标表示
x1 = Merton1.points2D[0][:,Merton1.corr[ndx,0]]
x1 = np.vstack( (x1,np.ones(x1.shape[1])) )
x2 = Merton1.points2D[1][:,Merton1.corr[ndx,1]]
x2 = np.vstack( (x2,np.ones(x2.shape[1])) )
#计算F
F = compute_fundamental(x1,x2)
#计算极点
e = compute_epipole(F)
#绘制图像
plt.figure()
plt.imshow(Merton1.im1)
#分别绘制每条线,这样会绘制出很漂亮的颜色
for i in range(5):
plot_epipolar_line(Merton1.im1,F,x2[:,i],e,False)
plt.axis('off')
plt.figure()
plt.imshow(Merton1.im2)
#分别绘制每个点,这样会绘制出和线同样的颜色
for i in range(5):
plt.plot(x2[0,i],x2[1,i],'o')
plt.axis('off')
plt.show()
运行结果如图所示:
5.2 照相机和三维结构的计算
5.2.1 三角剖分
对于两个照相机和
的视图,三维实物点X的投影点为
和
(这里用齐次坐 标表示),照相机方程定义了下列关系:
我们可以通过SVD算法来得到三维点的最小二乘估值。下面的函数用于计算一个点对的最小二乘三角剖分。
import Merton1
import numpy as np
import matplotlib.pyplot as plt
import sfm
from mpl_toolkits.mplot3d import axes3d
# 前两个视图中点的索引
ndx = (Merton1.corr[:,0]>=0) & (Merton1.corr[:,1]>=0)
# 获取坐标,并用齐次坐标表示
x1 = Merton1.points2D[0][:,Merton1.corr[ndx,0]]
x1 = np.vstack( (x1,np.ones(x1.shape[1])) )
x2 = Merton1.points2D[1][:,Merton1.corr[ndx,1]]
x2 = np.vstack( (x2,np.ones(x2.shape[1])) )
Xtrue = Merton1.points3D[:,ndx]
Xtrue = np.vstack( (Xtrue,np.ones(Xtrue.shape[1])) )
# 检查前三个点
Xest = sfm.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('equal')
plt.show()
上面的代码首先利用前两个视图的信息来对图像点进行三角剖分,然后把前三个图像点的齐次坐标输出到控制台,最后绘制出恢复的最接近三维图像点。输出到控制台的信息如下所示:
5.2.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.")
# 创建用于计算DLT解的矩阵
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 = np.linalg.svd(M)
return V[-1,:12].reshape((3,4))
import sfm, camera
import Merton1
import numpy as np
import matplotlib.pyplot as plt
corr = Merton1.corr[:,0] # 视图 1
ndx3D = np.where(corr>=0)[0] # 丢失的数值为-1
ndx2D = corr[ndx3D]
# 选取可见点,并用齐次坐标表示
x = Merton1.points2D[0][:,ndx2D] # 视图 1
x = np.vstack( (x,np.ones(x.shape[1])) )
X = Merton1.points3D[:,ndx3D]
X = np.vstack( (X,np.ones(X.shape[1])) )
# 估计P
Pest = camera.Camera(sfm.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(im1)
plt.plot(x[0], x[1], 'bo')
plt.plot(xest[0], xest[1], 'r.')
plt.axis('off')
plt.show()
实验效果如下:
5.2.3 由基础矩阵计算照相机矩阵
在两个视图的场景中,照相机矩阵可以由基础矩阵恢复出来。假设第一个照相机矩阵归一化为P1 =[I|0],现在我们需要计算出第二个照相机矩阵P2 。研究分为两类, 未标定的情况和已标定的情况。
1.为标定的情况——投影重建
在没有任何照相机内参数知识的情况下,照相机矩阵只能通过射影变换恢复出来。 也就是说,如果利用照相机的信息来重建三维点,那么该重建只能由射影变换计算 出来(你可以得到整个投影场景中无畸变的重建点)
在无标定的情况下,第二个照相机矩阵可以使用一个(3×3)的射影变换得 出。一个简单的方法是:
其中,e是左极点,满足,
是反对称矩阵。请注意,使用该矩阵做出的三角形剖分很有可能会发生畸变,如倾斜的重建。
2.已标定的情况——度量重建
在已经标定的情况下,重建会保持欧式空间中的一些度量特性(除了全局的尺度参 数)。给定标定矩阵K,我们可以将它的逆作用于图像点
,因此,在新的图像 坐标系下,照相机方程变为:
在新的图像坐标系下,点同样满足之前的基础矩阵方程:
从本质矩阵恢复出的照相机矩阵中存在度量关系,但有四个可能解。因为只有一个 解产生位于两个照相机前的场景,所以我们可以轻松地从中选出来。
5.3 多视图重建
让我们来看,如何使用上面的理论从多幅图像中计算出真实的三维重建。由于照相机的运动给我们提供了三维结构,所以这样计算三维重建的方法通常称为SfM(从运动中恢复结构)。
假设照相机已经标定,计算重建可以分为下面4个步骤:
(1) 检测特征点,然后在两幅图像间匹配;
(2) 由匹配计算基础矩阵;
(3) 由基础矩阵计算照相机矩阵;
(4) 三角剖分这些三维点。
5.3.1 稳健估计基础矩阵
类似于稳健计算单应性矩阵(3.3节),当存在噪声和不正确的匹配时,我们需要估 计基础矩阵。和前面的方法一样,我们使用RANSAC方法,只不过这次结合了八点算法。
class RansacModel(object):
""" 用从 http://www.scipy.org/Cookbook/RANSAC 下载的 ransac.py 计算基础矩阵的类 """
def __init__(self,debug=False):
self.debug = debug
def fit(self,data):
""" 使用选择的 8 个对应计算基础矩阵"""
# 转置,并将数据分成两个点集
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):
""" 计算所有对应的 x^T F x,并返回每个变换后点的误差"""
#转置,并将数据分成两个点集
data = data.T
x1 = data[:3]
x2 = data[3:]
# 将Sampson 距离用作误差度量
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):
""" 使用归一化的八点算法,由对应点(x1,x2 3×n 的数组)计算基础矩阵"""
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):
""" 使用 RANSAN 方法(ransac.py,来自http://www.scipy.org/Cookbook/RANSAC),
从点对应中稳健地估计基础矩阵F
输入:使用齐次坐标表示的点x1,x2(3×n的数组)"""
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,以及正确点的索引,这样可以知道哪些匹配和F 矩阵是一致的。
5.3.2 三维重建示例
import homography
import sfm
import sift2
import numpy as np
from PIL import Image
import matplotlib.pyplot as plt
import camera
# 标定矩阵
K = np.array([[2394,0,932],[0,2398,628],[0,0,1]])
# 载入图像,并计算特征
im1 = np.array(Image.open('alcatraz1.jpg'))
sift2.process_image('alcatraz1.jpg','im1.sift')
l1,d1 = sift2.read_features_from_file('im1.sift')
im2 = np.array(Image.open('alcatraz2.jpg'))
sift2.process_image('alcatraz2.jpg','im2.sift')
l2,d2 = sift2.read_features_from_file('im2.sift')
# 匹配特征
matches = sift2.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(np.linalg.inv(K),x1)
x2n = np.dot(np.linalg.inv(K),x2)
# 使用RANSAC方法估计E
model = sfm.RansacModel()
E,inliers = sfm.F_from_ransac(x1n,x2n,model)
# 计算照相机矩阵(P2是4个解的列表)
P1 = np.array([[1,0,0,0],[0,1,0,0],[0,0,1,0]])
P2 = sfm.compute_P_from_essential(E)
# 选取点在照相机前的解
ind = 0
maxres = 0
infront = 0
for i in range(4):
# 三角剖分正确点,并计算每个照相机的深度
X = sfm.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 = sfm.triangulate(x1n[:, inliers], x2n[:, inliers], P1, P2[ind])
X = X[:, infront]
# 绘制三维图像
from mpl_toolkits.mplot3d import axes3d
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot(-X[0], X[1], X[2], 'k.')
plt.axis('off')
# 绘制X的投影
# 绘制三维点
cam1 = camera.Camera(P1)
cam2 = camera.Camera(P2[ind])
x1p = cam1.project(X)
x2p = cam2.project(X)
# 反K归一化
x1p = np.dot(K, x1p)
x2p = np.dot(K, x2p)
plt.figure()
plt.imshow(im1)
plt.gray()
plt.plot(x1p[0], x1p[1], 'o')
plt.plot(x1[0], x1[1], 'r.')
plt.axis('off')
plt.figure()
plt.imshow(im2)
plt.gray()
plt.plot(x2p[0], x2p[1], 'o')
plt.plot(x2[0], x2[1], 'r.')
plt.axis('off')
plt.show()
5.3.3 多视图的扩展示例
1.多视图
当我们有同一场景的多个视图时,三维重建会变得更准确,包括更多的细节信息。
对于视频序列,我们可以使用时序关系,在连续的帧对中匹配特征。相对方位需要 从每个帧对增量地加入下一个帧对,该方法非常有效,同时可以使用跟踪有效地找到对应点。一个问题是误差会在加入的视图中积累。
对于静止的图像,一个办法是找到一个中央参考视图,然后计算与之有关的所有其 他照相机矩阵。另一个办法是计算一个图像对的照相机矩阵和三维重建,然后增量 地加入新的图像和三维点。
2.光束法平差
以三维重建为例,恢复出的点的位置和由估计 的基础矩阵计算出的照相机矩阵都存在误差。在多个视图的计算中,这些误差会进 一步累积。因此,多视图重建的最后一步,通常是通过优化三维点的位置和照相机 参数来减少二次投影误差。该过程称为光束法平差。
3.自标定
在未标定照相机的情形中,有时可以从图像特征来计算标定矩阵。该过程称为自标 定。根据在照相机标定矩阵参数上做出的假设,以及可用的图像数据的类型(特征 匹配、平行线、平面等),我们有很多不同的自标定算法。