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

大泽泽的小可爱

已于 2023-07-02 00:50:29 修改

阅读量3.2k

点赞数 8

分类专栏：图像匹配文章标签：计算机视觉 opencv 人工智能无人机目标跟踪

于 2023-06-08 18:31:31 首次发布

本文链接：https://blog.csdn.net/zhangzhao147/article/details/131101665

版权

图像匹配专栏收录该内容

3 篇文章 1 订阅

订阅专栏

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

一.摘要

SIFT（Scale-Invariant Feature Transform）算法是由David Lowe于1999年提出的一种用于图像处理和计算机视觉中的特征提取和匹配方法。它在航拍图像匹配中具有重要的意义，主要体现在以下几个方面：

尺度不变性：航拍图像通常具有大范围的尺度变化，例如拍摄距离目标较远或较近的情况。SIFT算法通过在不同尺度下检测关键点和描述图像特征，具有尺度不变性，可以有效应对航拍图像的尺度变化。
旋转不变性：航拍图像中的目标可能以不同的姿态出现，例如建筑物在不同方向上的拍摄角度。SIFT算法通过为关键点分配主方向，并构建与旋转无关的特征描述子，具有旋转不变性，可以准确匹配具有不同姿态的目标。
鲁棒性：航拍图像受到多种因素的影响，如光照变化、阴影、云层等。SIFT算法通过在关键点周围计算局部梯度方向直方图构建特征描述子，具有一定的鲁棒性，能够应对光照变化和部分遮挡的情况。
高匹配精度：SIFT算法通过计算特征描述子之间的相似性或距离进行特征匹配，具有较高的匹配精度。在航拍图像匹配中，SIFT算法可以识别和匹配具有相似特征的地物或目标，用于航拍图像的配准、拼接、地图制作等任务。

综上所述，SIFT算法在航拍图像匹配中具有尺度不变性、旋转不变性、鲁棒性和高匹配精度的特点，能够有效地提取和匹配航拍图像中的特征，对于航拍图像的处理和分析具有重要的意义。

本文分为两部分，第一部分详细讲解了SIFT的原理，第二部分利用openCV-python内置的SIFT模块进行图像匹配

二.SIFT算法的原理

SIFT算法的主要步骤如下：

尺度空间极值检测（Scale-space extrema detection）：SIFT算法首先通过在不同尺度上对图像进行高斯平滑来构建尺度空间（scale space）。然后，在尺度空间中检测出关键点，这些关键点对于不同尺度下的图像特征是稳定的。
关键点定位（Keypoint localization）：在尺度空间的极值点检测之后，对每个检测到的极值点进行精确定位，以获得更准确的关键点位置。该步骤使用了DoG（Difference of Gaussians）图像金字塔，通过计算高斯平滑图像的差异来检测局部极值点。
方向分配（Orientation assignment）：为每个关键点分配一个主方向，以提高特征描述子的旋转不变性。在该步骤中，通过对关键点周围的图像梯度方向进行统计，确定主要的梯度方向，并将其作为关键点的方向。
特征描述（Feature description）：在这一步骤中，使用关键点周围的图像区域计算特征描述子，以描述关键点周围的图像结构。SIFT算法使用了局部图像梯度的方向直方图来构建描述子，该描述子具有尺度不变性和旋转不变性。
特征匹配（Feature matching）：通过计算两幅图像中的特征描述子之间的距离或相似性来进行特征匹配。常用的方法是计算描述子之间的欧氏距离或余弦相似度，然后根据距离或相似度进行特征匹配。

1.尺度空间极值检测 &关键点定位

首先要明确，ORB是对图像灰度进行检测，得到一系列的角点，与之不同的额是，sift检测出来用于特征点匹配的特征点，是一个个具有尺度不变性的区域：

上图中的红色圆圈就是一个个具有尺度不变性的区域，后续的特征点匹配，也是对两张图的尺度不变区域进行匹配，那么什么是尺度不变性呢？

尺度不变性&尺度空间

尺度不变性是指在不同尺度下，对象或特征的外观和结构保持不变的属性。
对于SIFT算法，就是对于图像中的某个点，在这个点周围画圆，作为分析区域，希望找到一个 $f$ 函数，与不同尺度的该区域的图像相乘，响应结果能够有一个极值。
如果说响应函数存在一个极值，而不是平坦的，那就可以说这个区域具有尺度不变性（也可以说是一种极值特性）
下图横坐标为针对该区域，图像尺度的不同尺度，纵坐标是响应结果。
在这里插入图片描述
(Witkin, 1983)的那篇论文关于这一块的描述是：通过在所有可能的尺度上搜索稳定的特征，利用尺度的连续函数即尺度空间，可以实现对图像的尺度变化不变的位置的检测。
$\sigma)=G(x, y, \sigma) * I(x, y)$
$I (x, y)$ 是输入图像, $G$ 是一个变尺度高斯函数因此，图像的尺度空间定义为变尺度高斯函数 $G (x, y ， σ)$ 与输入图像 $I (x, y)$ 的卷积函数 $y,\sigma)$ , $G$ 的表达式为：
$\sigma)=\frac{1}{2 \pi \sigma^2} e^{-\left(x^2+y^2\right) / 2 \sigma^2}$
$\sigma$ 为尺度空间因子, $\sigma$ 值越小表示图像被平滑的越少，相应的尺度就越小,而当 $\sigma$ 较大时，图像的平滑效果会更好，但可能会导致细节的丢失。效果如下图所示，尺度从左到右依次增大。相应的图片也变得更粗糙，（可能得放大看才能看出差别）代码在文末 在这里插入图片描述
不同的 $\sigma$ 对应不同的尺度，通过调整 $\sigma$ 可以获得图片任意一点（x，y）在不同尺度下的响应（也就是所谓的尺度空间）
也就是说对于任意一个点，都可以获得其与不同 $G(\sigma)$ 的卷积后的响应，这样对每个点（x,y）都可以画出一个纵轴为响应，横轴为尺度的图片。
在这里插入图片描述

对于任意一点（x,y）,只要不同尺度下的相应图中存在峰值，那么这个点就是我们想要的点。

在整张图片上重复上述步骤，可以获得所有具有尺度不变性的点,如果该点是具有尺度不变性的点，以该点为圆心画一个半径为 $\sqrt(2)*尺度$ 的圆（这就是第一张图中蝴蝶周围圆圈的由来）。

接下来讲述如何利用尺度不变性的原理，通过构建高斯金字塔来实现关键点检测

高斯金字塔

SIFT用来了一种比较高效的方法：DOG高斯差分金字塔
在这里插入图片描述
首先，用不同的尺度的高斯核函数构造尺度空间（图左边所示）
然后，将将较小尺度的图像从较大尺度的图像中减去，得到差分图像（图右边所示）
$\begin{aligned} D(x, y, \sigma) & =(G(x, y, k \sigma)-G(x, y, \sigma)) * I(x, y) \\ & =L(x, y, k \sigma)-L(x, y, \sigma)\end{aligned}$
接下来，在差分金字塔中的每个尺度上，检测局部极值点作为特征点的候选。对每个像素点，在其相邻的3×3×3的邻域（包括当前尺度、上下尺度、左右尺度）内比较其值与邻域内所有其他像素点（26个）的值，找到局部极值点。以确保在尺度空间和二维图像空间都检测到极值点。
在这里插入图片描述

最后，对于检测到的局部极值点，由于前面构建的尺度空间是离散的：
在这里插入图片描述

所以下通过拟合二次曲线来精确定位其准确的位置和尺度。拟合过程使用每个极值点及其相邻像素点的响应值进行插值，得到更精确的特征点位置。

至此就可以检测出图像中的特征点，假设该特征点的尺度为Q,特征点周围画一个圆圈（半径为 $\sqrt(2)$ *S）

懒得画图了，还是展示这张图吧

2.方向分配

在SIFT算法的方向分配步骤中，对于每个关键点，可以按照以下步骤进行方向分配：

确定圆的半径：计算关键点所在的高斯图像的尺度的1.5倍，得到圆的半径。
统计梯度方向和梯度幅值：以关键点为圆心，以半径为半径的圆内的所有像素点。对于每个像素点，计算其梯度方向和梯度幅值。可以使用Sobel等算子计算图像的梯度。
高斯加权：对于圆内的每个像素点，根据其到关键点的距离应用高斯加权。距离关键点越近的像素点，其梯度方向和梯度幅值所占的权重越大。
构建梯度方向直方图：将梯度方向划分为若干个区间，例如36个区间。对于每个像素点，根据其梯度方向的值和加权后的梯度幅值，将其贡献到对应的区间上。
选择主方向：从梯度方向直方图中选择具有最大值的方向作为关键点的主方向。

通过以上步骤，可以为每个关键点分配一个主方向，用于后续的特征描述子计算。这样可以使特征描述子对图像的旋转具有不变性，提高特征匹配的准确性和稳定性。
在这里插入图片描述

3.特征描述

这部分参考的是：
SIFT算法
想要了解更深入的数学原理可以看这篇文章
尺度不变特征转换-SIFT

根据前面得到的特征点的主方向，将坐标轴旋转到主方向，从而实现了旋转不变性
在这里插入图片描述

描述符是与特征点所在的尺度有关的，所以描述特征点是需要在该特征点所在的高斯尺度图像上进行的，在高斯尺度图像上，以特征点为中心，将其附近邻域划分为dxd个子区域，论文中取d=4。σ为相对于特征点所在的高斯金字塔的组的基准层图像的尺度

如下图将区域划分为4×4的子块，对每一子块进行8个方向的直方图统计操作，获得每个方向的梯度幅值，总共可以组成128维描述向量.

在这里插入图片描述

所以SIFT算法的整体流程如下图展示
在这里插入图片描述

概括来说就是：
首先对两张输入的图构建尺度金字塔——>进行极值检测——>得到特征点在二维图像空间的位置，以及尺度空间的位置 $\sigma$ ——>特征点所在的高斯尺度图像上，获取特征点的主方向——>利用尺度信息进行区域归一化——>利用主方向消除旋转歧义——>计算外观直方图，获得描述子

4 .特征匹配`

分别对参考图，和实时图建立特征点的描述子集合。采用欧式距离来度量相似性。
设d1为实时图中的点A的描述子和参考图中与A欧式距离最近的点B之间的距离，d2为参考图中与A欧式距离次近的点Q和A之间的距离。设定了一个阈值T。
若要点A和匹配上，则需要
$d 1 / d 2 < T$
关键点的匹配可以用暴力匹配法，但是计算成本太高，一般采用kd树的数据结构来完成搜索。
比较快的匹配方法有FKNN（Fuzzy K-Nearest Neighbors）和KNN（ K-Nearest Neighbors）等。
kd树也可以融入到FKNN和KNN来加速匹配。

关于kd树可以看这个博主讲的，搜了半天感觉他讲的最清晰：
kd树详解

三.代码

1. 无人机航拍图像匹配

(1)导入必要的库加载图片

import cv2
import numpy as np
import time
path_1='C:/Users/22812/zz/opencv_match/H1_match/m50.jpg'
path_2='C:/Users/22812/zz/opencv_match/H1_match/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）特征提取

    # 创建 SIFT 特征提取器
sift = cv2.xfeatures2d.SIFT_create()

# 特征点检测和描述符提取
start_time = time.time()
keypoints1, descriptors1 = sift.detectAndCompute(image1, None)
keypoints2, descriptors2 = sift.detectAndCompute(image2, None)
end_time = time.time()

（3）特征匹配

# 创建 KNN 匹配器
#matcher = cv2.BFMatcher()#暴力检测要50多s
matcher = cv2.FlannBasedMatcher()
# 特征点匹配
start_time_matching = time.time()
matches = matcher.knnMatch(descriptors1, descriptors2, k=2)  # KNN 匹配，返回最近的两个匹配点
end_time_matching = time.time()

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

（4）滤波

ratio_threshold = 0.8
good_matches = []
start_time_filter = time.time()
for m, n in matches:
    if m.distance < ratio_threshold * n.distance:
        good_matches.append(m)
end_time_filter = time.time()
print("滤波花费的时间:", end_time_filter - start_time_filter, "秒")

用KNN匹配法的匹配结果在matches里，matches包含了一对最邻近匹配点，和一对次邻近匹配点。
m.distance < ratio_threshold * n.distance中m.distance越小，表明最邻近匹配点差异越小，越可靠，所以在对精度要求高的时候可以将ratio_threshold=0.8调至更小，一般最小0.4，再小有可能导致匹配点不够4对，不能求解转移矩阵H.

（5）计算转移矩阵

# 提取匹配的特征点的坐标
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)

（6）误差计算

# 计算每对匹配点的误差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)

（7）可视化

# result = cv2.drawMatches(color_image1, keypoints1, color_image2, keypoints2, good_matches, None)
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()

2.高斯核函数的代码

import numpy as np
import cv2
import matplotlib.pyplot as plt

# 定义高斯核函数
def gaussian_kernel(size, sigma):
    kernel = np.fromfunction(lambda x, y: (1 / (2 * np.pi * sigma**2)) * np.exp(-((x - size//2)**2 + (y - size//2)**2) / (2 * sigma**2)), (size, size))
    return kernel / np.sum(kernel)

# 读取输入图像
input_image = cv2.imread("input_image.jpg", 0)  # 以灰度图像方式读取
input_image = input_image.astype(np.float32)  # 转换为float32类型

# 定义高斯核参数
kernel_size = 11  # 高斯核尺寸
sigmas = [1.0, 1.5, 2.0]  # 不同的高斯核标准差

# 创建子图
plt.figure(figsize=(15, 5))

# 显示原图像
plt.subplot(1, 4, 1)
plt.imshow(input_image, cmap='gray')
plt.title('Original Image')
plt.axis('off')

# 遍历不同的sigma值
for i, sigma in enumerate(sigmas):
    # 创建高斯核
    gaussian_filter = gaussian_kernel(kernel_size, sigma)

    # 进行卷积操作
    convolved_image = cv2.filter2D(input_image, -1, gaussian_filter)

    # 显示模糊图像
    plt.subplot(1, 4, i+2)
    plt.imshow(convolved_image, cmap='gray')
    plt.title(f'Sigma = {sigma}')
    plt.axis('off')

# 调整子图之间的间距
plt.tight_layout()

# 显示图像
plt.show()

参考文献

[1] Lowe D G. Distinctive image features from scale-invariant keypoints[J]. International journal of computer vision, 2004, 60: 91-110.
[2]Lindeberg T. Feature detection with automatic scale selection[J]. International journal of computer vision, 1998, 30(2): 79-116.
[3] 部分图来自深蓝学院网课的ppt