基于ORB角点+RANSAC算法实现图像全景拼接

【更新2023/12/7】博客涉及代码部分已开源在github，需要自取，您的star⭐就是对我最大的支持💖：https://github.com/Scienthusiasts/cv-exp/tree/main/exp3

在上一篇博客基于SIFT特征点检测实现图像匹配后，本篇博客将在之前的基础之上实现图像拼接。由于我本人买了盗版的书，书上的字迹不清晰看得难受之外，再加上博主水平有限，不太看得懂也不太想看书上的代码，所以在本篇博客，博主将按照自己理解的思路，手写pipeline还原图像拼接的大致流程。如有任何不严谨，代码冗余或者错误，欢迎在评论区指出，谢谢。

数据采集

本次实战的数据共13张均来自于集美大学诚毅学院(，由本人手持智能手机拍摄获得。

基于ORB角点的图像特征提取与匹配

传统图像拼接算法的第一步无非就是提取两幅图像之间的特征点，然后基于这些特征点进行匹配。

因此我们可以基于上一篇博客的思路，基于SIFT算法实现特征点提取。但是成熟的角点算法多种多样，我们完全可以尝试着换一种角点检测算法作为pipeline的第一步。

ORB的全称是Orinted FAST+Rotated BRIEF 即加了BRIEF描述子的FAST角点检测。对于ORB角点检测算法的详细介绍，感兴趣的小伙伴可以参考我之前的一篇博客【理解】ORB特征提取与ORBSLAM特征匹配简要剖析，这里便不再赘述。值得一提的是，FAST角点检测以快著称，提取速度能够达到SIFT特征的300倍。同样的，BRIEF描述子基于汉明距离计算相似度，相比余弦距离也要快些。

ORB特征提取与匹配的主要方法直接调用OpenCV封装好的函数，代码的大致思路同上一篇博客，没啥好说的。总的需要封装两个函数：

一个就是将提取的角点信息进行格式化，输出为numpy数组：

# 提取图像ORB特征并转化为numpy数组
def extraORBfromImg(ORB, img):
    keypoints, desc = ORB.detectAndCompute(img, mask=None) # 关键点检测
    # 特征点信息
    axis = np.array([kp.pt for kp in keypoints]) # 特征点图像坐标
    scale = np.array([kp.octave+1 for kp in keypoints])# 特征点尺度(在哪一层金字塔)
    direct = np.array([kp.angle*np.pi/180 for kp in keypoints]) # 特征点方向(弧度)
    # 拼接
    infos = np.array([scale,direct]).T
    cors_info = np.hstack([axis,infos])
    return cors_info, desc

另一个就是进行图像间的角点匹配函数：

# ORB特征BRIEF描述子匹配
def ORBMatch(BF, desc1, desc2):
    matches = BF.match(desc1, desc2)
    dist = np.array([mc.distance for mc in matches])
    idx1 = np.array([mc.trainIdx for mc in matches])
    idx0 = np.array([mc.queryIdx for mc in matches])
    idx = np.array([idx0,idx1]).T
    # 匹配点筛选，当描述子之间的距离大于两倍的最小距离时，认为匹配有误
    min_dist = min(dist)
    filte_idx = np.where(dist <= max(2 * min_dist, 30))[0]
    return dist[filte_idx], idx[filte_idx,:]

由于ORBMatch函数返回的是所有点的匹配距离和匹配索引，因此必要时我们还需要一个函数计算匹配点对的坐标：

# 获取图像对匹配点的坐标(一对)
def findMatchCord(match_idx, cors1, cors2):
    left = cors1[match_idx[:,0], :2]
    right = cors2[match_idx[:,1], :2]
    return np.hstack([left, right])

这里放几张ORB角点的匹配结果，额，本来还想找一些误匹配点举例子的，但发现匹配问题不大，效果还行🤣：

盗版书：

校园场景：

标定板(可以发现FAST算法同Harris一样都倾向于检测角点，而SIFT算法倾向于检测斑点)：

基于RANSAC算法剔除误匹配点

得到了角点之后，通常我们无法直接根据匹配坐标直接利用最小二乘方法计算单应矩阵，这是因为图像中可能存在误匹配点，这些点会严重的影响单应矩阵的准确性。

这里就需要利用RANSAC算法来进行匹配点的筛选，(RANSAC全称叫做RAndom SAmple Consensus,随机一致性采样)

算法的基本的假设就是数据是由“内点”和“外点”组成的。“内点”就是组成模型参数的数据，“外点”就是不适合模型的数据。同时RANSAC假设：在给定一组含有少部分“内点”的数据，存在一个程序可以估计出符合“内点”的模型。但这是一种不确定算法，它只能在一种概率下产生结果，并且这个概率会随着迭代次数的增加而加大。

算法流程

RANSAC是通过反复选择数据集去估计出模型，一直迭代到估计出认为比较好的模型。
具体的实现步骤可以分为以下几步：

1. 选择出可以估计出模型的最小数据集；(对于直线拟合来说就是两个点，对于计算Homography矩阵就是4个点)

2. 使用这个数据集来计算出数据模型；

3. 将所有数据带入这个模型，计算出“内点”的数目；(累加在一定误差范围内的适合当前迭代推出模型的数据)

4. 比较当前模型和之前推出的最好的模型的“内点“的数量，记录最大“内点”数的模型参数和“内点”数；

5. 重复1-4步，直到迭代结束或者当前模型已经足够好了(“内点数目大于一定数量”)。

手写RANSAC主函数：

算法的输入是图像的匹配点对，输出为去除掉误匹配的剩余点对

def RANSAC(match_pts):
    pts_num = match_pts.shape[0]
    det_M = 0
    update_match_pts = []
    max_satisfy_rate = 0
    # 最大迭代次数100
    for i in range(100):
        det_M = 0
        while(det_M <= 0.1):
            # 随机选取4对点
            rand4 = np.random.randint(pts_num, size=4)  
            # 基于这4对点计算单应性矩阵
            M, det_M = Homographyfrom4Pts(match_pts[rand4,:])

        # 添加齐次坐标
        homo_pts1 = np.insert(match_pts[:,:2], 2, values=np.ones((1, pts_num)), axis=1).T
        # 重投影齐次坐标
        homo_pts2_hat = (M @ homo_pts1).T
        # 重投影坐标
        pts2_hat = (homo_pts2_hat / homo_pts2_hat[:, 2].reshape(-1,1))[:,:2]
        # 计算误差
        error_matrix = np.sum((match_pts[:,2:4] - pts2_hat)**2, axis=1)
        satisfy_rate = sum(error_matrix < 10) / pts_num
        # 若重投影正确率大于当前最大值, 更新认为是正确的匹配点
        if(satisfy_rate > max_satisfy_rate):
            max_satisfy_rate = satisfy_rate
            update_match_pts =  match_pts[error_matrix < 10]
        # 若重投影正确率大于阈值, 直接返回结果
        if(satisfy_rate > 0.75):
            return update_match_pts
    return update_match_pts

其中，主函数中的Homographyfrom4Pts()用于随机选取的4个点计算单应矩阵：

# 只利用4对点计算单应性矩阵:
def Homographyfrom4Pts(pair_points):
    pt1 = pair_points[:,:2].astype(np.float32)
    pt2 = pair_points[:,2:4].astype(np.float32)
    # 可能的问题出在三点共线或者两点重合的情况，导致误差巨大
    M = cv2.getPerspectiveTransform(pt1, pt2)
    # 返回的行列式用于辅助检查 M 是否正确
    return M, np.linalg.det(M)

值得一提的是，单应矩阵的计算至少需要四个点对也就是八个参数，倘若四个点对中存在重合或共线的点对时，单应矩阵的计算结果就会出错(存在线性无关组)，我们可以通过判断点对组成矩阵的行列式是否等于零来判断点对的选取。或者另一个思路(博客采取这一个)，通过计算M的行列式，来判断单应变换的合理性(一般情况下，单应矩阵应该是满秩的，并且行列式不为负(翻转变换))但更多的情况是，完全共线或者重合的点对几乎不存在，因此我们可以**设置一个阈值，单应矩阵M的行列式小于该阈值就直接认为是错误的(实验中阈值设置为0.1)，并重新选点计算。**而不需要再计算重投影误差。

实验结果：

存在重合的点对：

存在共线的点对：

正确选取的点对：

最后一张结果计算得到了相对正确的映射关系，然而，RANSAC算法考虑到，仅仅四个点对参与计算仍然会存在一定的误差，同时也无法仅通过行列式保证这些匹配点对的映射就是正确的。因此最好的方法是通过这四个点对得到的单应矩阵计算所有匹配点的重投影坐标，筛选出那些误差最小的点对。同时不断迭代，当误差最小的点对占总点对的比值最大时(或大于一个阈值)，就认为得到了正确的映射。最后再利用那些误差最小的点对重新利用最小二乘，计算出噪声较小的单应矩阵。

最小二乘方法，在OpenCV中有现成的函数接口，直接无脑调用就行了。由于接口同样包含了RANSAC方法，为了实验的真实性，这里设置为不调用：

def calc_homography(match_pts):
    src_pts = match_pts[:,:2]
    dst_pts = match_pts[:,2:4]
    # 第三个参数可以使用cv2.RANSAC, 0表示所有点参与计算
    # 第四个参数表示可容忍的重投影误差,范围(1,10)
    # 返回的参数还包含一个mask,参与计算的点数，这里没什么用
    M, _ = cv2.findHomography(src_pts, dst_pts, 0, 10)

    return M

实验结果(绿点为真值，黄点为重投影点,阈值:MSE<10）：

重投影匹配率:0.26

重投影匹配率:0.48

重投影匹配率:0.74

重投影匹配率:1.0

利用最小二乘优化得到的最终映射结果:

值得注意的是，单应变换的基本假设是变换内容属于同一平面上在不同视角下产生的透视变换，然而真实场景(并不满足同一平面的假设，景物必然存在纵深)下必然存在丰富的景深信息，两张图像间存在的视差关系即使是最精确的单应变换也可能无法使得映射后的图像完全重合，因此单纯的单应变换必然也忽略了景深信息。

映射结果可视化

别以为图像全景拼接的难点只在单应矩阵的计算上，事实上可视化出结果同样是具有挑战性的

就拿两张图像的拼接举例，由于透视变换的不规则，上图的拼接结果损失了边角的一些图像信息，因此好的可视化函数需要保留映射图像的所有信息，一个直接的方法就是在M上添加坐标平移，将图像映射到区域的正中心。同时为了保证两张图像映射的位置一致，还需要在某些区域进行黑色填充。除此之外，填充后还需保证两张图像的尺寸一致：

将可视化方法封装成函数(我实现起来可能有些繁琐，欢迎批评指导):

# 可视化图像对映射效果
def homography_trans(M, img1, img2):
    # out_img 第一张图像映射到第二张
    x_min, x_max, y_min, y_max, M2 = calc_border(M, img1.shape)
    # 透视变换+平移变换(使得图像在正中央)
    M = M2 @ M
    out_img = cv2.warpPerspective(img1, M, (round(x_max)-round(x_min), round(y_max)-round(y_min)))
    # 调整两张图像位姿一致：
    # x方向
    out_img_blank_x = np.zeros((out_img.shape[0], abs(round(x_min)), 3)).astype(np.uint8)
    img2_blank_x = np.zeros((img2.shape[0], abs(round(x_min)), 3)).astype(np.uint8)
    if(x_min>0):
        print(1)
        out_img = cv2.hconcat((out_img_blank_x, out_img))
    if(x_min<0):
        print(2)
        img2 = cv2.hconcat((img2_blank_x, img2))
    # y方向
    out_img_blank_y = np.zeros((abs(round(y_min)), out_img.shape[1], 3)).astype(np.uint8)
    img2_blank_y = np.zeros((abs(round(y_min)), img2.shape[1], 3)).astype(np.uint8)
    if(y_min>0):
        print(3)
        out_img = cv2.vconcat((out_img, out_img_blank_y))
    if(y_min<0):
        print(4)
        img2 = cv2.vconcat((img2_blank_y, img2))
    # 调整两张图像尺度一致：
    if(img2.shape[0]<out_img.shape[0]):
        blank_y = np.zeros((out_img.shape[0]-img2.shape[0], img2.shape[1], 3)).astype(np.uint8)
        img2 = cv2.vconcat((img2, blank_y)) 
    else:
        blank_y = np.zeros((img2.shape[0]-out_img.shape[0], out_img.shape[1], 3)).astype(np.uint8)
        out_img = cv2.vconcat((out_img, blank_y)) 
    if(img2.shape[1]<out_img.shape[1]):
        blank_x = np.zeros((img2.shape[0], out_img.shape[1]-img2.shape[1], 3)).astype(np.uint8)
        img2 = cv2.hconcat((img2, blank_x)) 
    else:
        blank_x = np.zeros((out_img.shape[0], img2.shape[1]-out_img.shape[1], 3)).astype(np.uint8)
        out_img = cv2.hconcat((out_img, blank_x))        

    # 叠加
    result = addMatches(out_img, img2)
 
    mask = 255*np.ones(result.shape).astype(np.uint8)
    gray_res = cv2.cvtColor(result, cv2.COLOR_BGR2GRAY)
    mask[gray_res==0]=0

    cv2.imwrite('mask.jpg',mask)
    result[result==0]=255
    cv2.imwrite('result.jpg',result)
    return result, out_img

其中的calc_border()方法计算得到变换后图像的四角坐标(不规则的),以便能够截取完整的图像

# 计算边界
def calc_border(M, shape):
    w, h = shape[1], shape[0]
    pt1 = np.array([[0,0],[w,0],[0,h],[w,h]]).astype(np.float32)
    original_border = np.c_[pt1,[1,1,1,1]]
    #计算透视变换后的图像四个角的坐标
    perspected_border = M @ original_border.T
    perspected_border = perspected_border / perspected_border[2,:]
    x_min = min(perspected_border[0,:])
    x_max = max(perspected_border[0,:])
    y_min = min(perspected_border[1,:])
    y_max = max(perspected_border[1,:])
    pt2 = np.array([[-x_min,-y_min],[w-x_min,-y_min],[-x_min,h-y_min],[w-x_min,h-y_min]]).astype(np.float32)
    # 平移变换(将图像平移至正中央防止负坐标变换后被遮挡的情况)
    M2 = cv2.getPerspectiveTransform(pt1, pt2)
    return x_min, x_max, y_min, y_max, M2

然后是图像拼接后重叠区域的处理，这里直接无脑拼，重叠区域均来自同一张图:

# 图像拼接(无脑)
def addMatches(img1, img2):
    img3 = img1[:,:,:]
    img3[img2==0]=img1[img2==0]
    img3[img2!=0]=img2[img2!=0]
    return img3

部分拼接结果可视化

存在的问题:

两张图像的拼接可能还看不出什么，由于透视本身存在的近大远小的特点，拼接后的图像一般会呈现出强烈的透视关系，越到后面图像扭曲的效果越明显，这种现象是正常的但又十分影响观感，就像下面这样:

多张图像的拼接结果：

一种简单且有效的处理方法是对原始图像添加径向畸变(桶形畸变)(属于相机畸变模型的一种)，进行桶形变换后能够缩减边缘图像匹配的发散，尽量让匹配点的纵坐标的差值减小，缩减透视变换后的变形。

为了方便解释直接搬运别人博客的两张图了，十分的直观：

这里参考博客:图像拼接基本步骤

相机畸变模型还不太了解，就先不搞了😓

如何解决拼接图像不规则边界的现象

拼接图像的矩形化最主要的目的是为了解决不规则边界的问题(通常在畸变图像的拼接后处理)，但某种程度上也能解决多图拼接边缘发散的问题。

目前学术上的主流方法来自于Kaiming He(何恺明就是他)于2013年提出的 rectangling方法(基于seam-carving方法)(说实话这部分内容本人没有进行充分的了解，所以还请各位自行查阅相关文献)

论文链接:Rectangling Panoramic Images via Warping

然而最近(上个月的)国内的一项最新研究表明，在传统方法的基础之上，利用深度网络来预测传统方法中mesh的模式，也能够很好的实现不规则边界拼接图像的矩形化(这部分内容本人也没有进行充分的了解，所以还请各位自行查阅相关文献):

论文链接:Deep Rectangling for Image Stitching: A Learning Baseline

Github链接：https://github.com/nie-lang/DeepRectangling

基于好奇，我在此基础之上配了环境大致看了下这项最新成果的矩形化效果，效果非常的amazing啊，唯二的缺点在于，目前只能矩形化两幅图像的拼接结果，还有就是网络的输出的size是固定的，可能会对原始场景带来一定的变形(也有可能是可以更改网络的参数，这部分代码没细看)

环境的配置比较简单，就不细说了，下面贴几张用自己数据跑的效果：

没了