【图像配准】多图配准/不同特征提取算法/匹配器比较测试

前言

本文首先完成之前专栏前置博文未完成的多图配准拼接任务，其次对不同特征提取器/匹配器效率进行进一步实验探究。

各类算法原理简述

看到有博文[1]指出，在速度方面SIFT<SURF<BRISK<FREAK<ORB，在对有较大模糊的图像配准时，BRISK算法在其中表现最为出色，后面考虑选取其中SIFT、BRISK、ORB三种算法进行验证。
在此之前，先对后续算法的原理做一些初步了解。

SIFT算法

在前文【图像配准】SIFT算法原理及二图配准拼接已经对此做过分析，这里不作赘述。

BRISK算法

BRISK算法是2011年ICCV上《BRISK:Binary Robust Invariant Scalable Keypoints》文章中，提出来的一种特征提取算法。

BRISK算法通过利用简单的像素灰度值比较，进而得到一个级联的二进制比特串来描述每个特征点，之后采用了邻域采样模式，即以特征点为圆心，构建多个不同半径的离散化Bresenham同心圆，然后再每一个同心圆上获得具有相同间距的N个采样点。

更详细的内容可参考文献[3]对论文的解读。

ORB算法

ORB(Oriented FAST and rotated BRIEF)是OpenCV实验室开发的一种特征检测与特征描述算法，将 FAST 特征检测与 BRIEF 特征描述结合并进行了改进，具有尺度不变性和旋转不变性，对噪声有较强的抗干扰能力[4]。

ORB算法在图像金字塔中使用FAST算法检测关键点，通过一阶矩计算关键点的方向，使用方向校正的BRIEF生成特征描述符。

更详细的内容可参考文献[4]。

AKAZE算法

Alcantarilla等人提出了AKAZE(Accelerated-KAZE)算法，即加速KAZE算法，加速了非线性尺度空间的构造，效率较KAZE有所提升，以各向异性的非线性滤波来构造尺度空间，将整个尺度空间进行分割，利用局部自适应分级获得细节和噪声，保留较多的边缘细节信息，但该算法关键点检测能力不足，且鲁棒性不强[5]。

多图配准

无论何种算法，图像配准无非是这样几个步骤->图像灰度化->提取特征->构建匹配器->计算变换矩阵->图像合并。
那么多图配准，实际上可以分解为多个双图配准。
以下代码主要参考了这个仓库：https://github.com/799034552/concat_pic

下面按处理顺序对各部分内容进行分块拆解：

图像读取

首先是读取图像再进行灰度化转换。
这里进行了一个判断，判断传入的是否是图像的文本路径，这一步主要是为了后面多图拼接的便利性，因为后面多图拼接会把拼接好的部分图像直接放在内存中，这里若不是路径，就直接赋值给变量，相当于用整张大图去和另外一张小图去做拼接。

# 读取图像-转换灰度图用于检测
# 这里做一个文本判断是为了后面多图拼接处理
if isinstance(path2, str):imageA = cv2.imread(path2)
else:imageA = path2
if isinstance(path1, str):imageB = cv2.imread(path1)
else:imageB = path1
imageA_gray = cv2.cvtColor(imageA, cv2.COLOR_BGR2GRAY)
imageB_gray = cv2.cvtColor(imageB, cv2.COLOR_BGR2GRAY)

构建特征提取器

OpenCV对各种算法都进行了较好的封装，这里主要对比测试了sift，brisk，orb，akaze这几种算法，所用opencv-python版本为4.7.0，值得注意的是，OpenCV4以后的版本，cv2.SURF_create()无法使用，只能用老版本的cv2.xfeatures2d.SURF_create()来实现SURF，因此这里没有对SURF算法进行比较测试。

# 选择特征提取器函数
def detectAndDescribe(image, method=None):if method == 'sift':descriptor = cv2.SIFT_create()elif method == 'surf':descriptor = cv2.xfeatures2d.SURF_create()  # OpenCV4以上不可用elif method == 'brisk':descriptor = cv2.BRISK_create()elif method == 'orb':descriptor = cv2.ORB_create()elif method == 'akaze':descriptor = cv2.AKAZE_create()(kps, features) = descriptor.detectAndCompute(image, None)return kps, features

提取特征/特征匹配

# 提取两张图片的特征
kpsA, featuresA = detectAndDescribe(imageA_gray, method=feature_extractor)
kpsB, featuresB = detectAndDescribe(imageB_gray, method=feature_extractor)
# 进行特征匹配
if feature_matching == 'bf':matches = matchKeyPointsBF(featuresA, featuresB, method=feature_extractor)
elif feature_matching == 'knn':matches = matchKeyPointsKNN(featuresA, featuresB, ratio=0.75, method=feature_extractor)if len(matches) < 10:return None, None

这里比较了两种匹配器，一种是暴力匹配器(BFMatcher)，函数接口为cv2.BFMatcher，主要有下面两个参数可以设置：

• normType：距离类型，可选项，默认为欧式距离NORM_L2。
  • NORM_L1：L1范数，曼哈顿距离。
  • NORM_L2：L2范数，欧式距离。
  • NORM_HAMMING：汉明距离。
  • NORM_HAMMING2：汉明距离2，对每2个比特相加处理。
• crossCheck：交叉匹配选项，可选项，默认为False，若为True，即两张图像中的特征点必须互相都是唯一选择

注：对于SIFT、SURF描述符，推荐选择欧氏距离L1和L2范数；对于ORB、BRISK、BRIEF描述符，推荐选择汉明距离NORM_HAMMING；对于ORB描述符，当WTA_K=3或4时，推荐使用汉明距离NORM_HAMMING2。

对于该函数更详细的内容，可参考博文[6]。

另一个是FLANN匹配器，Flann-based matcher 使用快速近似最近邻搜索算法寻找，FlannBasedMatcher接受两个参数：index_params和search_params：

• index_params：可用不同的数值表示不同的算法，有下表这些可选项(表中数据来源文章[7])

• search_params(int checks=32, float eps=0, bool sorted=true)
  checks为int类型，是遍历次数，一般只改变这个参数

# 创建匹配器
def createMatcher(method, crossCheck):"""不同的方法创建不同的匹配器参数，参数释义BFMatcher：暴力匹配器NORM_L2-欧式距离NORM_HAMMING-汉明距离crossCheck-若为True，即两张图像中的特征点必须互相都是唯一选择"""if method == 'sift' or method == 'surf':bf = cv2.BFMatcher(cv2.NORM_L2, crossCheck=crossCheck)elif method == 'orb' or method == 'brisk' or method == 'akaze':# 创建BF匹配器# bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=crossCheck)index_params = dict(algorithm=1, trees=5)search_params = dict(checks=50)# 创建Flann匹配器bf = cv2.FlannBasedMatcher(index_params, search_params)return bf

二者的区别在于BFMatcher总是尝试所有可能的匹配，从而使得它总能够找到最佳匹配，这也是Brute Force（暴力法）的原始含义。
而FlannBasedMatcher中FLANN的含义是Fast Library forApproximate Nearest Neighbors，从字面意思可知它是一种近似法，算法更快但是找到的是最近邻近似匹配，所以当我们需要找到一个相对好的匹配但是不需要最佳匹配的时候往往使用FlannBasedMatcher。当然也可以通过调整FlannBasedMatcher的参数来提高匹配的精度或者提高算法速度，但是相应地算法速度或者算法精度会受到影响[8]。

特征匹配也有两种方式，可以直接进行暴力检测，也可以采用KNN进行检测，不同检测方式的代码如下：

# 暴力检测函数
def matchKeyPointsBF(featuresA, featuresB, method):start_time = time.time()bf = createMatcher(method, crossCheck=True)best_matches = bf.match(featuresA, featuresB)rawMatches = sorted(best_matches, key=lambda x: x.distance)print("Raw matches (Brute force):", len(rawMatches))end_time = time.time()print("暴力检测共耗时" + str(end_time - start_time))return rawMatches# 使用knn检测函数
def matchKeyPointsKNN(featuresA, featuresB, ratio, method):start_time = time.time()bf = createMatcher(method, crossCheck=False)# rawMatches = bf.knnMatch(featuresA, featuresB, k=2)# 上面这行在用Flann时会报错rawMatches = bf.knnMatch(np.asarray(featuresA, np.float32), np.asarray(featuresB, np.float32), k=2)matches = []for m, n in rawMatches:if m.distance < n.distance * ratio:matches.append(m)print(f"knn匹配的特征点数量:{len(matches)}")end_time = time.time()print("KNN检测共耗时" + str(end_time - start_time))return matches

计算视角变换矩阵/透视变换

匹配完关键点后，就可以计算视角变换矩阵，然后一幅图不动，另一幅图进行透视变换，这里的具体方式和前文较为类似。

# 计算视角变换矩阵
def getHomography(kpsA, kpsB, matches, reprojThresh):start_time = time.time()# 将各关键点保存为ArraykpsA = np.float32([kp.pt for kp in kpsA])kpsB = np.float32([kp.pt for kp in kpsB])# 如果匹配点大于四个点，再进行计算if len(matches) > 4:# 构建出匹配的特征点ArrayptsA = np.float32([kpsA[m.queryIdx] for m in matches])ptsB = np.float32([kpsB[m.trainIdx] for m in matches])# 计算视角变换矩阵(H, status) = cv2.findHomography(ptsA, ptsB, cv2.RANSAC, reprojThresh)end_time = time.time()print("透视关系计算共耗时" + str(end_time - start_time))return matches, H, statuselse:return NoneM = getHomography(kpsA, kpsB, matches, reprojThresh=4)if M is None:print("Error!")(matches, H, status) = M# 将图片A进行透视变换result = cv2.warpPerspective(imageA, H, ((imageA.shape[1] + imageB.shape[1]) * 2, (imageA.shape[0] + imageB.shape[0]) * 2))resultAfterCut = cutBlack(result)

图片黑边裁剪

在做透视变换时，往往会采取一个比较大的背景，以确保图片能够不遗漏的拼接上去，比如这里图片的尺寸设定为(imageA.shape[1] + imageB.shape[1]) * 2, (imageA.shape[0] + imageB.shape[0]) * 2)，这样会产生一些背景黑边，需要进行裁切。

之前的文章提到过一种通过膨胀方式来找到最大内接矩形，这里的代码处理方式更为巧妙，直接采用像素点搜索的方式，找到图像的最大外接矩形。

# 去除图像黑边
def cutBlack(pic):rows, cols = np.where(pic[:, :, 0] != 0)min_row, max_row = min(rows), max(rows) + 1min_col, max_col = min(cols), max(cols) + 1return pic[min_row:max_row, min_col:max_col, :]

图片位置检查

由于无法提前知道两张图片的位置关系，对于透视变换，可能图片会映射到整个选取区域的左边，这样的话，无法正常显示图片，因此，要对透视变换后的图片进行面积检查，如果比原来的图片面积小太多，就用另一张图片来进行透视变换[9]。

if np.size(resultAfterCut) < np.size(imageA) * 0.95:print("图片位置不对,将自动调换")# 调换图片kpsA, kpsB = swap(kpsA, kpsB)imageA, imageB = swap(imageA, imageB)if feature_matching == 'bf':matches = matchKeyPointsBF(featuresB, featuresA, method=feature_extractor)elif feature_matching == 'knn':matches = matchKeyPointsKNN(featuresB, featuresA, ratio=0.75, method=feature_extractor)if len(matches) < 10:return None, NonematchCount = len(matches)M = getHomography(kpsA, kpsB, matches, reprojThresh=4)if M is None:print("Error!")(matches, H, status) = Mresult = cv2.warpPerspective(imageA, H,((imageA.shape[1] + imageB.shape[1]) * 2, (imageA.shape[0] + imageB.shape[0]) * 2))

图像融合

图像融合这里处理得也比较巧妙，对图片接壤部分选取最大值，这样确保了色调的统一性。

# 合并图片-相同的区域选取最大值，从而实现融合result[0:imageB.shape[0], 0:imageB.shape[1]] = np.maximum(imageB, result[0:imageB.shape[0], 0:imageB.shape[1]])result = cutBlack(result)  # 结果去除黑边

多图拼接

最后是拼接多幅图像，反复调用拼接双图即可。

#  合并多张图
def handleMulti(*args):l = len(args)assert (l > 1)# isHandle用于标记图片是否参与合并isHandle = [0 for i in range(l - 1)]nowPic = args[0]args = args[1:]for j in range(l - 1):isHas = False  # 在一轮中是否找到matchCountList = []resultList = []indexList = []for i in range(l - 1):if isHandle[i] == 1:continueresult, matchCount = handle(nowPic, args[i])if not result is None:matchCountList.append(matchCount)  # matchCountList存储两图匹配的特征点resultList.append(result)indexList.append(i)isHas = Trueif not isHas:  # 一轮找完都没有可以合并的return Noneelse:index = matchCountList.index(max(matchCountList))nowPic = resultList[index]isHandle[indexList[index]] = 1print(f"合并第{indexList[index] + 2}个")return nowPic

完整代码

utils.py

import cv2
import numpy as np
import time# 选择特征提取器函数
def detectAndDescribe(image, method=None):if method == 'sift':descriptor = cv2.SIFT_create()elif method == 'surf':descriptor = cv2.xfeatures2d.SURF_create()  # OpenCV4以上不可用elif method == 'brisk':descriptor = cv2.BRISK_create()elif method == 'orb':descriptor = cv2.ORB_create()elif method == 'akaze':descriptor = cv2.AKAZE_create()(kps, features) = descriptor.detectAndCompute(image, None)return kps, features# 创建匹配器
def createMatcher(method, crossCheck):"""不同的方法创建不同的匹配器参数，参数释义BFMatcher：暴力匹配器NORM_L2-欧式距离NORM_HAMMING-汉明距离crossCheck-若为True，即两张图像中的特征点必须互相都是唯一选择"""if method == 'sift' or method == 'surf':bf = cv2.BFMatcher(cv2.NORM_L2, crossCheck=crossCheck)elif method == 'orb' or method == 'brisk' or method == 'akaze':# 创建BF匹配器# bf = cv2.BFMatcher(cv2.NORM_HAMMING, crossCheck=crossCheck)index_params = dict(algorithm=1, trees=5)search_params = dict(checks=50)# 创建Flann匹配器bf = cv2.FlannBasedMatcher(index_params, search_params)return bf# 暴力检测函数
def matchKeyPointsBF(featuresA, featuresB, method):start_time = time.time()bf = createMatcher(method, crossCheck=True)best_matches = bf.match(featuresA, featuresB)rawMatches = sorted(best_matches, key=lambda x: x.distance)print("Raw matches (Brute force):", len(rawMatches))end_time = time.time()print("暴力检测共耗时" + str(end_time - start_time))return rawMatches# 使用knn检测函数
def matchKeyPointsKNN(featuresA, featuresB, ratio, method):start_time = time.time()bf = createMatcher(method, crossCheck=False)# rawMatches = bf.knnMatch(featuresA, featuresB, k=2)# 上面这行在用Flann时会报错rawMatches = bf.knnMatch(np.asarray(featuresA, np.float32), np.asarray(featuresB, np.float32), k=2)matches = []for m, n in rawMatches:if m.distance < n.distance * ratio:matches.append(m)print(f"knn匹配的特征点数量:{len(matches)}")end_time = time.time()print("KNN检测共耗时" + str(end_time - start_time))return matches# 计算视角变换矩阵
def getHomography(kpsA, kpsB, matches, reprojThresh):start_time = time.time()# 将各关键点保存为ArraykpsA = np.float32([kp.pt for kp in kpsA])kpsB = np.float32([kp.pt for kp in kpsB])# 如果匹配点大于四个点，再进行计算if len(matches) > 4:# 构建出匹配的特征点ArrayptsA = np.float32([kpsA[m.queryIdx] for m in matches])ptsB = np.float32([kpsB[m.trainIdx] for m in matches])# 计算视角变换矩阵(H, status) = cv2.findHomography(ptsA, ptsB, cv2.RANSAC, reprojThresh)end_time = time.time()print("透视关系计算共耗时" + str(end_time - start_time))return matches, H, statuselse:return None# 去除图像黑边
def cutBlack(pic):rows, cols = np.where(pic[:, :, 0] != 0)min_row, max_row = min(rows), max(rows) + 1min_col, max_col = min(cols), max(cols) + 1return pic[min_row:max_row, min_col:max_col, :]# 交换
def swap(a, b):return b, a

main.py

from utils import *# 合并两张图（合并多张图基于此函数）
def handle(path1, path2):# 超参数-选择具体算法feature_extractor = 'brisk'feature_matching = 'knn'# 读取图像-转换灰度图用于检测# 这里做一个文本判断是为了后面多图拼接处理if isinstance(path2, str):imageA = cv2.imread(path2)else:imageA = path2if isinstance(path1, str):imageB = cv2.imread(path1)else:imageB = path1imageA_gray = cv2.cvtColor(imageA, cv2.COLOR_BGR2GRAY)imageB_gray = cv2.cvtColor(imageB, cv2.COLOR_BGR2GRAY)# 提取两张图片的特征kpsA, featuresA = detectAndDescribe(imageA_gray, method=feature_extractor)kpsB, featuresB = detectAndDescribe(imageB_gray, method=feature_extractor)# 进行特征匹配if feature_matching == 'bf':matches = matchKeyPointsBF(featuresA, featuresB, method=feature_extractor)elif feature_matching == 'knn':matches = matchKeyPointsKNN(featuresA, featuresB, ratio=0.75, method=feature_extractor)if len(matches) < 10:return None, None# 计算视角变换矩阵matchCount = len(matches)M = getHomography(kpsA, kpsB, matches, reprojThresh=4)if M is None:print("Error!")(matches, H, status) = M# 将图片A进行透视变换result = cv2.warpPerspective(imageA, H, ((imageA.shape[1] + imageB.shape[1]) * 2, (imageA.shape[0] + imageB.shape[0]) * 2))resultAfterCut = cutBlack(result)# 查看裁剪完黑边后的图片# cv2.imshow("resultAfterCut", resultAfterCut)# cv2.waitKey(0)if np.size(resultAfterCut) < np.size(imageA) * 0.95:print("图片位置不对,将自动调换")# 调换图片kpsA, kpsB = swap(kpsA, kpsB)imageA, imageB = swap(imageA, imageB)if feature_matching == 'bf':matches = matchKeyPointsBF(featuresB, featuresA, method=feature_extractor)elif feature_matching == 'knn':matches = matchKeyPointsKNN(featuresB, featuresA, ratio=0.75, method=feature_extractor)if len(matches) < 10:return None, NonematchCount = len(matches)M = getHomography(kpsA, kpsB, matches, reprojThresh=4)if M is None:print("Error!")(matches, H, status) = Mresult = cv2.warpPerspective(imageA, H,((imageA.shape[1] + imageB.shape[1]) * 2, (imageA.shape[0] + imageB.shape[0]) * 2))# 合并图片-相同的区域选取最大值，从而实现融合result[0:imageB.shape[0], 0:imageB.shape[1]] = np.maximum(imageB, result[0:imageB.shape[0], 0:imageB.shape[1]])result = cutBlack(result)  # 结果去除黑边return result, matchCount#  合并多张图
def handleMulti(*args):l = len(args)assert (l > 1)# isHandle用于标记图片是否参与合并isHandle = [0 for i in range(l - 1)]nowPic = args[0]args = args[1:]for j in range(l - 1):isHas = False  # 在一轮中是否找到matchCountList = []resultList = []indexList = []for i in range(l - 1):if isHandle[i] == 1:continueresult, matchCount = handle(nowPic, args[i])if not result is None:matchCountList.append(matchCount)  # matchCountList存储两图匹配的特征点resultList.append(result)indexList.append(i)isHas = Trueif not isHas:  # 一轮找完都没有可以合并的return Noneelse:index = matchCountList.index(max(matchCountList))nowPic = resultList[index]isHandle[indexList[index]] = 1print(f"合并第{indexList[index] + 2}个")return nowPicif __name__ == "__main__":start_time_all = time.time()# 传入图片路径列表，既可以处理两张，也可以处理多张result = handleMulti("./input/foto2B.jpg", "./input/foto2A.jpg")if not result is None:cv2.imwrite("output.tif", result[:, :, [0, 1, 2]])else:print("没有找到对应特征点,无法合并")end_time_all = time.time()print("共耗时" + str(end_time_all - start_time_all))

实验结果

拿原仓库中的两张无人机图片进行了测试，拼接效果如下：

两张原图：

拼接后的图像：

此外，我选取了更大分辨率(4k x 7k)的图像进行拼接测试，比较不同算法的所用时间，结果如下表所示：

特征提取算法	匹配器	特征点个数	时间(s)
sift	bf	14438	463
brisk	bf	9648	31.83
orb	bf	109	20.57
akaze	bf	4872	26.58
brisk	flann	5000	24.71
orb	flann	50	22.02

结论

经过此番实验，可以发现：

1. 从速度上来说orb算法是最快的，比sift这种古老的算法快了一个数量级。但是通过观察生成的图像质量会发现，orb的图像会比较模糊，拼接质量不如其它算法高，增加速度的同时会牺牲部分质量。
2. akaze算法速度和质量和brisk相差不大
3. flann匹配器比bf匹配器通常情况下速度更快

因此，后续实验可以首选brisk算法+flann匹配器的组合方式。
另外说明，上面这些实验参数并没有针对性的进行调参，基本使用默认参数；若进行调优，可能会结果会发生一定变化。

Todo

• 此示例中，默认图像位置是未知的，而在遥感图像中，可以通过gps坐标来确定图像的大致方位，后续考虑引进gps坐标，构建图像排布坐标系，从而加快配准速度。
• 此示例中，多图拼接是直接用大图和小图去做配准，效率并不是太高。后续可能可以结合gps信息，从大图中挖出一部分小图来做配准。

参考文献

[1] BRISK特征点描述算法详解 https://blog.csdn.net/weixin_41063476/article/details/90407916
[2] 基于视觉的特征匹配算法（持续更新）https://zhuanlan.zhihu.com/p/147325381?ivk_sa=1024320u
[3] BRISK算法学习笔记 https://blog.csdn.net/weixin_40196271/article/details/84143545
[4] 【OpenCV 例程 300篇】246. 特征检测之ORB算法 https://blog.csdn.net/youcans/article/details/128033070
[5] 一种无人机滑坡遥感影像的快速匹配算法 https://mp.weixin.qq.com/s?__biz=MzIxOTY4NDQ1MA==&mid=2247493784&idx=1&sn=65676fc368e6b4fa62c965a996f956ef
[6]【OpenCV 例程 300篇】251. 特征匹配之暴力匹配 https://blog.csdn.net/youcans/article/details/128253435
[7] 以图搜图–基于FLANN特征匹配 https://zhuanlan.zhihu.com/p/520575652
[8] [opencv] BF匹配器和Flann匹配器 https://blog.csdn.net/simonyucsdy/article/details/112682566
[9] 基于opencv的图像拼接 https://blog.csdn.net/qq_30111427/article/details/121127233