python图像处理笔记-十一-多视图重建与立体图像

python图像处理笔记-十一-多视图重建与立体图像

多视图重建

由于照相机运动给我们提供了三维结构,所以这样计算三维重建的方法通常称作SFM (Structure from Motion,从运动中恢复结构)。我们假设摄像机已经标定,计算重建可以分为下面四个步骤:

  • 检测特征点,在两幅图像中匹配

  • 由匹配计算出基础矩阵

  • 由基础矩阵计算照相机矩阵

  • 三角形剖分这些三维点

我们前面已经把者四个东西都做过了,但是当图像间的点包含不正确的匹配关系时,需要一个文集爱你方法来估计矩阵。

稳健估计基础矩阵

类似于稳健计算单应性矩阵,存在噪声和不正确匹配的时候,我们需要估计基础矩阵,和前面的方法一样,使用RANSAC,不过这次结合了八点法。

我们写一段代码来实现它:

1
# 归一化的八点法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 = np.compute_fundamental(x1,x2)    # 反归一化    F = np.dot(T1.T,np.dot(F,T2))    return F/F[2,2]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):        """        计算所有对应的x^FFx,并返回误差        """        # 转置,将数据分为两个点集        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[1] ** 2 + Fx2[0] ** 2        err = (np.diag(np.dot(F,x2))) ** 2 /denom        return errdef 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']

三维重建实例:

老规矩,我们在来写个demo跑一下数据:

1
from PCV.geometry import homographyfrom PCV.geometry import sfmfrom PCV.localdescriptors import siftfrom numpy import *from pylab import *from scipy import linalgfrom pylab import *from PIL import Image# 标定矩阵K = array([[2394,0,932],[0,2398,628],[0,0,1]])# 载入图像,并计算特征im1 = array(Image.open('alcatraz1.jpg'))sift.process_image('alcatraz1.jpg','im1.sift')l1,d1 = sift.read_features_from_file('im1.sift')im2 = array(Image.open('alcatraz2.jpg'))sift.process_image('alcatraz2.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 = dot(inv(K),x1)x2n = dot(inv(K),x2)# 使用RANSAC 方法估计Emodel = sfm.RansacModel()E,inliers = sfm.F_from_ransac(x1n,x2n,model)# 计算照相机矩阵(P2 是4 个解的列表)P1 = array([[1,0,0,0],[0,1,0,0],[0,0,1,0]])P2 = sfm.compute_P_from_essential(E)# 选取点在照相机前的解ind = 0maxres = 0for i in range(4):  # 三角剖分正确点,并计算每个照相机的深度  X = sfm.triangulate(x1n[:,inliers],x2n[:,inliers],P1,P2[i])  d1 = dot(P1,X)[2]  d2 = 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 axes3dfig = figure()ax = fig.gca(projection='3d')ax.plot(-X[0],X[1],X[2],'k.')axis('off')# 绘制X 的投影from PCV.geometry import camera# 绘制三维点cam1 = camera.Camera(P1)cam2 = camera.Camera(P2[ind])x1p = cam1.project(X)x2p = cam2.project(X)# 反K 归一化x1p = dot(K,x1p)x2p = dot(K,x2p)figure()imshow(im1)gray()plot(x1p[0],x1p[1],'o')plot(x1[0],x1[1],'r.')axis('off')figure()imshow(im2)gray()plot(x2p[0],x2p[1],'o')plot(x2[0],x2[1],'r.')axis('off')show()

最后跑出来的结果如下:

立体图像

多视图成像的特殊例子是:立体视觉,它使用两台只有水平偏移的相机观测同一场景。

假设两张图片经过了矫正,那么对应点的寻找限制在图像的同一行上。一旦找到对应点,由于深度和偏移是成正比的,所以就可以通过水平偏移来计算深度:

\[ Z = \frac{fb}{x_l - x_r} \]

其中,f是矫正后图像的焦距,b是两个照相机中心之间的距离,\(x_l,x_r\)是左右两幅图像中,对应点的x坐标。分开摄相机中心的距离成基线,如下图所示:

立体重建是从两幅图像中恢复出深度图的work。

计算视差图

作者介绍了一种非常之暴力的方法,就是计算每一点取所有深度的代价,然后对于每个点取代价最小的选项,其中,度量在左图中取\(I_1\)点,右图中取\(I_2\)点的时候的代价的函数如下:

\[ ncc(I_1,I_2)=\frac{\Sigma_x((I_1(x)-\mu_1)(I_2(x)-\mu_2))}{\sqrt{\Sigma_x(I_1(x)-\mu_1)^2\Sigma_x(I_2(x)-\mu_2)^2}} \]

大家应该是学过概率论的,下面再普及几个概念(以听懂为目的,):

  • 数学期望:

    • \(EX = \Sigma P_ix_i\),相当于每件事发生的概率与发生时的值的乘积。
  • 方差:

    • \(DX = (EX^2)-E(X^2)\)
  • 协方差:

    • \(Cov(X,Y) = E[(X-\mu_x)(Y-\mu_y)]\)

    • 协方差越接近0,就说明两个量越相互独立

    • 协方差是正数时,说明两个量正相关,反之则说明他们负相关

  • 相关系数:

    • 相关系数代表的是两个量之间的余弦值,这个值越接近0,就说明越不相关,反之则越相关。

所以上面的式子实际上在计算两个像素点附近的协方差,对图像周围的区域计算协方差能够计算出该点附近范围的相似性。

这个很简单,之前我做过类似的,所有就不在相同的事情上浪费时间了。