CN104599258A - 一种基于各向异性特征描述符的图像拼接方法 - Google Patents

Abstract

本发明公开了一种基于各向异性特征描述符的图像拼接方法，适用于拼接多幅具有一定重叠关系的图像。本发明的方法包括如下步骤：A、检测参考图像和待配准图像的特征点，求取特征点的主方向；B、采用各向异性的点对点抽样模型，构成多组二值测试，得到特征描述符；C、采用最近邻方法对特征描述符进行匹配，利用PROSAC方法求取单应性矩阵；D、通过求解误差函数获得光照增益补偿矩阵，利用多波段融合获得过渡自然的全景图。本发明既能够准确地配准不同视角、视点拍摄的图像，获得清晰自然的宽视角的场景图像，又具有简单的复杂度与较快的运行速度，从而为监控***或遥感***提供很好的应用价值。

Description

一种基于各向异性特征描述符的图像拼接方法

技术领域

本发明属于图像信息处理的技术领域，本发明是一种基于各向异性特征描述符的图像拼接方法，适用于拼接在不同视角或视点拍摄的多幅具有一定重叠关系的图像。

背景技术

图像拼接已经成为一个日益流行的研究领域，在户外监控***、医学图像分析、遥感图像处理等领域得到了广泛的应用。当使用普通相机来获取宽视野的场景图像时，必须通过调节焦距才可以获得完整的场景，但这样会以损失场景图像的分辨率为代价。而利用价钱昂贵、操作复杂的广角镜头和扫描式相机可以解决视角不足这一问题，但是广角镜头的边缘很容易产生扭曲变形。图像拼接技术是将几路普通图像或视频图像进行无缝拼接，得到更宽视角的场景图像，从而可把普通相机拍摄的多幅不同视角的图片拼接为全景图。图像拼接的目的就是提供一种自动匹配的方法，将具有一定重叠区域的多幅图片合成一副广视角图片，来扩大视区的范围，具有重要的现实意义。

图像获取、图像预处理、图像配准和图像融合这四步是图像拼接的整体流程。其中图像配准是图像拼接的核心部分，其目标是找出几幅重叠图像之间的相对变换关系，直接影响整个***的成功率和运行速度。

基于特征点匹配的图像拼接算法是目前图像拼接算法的研究热门。Harris算法是最早提出的特征点检测模型，该特征具有旋转不变性，对光照和噪声也具有很好的鲁棒性。David G low于2004年提出的尺度不变变换特征的算法对平移、旋转和缩放等操作均有较好的鲁棒性，同时对光照和噪声具有很高的抗干扰性能。

论文名：Automatic Panoramic Image Stitching using Invariant Features，期刊： International Journal of Computer Vision(IJCV)，年份：2007年。Brown等人提出了一种基于SIFT特征点匹配的全景图拼接方法，该算法利用Low等人提出的尺度不变特征，通过检测、描述和匹配SIFT特征点进行图像配准，然后利用增益误差函数估算出图像块的光照增益，对光圈或光线影响造成的亮度差问题进行补偿，最后利用多频段融合方法来消除拼接形成的缝合线，达到了较好的拼接效果。但由于SIFT算法运算量庞大，耗时长，导致该算法很难应用到实际的场合。

论文名：ORB:An Efficient Alternative To SIFT or SURF,会议：International Conference on Computer Vision(ICCV)，年份：2011年。Ethan等人采用了FAST作为特征点检测算子，提出了一种具有旋转不变性的二值特征描述符。特征点的主方向通过计算矩得到，在特征点附近随机选取若干点对，将点对点灰度值比较的二值测试组合成一个二进制串特征描述子。二值描述符之间的距离利用汉明距离计算。该特征点方法具有较快的速度。但是由于FAST特征点鲁棒性较差，当应用到图像拼接中时，受非重叠区域特征点以及其他外点的影响，匹配正确率急剧下降。而且当图像有变形扭曲时，各向同性的二值描述符匹配性能也受到很大的影响。

针对上述背景内容，研究一种快且有效的特征点配准方法，从而应用到全景图拼接中，具有重要的意义。

发明内容

本发明的目的是为了克服现有图像拼接算法的不足，提供一种基于各向异性特征描述符的图像拼接方法，既能够准确的提取和匹配特征，精确的实现图像配准，从而清晰自然地合并多幅具有一定重叠区域的图像，扩大图像视区的范围，又具有简单的复杂度与较快的运行速度，从而为监控***或遥感***提供很好的应用价值。

本发明提供的技术方案包括如下步骤：

A、检测图像中的特征点，求取特征点的主方向和Hessian阵；

B、采用各向异性的点对点抽样模型，抽取点对构成多组二值测试，最终组成特征描述符；

C、采用最近邻方法对特征描述符进行匹配，利用PROSAC方法去除错误匹配对，求得图像之间单应性变换矩阵；

D、利用重叠区域光照强度误差函数获得光照补偿增益矩阵；消除图片融合的缝合线，利用多波段融合获得过渡自然的全景图；

所述步骤A为：

A1、采用Hessian矩阵行列式近似值图像，构造高斯金字塔尺度空间。图像中某个像素点的Hessian矩阵定义为

$H(x,\mathrm{\σ})=\left[\begin{array}{cc}{L}_{\mathrm{xx}}(\hat{x},\mathrm{\σ})& L_{\mathrm{xy}}(\hat{x},\mathrm{\σ})\\ L_{\mathrm{xy}}(\hat{x},\mathrm{\σ})& L_{\mathrm{yy}}(\hat{x},\mathrm{\σ})\end{array}\right]---\left(1\right)$

其中，L_xx、L_xy、L_yy是选用二阶标准高斯函数作为滤波器，通过特定核间的卷积计算二阶偏导数。Bay等人提出用方框滤波近似代替二阶高斯滤波，用积分图像来加速卷积计算，并通过改变方框的大小形成不同尺度的图像金字塔。每个像素的Hessian矩阵行列式近似公式如下：

Δ(H)＝D_xxD_yy-(0.9D_xy)²             (2) 

其中D_xx、D_yy，D_xy分别为方框滤波模板同图像卷积后的二阶偏导近似值。

A2、特征点的搜索是通过在同一组内各相邻层之间比较完成的，在3×3×3的立体领域内进行非极大值抑制，然后在尺度空间中进行亚像素插值运算，得到精确的位置坐标；

A3、为保证特征点的旋转不变性，需计算特征点的主方向，以特征点为中心，统计邻域内点的Haar小波响应。

所述步骤B为：

B1、确定每个特征点的抽样模型，这里采用各向异性的点对点抽样模型。经典的二值描述符ORB、FREAK的抽样模型的定义如下：

Λ_i＝R_θ·Φ_i,               (3) 

Φ_i＝[r_i cosθ_i r_i sinθ_i]^T

其中R_θ为特征点的主方向，为保证特征点的旋转不变性。r_i，θ_i分别为随机抽样模型Φ_i中任一点i的半径和角度。当图像出现扭曲变形时，抽样模型必须跟随做出修正。Schmid等人在论文《Scale and Affine Invariant Interest Point Detectors》中证明了特征点附近的仿射模型可以通过与Hessian阵的平方根矩阵相乘得到纠正。

因此，某一特征点的随机抽样模型将被修正为：

Λ_i′＝H^-1/2·R_θ·Φ_i              (4) 

其中H^-1/2为任一特征点的Hessian阵的平方根矩阵的逆。

B2、特征点的二值描述符是由特征点附近的多组点对对比构成，即通过比较两点像素的强度形成一位二值测试，最终由多个二值测试组成二值描述符F： 其中N是描述符的长度，p_i,p_j表示一个点对，Τ(Λ′；p_i,p_j)是一位二值测试，其表达形式如下所示

$T({\mathrm{\&Lambda;}}^{\&prime;};p_i,p_j)=\left\{\begin{array}{cc}1& I({\mathrm{\&Lambda;}}^{\&prime;},p_i)<I({\mathrm{\&Lambda;}}^{\&prime;},p_j)\\ 0& \mathrm{otherwise}\end{array}\right.---\left(5\right)$

其中I(Λ′,p_i)和I(Λ′,p_j)为各向异性随机抽样模型Λ′上的随机取样点对p_i和p_j的强度。

所述步骤C为：

C1、由步骤B获得参考图像I₁与待配准图像I₂的特征描述符后，进行特性匹配，对于二值描述符，Hamming距离是一个理想的相似度测量工具。采用最近邻匹配法，对I₁中任一个特征点n_1i，I₁中与之汉明距离最小的两个特征点分别为n_2j、n_2j′，(对应距离为d_ij、d_ij′)，如果d_ij≤a^*d_ij′，则认为是匹配的点对。

C2、估算图像间的变换关系，采用的变换模型为单应性矩阵，它符合平面目标 在不同视点或视角成像的变换关系。变换关系如下：

$\left(\begin{array}{c}{x_i}^{\&prime;}\\ {y_i}^{\&prime;}\\ 1\end{array}\right)=\left(\begin{array}{ccc}{m}_{11}& m_{12}& m_{13}\\ m_{21}& m_{22}& m_{23}\\ m_{31}& m_{32}& 1\end{array}\right)\left(\begin{array}{c}{x}_i\\ y_i\\ 1\end{array}\right)---\left(6\right)$

其中(x_i,y_i)，(x_i′,y_i′)分别为图像和待配准图像上匹配点对的坐标。这里采用PROSAC(progressive sample consensus)算法来去除误匹配点，求取变换矩阵。

所述步骤D为：

D1、将图像进行分块，并且求出任意两块图像的重叠区域的平均光强，其数学表达式为：

$\stackrel{\&OverBar;}{I_i}=\frac{\underset{\mathrm{if}\left(\mathrm{inter}\sec t\right)}{\mathrm{\&Sigma;}}\sqrt{\sqrt{R}+\sqrt{G}+\sqrt{B}}}{N(i,j)}---\left(7\right)$

式中，N(i,j)为图像块i与图像块j相交区域的像素总数目，R、G、B为图像块i在相交区域任一点R、G、B三通道的像素值。建立误差函数如下：

$e=\frac{1}{2}\underset{i}{\mathrm{\&Sigma;}}\underset{j}{\mathrm{\&Sigma;}}(\underset{\mathrm{if}\left(\mathrm{inter}\sec t\right)}{\mathrm{\&Sigma;}}{(g_i{\stackrel{\&OverBar;}{I}}_i-g_j\stackrel{\&OverBar;}{I_j})}^2/{\mathrm{\&delta;}}_N^2+{(1-g_i)}^2/{\mathrm{\&delta;}}_g^2)---\left(8\right)$

式中，g_i为某一图像块的光照增益系数。δ_N和δ_g分别为亮度和增益系数的标准差。

D2、根据步骤D1，将对应区域乘以g_i进行增益补偿。然后设定层数N_bands，对每幅图像构建金字塔。具体流程如下：首先调整图像的宽和高，使之可被除尽，这样才能下采样N_bands次，然后再由最底部的图像上采样N_bands次。把对应层上采样和下采样的差放入金子塔相应层中。最后将各层金字塔叠加，得到完整的全景图。

本发明的有益效果：

(1)采用各向异性的特征描述符，对特征点处的点对点抽样模型进行合理正确的修正，在图片有扭曲变形的恶劣情况下，依然具备较好的描述和匹配性能。与FREAK、ORB等算法相比，本发明可以计算出准确的单应性变换矩阵。另外由于 各向异性的特征描述符是二值描述符，具备二值描述符计算快速、匹配算法复杂度低的优点。与SIFT、SURF等算法相比，本发明具有更快的速度。在实时性***中具有很好的实际应用价值。

(2)增加了光照补偿和多波段融合的方法，在下采样的情况下计算增益矩阵，该方法不仅计算快速而且非常有效，与加权平均的图像融合方法相比，图像交叉区域过渡自然，保留了更多的细节。

附图说明

图1为一种基于各向异性特征描述符的拼接流程示意图。

图2为计算特征点描述符的流程示意图。

图3为计算特征点描述符的点对点抽样模型。其中(a)和(b)分别为FREAK抽样模型和本发明采用的各向异性采样模型。

图4为PROSAC方法计算单应性矩阵的流程示意图。

图5为两幅待拼接图像及三种拼接方法的处理结果。

其中，(a)和(b)为参考图像和待配准图像；(c)和(d)检测的特征点；(e)为加权融合方法的拼接结果；(f)为本发明对(a)和(b)的拼接结果。

具体实施方式

下面结合具体实施例和附图详细阐述本发明。

A、检测参考图像和待配准图像(见附图5(a)、5(b)所示)的特征点，并计算特征点的主方向。

A1、对参考图像和待配准图像先做灰度化处理，采用Hessian矩阵行列式的近似值图像，构造高斯金子塔空间。具体操作见公式(1)(2)。

A2、在同一组内各相邻层3×3×3的空间邻域内搜索特征点，在尺度空间中进行亚像插值运算，得到特征点的准确位置坐标。

A3、以特征点为中心、6δ为半径做圆，计算圆形区域上任一点在x和y方向 上大小为4δ的Haar小波响应，其中δ为特征点所在尺度空间的尺度。最后以60°范围作为一个区域，遍历一周得到6个扇形区域，每个区域的响应相加得到新的矢量，把具有最大模值的矢量方向作为该特征点的主方向。由此求得的两幅图的特征点(部分特征点见附图5(c)(d))。

B、采用各向异性的抽样模型抽取点对构成二值描述符，方法的流程示意图见附图2所示：

B1、这里采用FREAK方法的7层Retina抽样模型(见图3(a))，按照公式(4)修正Retina抽样模型得到各向异性的抽样模型(示例见图3(b))。

B2、按照步骤B1求取的各向异性的点对点抽样模型，形成N组二值测试，最终组成二值描述符F。这里N取值为128，即二值描述符F的维数为128。

C、利用最近邻方法对特征描述符进行匹配，然后采用PROSAC方法去除误匹配点对，求出单应性矩阵。

C1、按照公式d_ij≤a^*d_ij′，比较I₁中任一特征描述符与I₁中距离最小和次最小的汉明距离，实验中取a^*＝0.7，满足则认为是匹配的点对。

C2、采用PROSAC方法去除误匹配对，方法的流程示意图见附图4所示。由此求得的图1(a)与图1(b)的单应性矩阵为：

$\mathrm{Homography}=\left(\begin{array}{ccc}0.666871& -0.012074& 636.078301\\ -0.082783& 0.949450& -4.201992\\ -0.000340& 0.000024& 1.000000\end{array}\right)$

D、采用光照增益补偿和多波段融合获得清晰而且过渡自然的全景图。

D1、将图像进行分块，在具体实施中，为了对算法加速，首先对图像下采样，将图像下采样到总像素面积S的尺度上，由大量实验经验值选取S为10⁵。将图像分成32×32的图像块，按照公式(7)计算平均光强，最后按照公式(8)求解误差函数。具体实施中，分别选取δ_N和δ_g为10和0.1。

D2、利用步骤D1求得的增益矩阵对图像进行增益补偿，然后对每幅图像构建 金字塔，设置金子塔层数为5。调整图像的宽和高，使其可被32整除。将图像下采样5次，然后由最底部的图像上采样5次，将对应层上采样和下采样的差放入金字塔中，将各层金字塔叠加，得到最终的全景图。实施步骤D1和D2前后的对比效果见附图5(e)和5(f)。

经过上述步骤，图5(f)为本发明对不同视角图像5(a)和5(b)的拼接结果。

上述实施例的实施平台为Windows 7(64位)操作***、处理器主频为3.2GHz、***内存4G的PC上，Microsoft Visual C++2010软件。图5(f)使用sift特征点配准与使用本发明各向异性特征描述符配准的两种方法对图5(a)和图5(b)的拼接结果。

对于图像大小为1000×562的图5(a)和5(b)，本发明的图像配准处理时间为0.279s，而sift与surf的图像配准处理时间分别为2.83s和0.86s。

本发明若移植于FPGA硬件平台，采取并行运算，则可进一步加速。

Claims (1)

1.一种基于各向异性特征描述符的图像拼接方法，其特征在于以下步骤：

A、检测参考图像和待配准图像的特征点，计算特征点的主方向和Hessian阵；

(1)对图像进行下采样初始化处理，将图像下采样到总像素面积0.6×10⁶的尺度上；然后对图像做灰度化处理，采用Hessian矩阵行列式的近似值图像构造高斯金字塔空间，表达式为：Δ(H)＝D_xxD_yy-(0.9D_xy)²；用积分图像来加速卷积计算；

(2)在同一组内各相邻层3×3×3的立体领域内搜索候选特征点；并在尺度空间上进行亚像素插值得到精确的位置坐标；

(3)以特征点为中心，计算半径6s邻域内的点的Haar小波响应；以60度范围作为一个区域，遍历一周得到6个扇形区域，每个区域的响应相加得到新的矢量，把具有最大模值的矢量方向作为该特征点的主方向；

B、采用各向异性的点对点抽样模型构造二值测试对，组成多维特征描述符；具体步骤如下：

(1)采用FREAK的7层Retina抽样模型，计算每个特征点i的Hessian矩阵和主方向，修正FREAK抽样模型，其表达式为：Λ_i′＝H^-1/2·R_θ·Φ_i，式中R_θ为特征点的主方向，H^-1/2为任一特征点的Hessian阵的平方根矩阵的逆；通过上述公式，将FREAK模型Λ_i＝R_θ·Φ_i修正为各向异性的抽样模型Λ_i′，在图像出现扭曲变形的恶劣情况下，该描述符依然具备很好的描述性能；

(2)在抽样模型Λ_i′上随机抽样一组点对，比较两点像素的强度值形成一位二值测试，其表达式为： $T({\mathrm{\&Lambda;}}^{\&prime;};p_i,p_j)=\left\{\begin{array}{cc}1& I({\mathrm{\&Lambda;}}^{\&prime;},p_i)<I({\mathrm{\&Lambda;}}^{\&prime;},p_j)\\ 0& \mathrm{otherwise}\end{array}\right.;$ 式中I(Λ′,p_i)和I(Λ′,p_j)为抽样模型Λ′上的随机取样点对p_i和p_j的强度；最终由512个二值测试组成二值描述符F： $F=\underset{i=1}{\overset{512}{\mathrm{\&Sigma;}}}{2}^{i-1}T({\mathrm{\&Lambda;}}^{\&prime;};p_i,p_j);$

C、采用最近邻方法对特征描述符进行匹配，采用PROSAC去除误匹配对，求取单应性矩阵；

(1)比较图像I₁中任一特征描述符与图像I₂中汉明距离最小和次最小的描述符，按照公式d_ij≤a^*d_ij′，实验中取a^*＝0.7，满足则认为是匹配的点对；

(2)采用PROSAC去除误匹配对，首先按照最小与次最小汉明距离比a^*将匹配数据排序，设置最大迭代次数和内外点的误差门限，在前n-1个数据中抽取m-1个数据和第n个数据组成样本计算单应性矩阵，内点个数大于设定的阈值则迭代结束，否则继续迭代，在前n个数据中抽取m-1个数据和第n+1个数据组成样本计算单应性矩阵和内点个数，直至满足内点个数大于设定阈值或者迭代次数超过阈值为止；

D、采用光照增益补偿和多波段融合获得清晰而且过渡自然的全景图；具体步骤如下：

(1)将图像下采样到总像素面积为0.1×10⁶的尺度上，将图像分成32×32的图像块，计算平均强度，其表达式为：式中N(i,j)为图像块i与图像块j相交区域的像素总数目；建立误差函数表达式为： $e=\frac{1}{2}\underset{i}{\mathrm{\&Sigma;}}\underset{j}{\mathrm{\&Sigma;}}(\underset{\mathrm{if}\left(\mathrm{inter}\sec t\right)}{\mathrm{\&Sigma;}}{(g_i{\stackrel{\&OverBar;}{I}}_i-g_j\stackrel{\&OverBar;}{I_j})}^2/{\mathrm{\&delta;}}_N^2+{(1-g_i)}^2/{\mathrm{\&delta;}}_g^2);$ 式中g_i为图像块i的光照增益系数；δ_N和δ_g分别为亮度和增益系数的标准差分别取10和0.1；

(2)求取的增益矩阵对参考图像和待配准图像进行增益补偿，然后进行多波段融合，具体流程如下，建立拉普拉斯金子塔，设定层数为5；首先调整图像的宽和高，使其可被32整除，向下采样5次，然后由最底部的图像向上采样5次，将对应层上下采样的差放入金字塔中；最后将5层金字塔叠加，得到最终的全景图。