无人机航拍图像匹配——ORB算法实践（含代码）

一.摘要
二.ORB算法的原理
- 1. FAST角点检测
- 2.BRIEF描述子
- - 描述子生成
- 3.方向计算（角点检测和描述子之间的步骤）
- 4.尺度金字塔
三.使用opencv-python内置的ORB模块进行图像匹配
参考文献

一.摘要

航拍图像通常涉及不同视角和尺度的变化，这导致相邻图像之间的特征点具有不同的尺度和方向，增加了匹配的困难。特征点提取和描述子生成算法需要具备尺度不变性和旋转不变性，以应对不同尺度和旋转变化的情况，并且航拍图像往往包含大规模的场景，其中可能存在大量的特征点和复杂的纹理。在处理大规模场景时，特征点匹配算法需要具备高效的计算能力和处理速度，以应对大量特征点的匹配需求。

比较通用的特征点提取和描述的算法是SIFT(Scale-Invariant Feature Transform)和ORB(Oriented FAST and Rotated BRIEF）。
SIFT算法是一种基于尺度不变特征的算法，能够在不同尺度和旋转下提取稳定的特征点，并生成描述子进行匹配，适用于各种复杂场景，但是实时性较差
ORB则具有较好的实时性，适用于航拍定位的各种场景

本文分为两部分，第一部分详细讲解ORB算法的原理，第二部分结合代码详细讲解如何使用opencv-python内置的ORB模块进行图像匹配

二.ORB算法的原理

图像匹配的流程是：
角点检测——>用描述子对角点周围图像区域的特征进行编码的向量或二进制表示——>对两张图中的描述子进行匹配——>离群点消除

ORB算法是一种结合了FAST**角点检测**和BRIEF**描述子的特征点提取和匹配算法，并且通过方向计算和构造尺度金字塔**使得其具有一定的旋转，尺度不变性
想要弄清楚ORB为什么快，就必须要知道FAST角点检测，和BRIEF描述子。

1. FAST角点检测

FAST的基本原理就是先在图片中取一个像素点P,其强度为 $I_p$ , 然后用周长为3个像素的圆来判定P点是否为角点，具体如图所示：
FAST原理图
顺时针为圆编号，这里编了16个序号，设置一个阈值 $T$ 如果连续有 N 个点的强度 $I$ 都大于 $I_p+T$ 或者小于 $I_p-T$ ，那么中心像素点P就可以被判定为角点。
为了提升速度，可以先检测 $I$ -1,5,9,13这四个点的强度，如果没有超过三个点强度 $I$ 大于 $I_p+T$ 或者小于 $I_p-T$ ，直接拒绝判定中心像素点P为角点。否则检查所有16个点再进行判定。
上述是判定某个点是否为角点的流程，之后遍历整张图片。

N一般取12，实际上N取9即可获得较好的效果

在SIFT算法中，为了实现旋转不变性，需要对每个关键点计算方向，这需要对关键点周围的像素进行梯度计算。FAST角点检测器使用简单的算法和查找表，避免了复杂的计算，因此速度较快。

2.BRIEF描述子

在ORB特征提取根据FAST角点检测算法检测出特征点之后，我们需要描述这些特征点的属性。对于这些特征点的描述算法，我们称之为特征点的描述子(Feature DescritorS)。ORB特征提取算法采用BRIEF描述子来描述这些角点的属性。

描述子生成

当使用BRIEF算法生成描述子时，具体的步骤如下：

关键点选择：首先，从检测到的关键点中选择一个关键点。
生成像素对：在该关键点周围的邻域内选择一组像素对。通常情况下，这组像素对是在固定位置上选择的，例如，可以使用一个预定义的像素对模板，其中包含要比较的像素对的相对位置。可以根据特定的需求和性能要求自定义这个模板。
灰度值比较：对于每对像素，比较其灰度值。将第一个像素的灰度值与第二个像素的灰度值进行比较。
生成二进制编码：根据比较结果，为每个像素对生成一个二进制位。如果第一个像素的灰度值大于第二个像素的灰度值，则将对应的二进制位设为1，否则设为0。
组合成描述子：将所有生成的二进制位组合成一个固定长度的二进制描述子。这个固定长度通常是根据像素对的数量来确定的。

需要注意的是，BRIEF算法是基于像素对的灰度值比较来生成二进制编码的。为了提高描述子的鲁棒性和区分度，可以采用一些优化技巧，例如使用高斯核函数加权的像素对选择、使用积分图像提高计算效率等。

总的来说就是，每个角点都会生成一个二进制的向量（brief描述子），用来后续的特征匹配。这些描述子具有快速匹配的优势，并且具有一定的旋转不变性和尺度不变性，可以用于特征匹配、目标跟踪、三维重建等计算机视觉任务中。

3.方向计算（角点检测和描述子之间的步骤）

1.对于每个检测到的角点，定义一个固定大小的邻域（例如，以关键点为中心的圆形区域）。
2.在邻域内计算图像的梯度。可以使用Sobel算子来计算梯度的幅值和方向。
3.对于邻域内的每个像素，计算其梯度方向，并统计这些方向的分布。
4.选择主要方向作为关键点的方向。通常，可以选择具有最大梯度幅值或与最大幅值相近的峰值方向作为主要方向。

在ORB算法中，方向计算主要用于两个方面：

关键点的主方向选择：在关键点检测阶段，ORB算法会计算每个关键点周围区域的图像梯度，并统计梯度方向的分布。然后，选择具有最大梯度幅值或与最大幅值相近的峰值方向作为关键点的主方向。这个主方向的选择是为了实现旋转不变性，并在描述子生成时提供参考方向。

描述子的旋转：在生成关键点的描述子时，ORB算法会根据关键点的主方向来调整描述子的旋转。通过将关键点周围的像素点根据主方向进行旋转，可以使得描述子具有旋转不变性。这样，在进行特征匹配时，即使关键点在不同旋转角度下，其描述子仍能准确匹配。

方向计算的目的是为了增强ORB算法的旋转不变性，并在描述子生成和匹配过程中提供方向的一致性。这样，ORB算法可以在不同尺度、旋转和视角变化下提取到稳定的特征，并实现更可靠的特征匹配和识别。（这也是FAST算法不具有的）

4.尺度金字塔

ORB使用高斯金字塔构建图像的尺度空间，通过不同的尺度对图像进行平滑和缩放。这是为了在不同尺度下检测关键点，以实现尺度不变性。(关于尺度金字塔，后续也会更新一篇相关文章，对其进行详细分析)
在尺度空间金字塔的每个层级上，通过比较关键点的响应值与其周围像素的值，来检测关键点。

三.使用opencv-python内置的ORB模块进行图像匹配

（1）导入必要的库，和载入图像，ORB是在灰度图像上进行处理

import cv2
import numpy as np
import time
path_1='xxxxxxxxxxxx/m200.jpg'
path_2='xxxxxxxxxxx/m100.jpg'

image1 = cv2.imread(path_1, 0)
image2 = cv2.imread(path_2, 0)
if image1 is None or image2 is None:
    print("图像文件读取失败")
else:
    print("图像文件读取成功")
# 创建彩色图像副本
color_image1 = cv2.imread(path_1)
color_image2 = cv2.imread(path_2)
if image1 is None or image2 is None:
    print("彩色图像文件读取失败")
else:
    print("彩色图像文件读取成功")

（2）角点检测，描述子生成

orb = cv2.ORB_create()
start_time = time.time()
keypoints1, descriptors1 = orb.detectAndCompute(image1, None)
keypoints2, descriptors2 = orb.detectAndCompute(image2, None)
end_time = time.time()

这里的time是为了计算特征提取耗费的时间
（3）特征匹配

# 暴力匹配
matcher = cv2.DescriptorMatcher_create(cv2.DESCRIPTOR_MATCHER_BRUTEFORCE_HAMMING)
# 特征点匹配
start_time_matching = time.time()
matches = matcher.match(descriptors1, descriptors2)
end_time_matching = time.time()

# 输出特征点检测花费的时间
print("特征点检测花费的时间:", end_time - start_time, "秒")
# 输出特征点匹配花费的时间
print("特征点匹配花费的时间:", end_time_matching - start_time_matching, "秒")

特征点检测花费的时间: 0.24335598945617676 秒
特征点匹配花费的时间: 0.0009953975677490234 秒

这是我的结果，我的机载相机拍照的分辨率是4000×3000，如果图片小一点可能会快很多
（4）初步匹配结果进行离群点去除（滤波滤除不好的匹配点）

ratio_threshold = 0.4
good_matches = []

start_time_filter = time.time()
for match in matches:
    if hasattr(match, 'distance'):
        m = match.distance
        if hasattr(match, 'trainIdx'):
            n = matches[match.trainIdx].distance
            if m < ratio_threshold * n:
                good_matches.append(match)
        else:
            good_matches.append(match)

end_time_filter = time.time()

print("滤波花费的时间:", end_time_filter - start_time_filter, "秒")

这是一个滤波方法，m代表匹配点之间的差异，当然越小越好，通过这个步骤可以剔除离群点。
ratio_threshold一般设置在0.4-0.8之间，如果对精度要求极高，可尽量小。当然，如果设置的太小，可能导致筛选后的匹配点数目稀少，如果少于4对，将会报错。
（4）转移矩阵H

# 提取匹配的特征点的坐标
src_points = np.float32([keypoints1[m.queryIdx].pt for m in good_matches]).reshape(-1, 1, 2)
dst_points = np.float32([keypoints2[m.trainIdx].pt for m in good_matches]).reshape(-1, 1, 2)

# 计算单应矩阵
H, _ = cv2.findHomography(src_points, dst_points, cv2.RANSAC, 3)
#H ,_= cv2.findHomography(src_points, dst_points)
# 输出单应矩阵H
print("Homography Matrix:")
print(H)

这是单应矩阵H，即转移矩阵 $u1,v1,1]^T=H[u2,v2,1]^T$
根据H,没可以获得第一张图中任意一点A(u1，v1)在第二张图中的坐标A’(u2，v2)

（5）误差

# 计算每对匹配点的误差L2
errors = []
right_point=[]
thrshold=10
for match in good_matches:
    kp1 = keypoints1[match.queryIdx]
    kp2 = keypoints2[match.trainIdx]
    pt1 = np.array([kp1.pt[0], kp1.pt[1], 1])
    pt2 = np.array([kp2.pt[0], kp2.pt[1], 1])
    pt1_transformed = np.dot(H, pt1)
    error = np.linalg.norm(pt1_transformed - pt2) ** 2
    if np.linalg.norm(pt1_transformed - pt2)<thrshold :
        right_point.append(np.linalg.norm(pt1_transformed - pt2))
    errors.append(error)

# 计算误差的和L1
L1 = sum(errors)

# 计算均方误差MSE
MSE = np.sqrt(L1) / len(errors)
num_right=len(right_point)
num_all=len(errors)
precision=num_right/num_all

print("MSE:", MSE)
print("总匹配点数:",num_all )
print("正确匹配点数:",num_right )
print("precision:", precision)

将第一张图的特征点通过H投影到第二张上，与其对应的匹配点对比，计算误差。
（6）匹配结果展示

result = cv2.drawMatches(image1, keypoints1, image2, keypoints2, good_matches, None,
                         matchColor=(0, 255, 0), singlePointColor=(0, 0, 255), flags=cv2.DrawMatchesFlags_DEFAULT)
# result = cv2.drawMatches(color_image1, keypoints1, color_image2, keypoints2, good_matches, None,
#                          matchColor=(0, 255, 0), singlePointColor=(0, 0, 255), flags=cv2.DrawMatchesFlags_DEFAULT)

# 调整线条宽度
line_thickness = 2
for match in good_matches:
    pt1 = (int(keypoints1[match.queryIdx].pt[0]), int(keypoints1[match.queryIdx].pt[1]))
    pt2 = (int(keypoints2[match.trainIdx].pt[0]), int(keypoints2[match.trainIdx].pt[1]))
    cv2.line(result, pt1, pt2, (0, 0, 255), thickness=line_thickness)

# 创建匹配结果显示窗口
cv2.namedWindow('Matches', cv2.WINDOW_NORMAL)
cv2.resizeWindow('Matches', 800, 600)
# 显示匹配结果
cv2.imshow('Matches', result)
cv2.waitKey(0)
cv2.destroyAllWindows()

在这里插入图片

参考文献

[1] Viswanathan D G. Features from accelerated segment test (fast)[C]//Proceedings of the 10th workshop on image analysis for multimedia interactive services, London, UK. 2009: 6-8.文章来源地址https://uudwc.com/A/woqN6