1 引言

多通道光谱成像设备在采集数据过程中,需要更换滤光片来调整成像通道参数。更换滤光片会造成多光谱相机位置的偏移,导致成像的多通道光谱图像之间产生几何畸变和偏移,并对后续的光谱图像像元分析产生影响。因此,必须对多通道光谱图像进行配准预处理校正。

配准方法主要有基于区域的配准和基于特征的配准^[1]。由于基于特征的配准方法计算量相对较小、稳健性强,对图像偏移和旋转等变化有较好的适应能力,已逐渐成为图像配准的主流方向^[2-3]。图像配准的关键在于寻找一种较好的特征点提取算法,目前常见的特征点提取方法有角点(Harris)特征、最小吸收核同值区(SUSAN)特征、尺度不变特征变换(SIFT)、快速稳健特征(SURF)等提取方法^[4-7]。Bay等^[8]提出的SURF相对于Harris与SUSAN特征具有尺度、旋转不变特性,在SIFT的基础上简化了图像金字塔分解的操作,对图像的卷积作近似处理,在特征点定位环节引入积分图的概念,使求解Hessian矩阵的计算量大大降低,运算时间低于基于SIFT特征方法的1/3^[9]。

多通道光谱相机成像为灰度图像,各通道反映不同波段的光谱特性,图像灰度值变化比较明显。相比于SIFT方法,SURF方法能够更好地反映光谱图像的特性,在亮度变化较大的情况下可以较好地提取特征点^[10]。但是使用SURF方法配准后的多通道光谱图像中存在边界区域像素值为零的无效区域,该区域会给后续基于多通道图像的光谱反射率重建和高精度复制带来计算误差。针对上述问题,本文提出了一种基于SURF融合最大子矩阵(MSM)的多通道光谱图像配准方法。利用SURF对多通道光谱图像进行初步配准,再提出最大子矩阵算法检测图像有效区域的最大内接矩形(MER),利用最大内接矩形分割初步配准结果中的有效信息区域,最大化地保留了有效区域。

2 特征提取与透视变换

2.1 特征提取

SURF算法对SIFT算法进行了改进,尺度不变性相对于Harris较好,时间复杂度相对SIFT较低,并且对图像亮度变化稳健性较强。SURF算法包括两个主要部分:特征点检测定位和生成特征点描述算子。

1) 特征点检测与定位。Hessian矩阵H(x,σ)是SURF算法的核心,是一个多元函数的二阶偏导数构成的方阵,描述了函数的局部曲率。SURF检测特征点之前,首先要计算积分图像I,利用图像I的Hessian矩阵完成特征点的选取与尺度变换,图像I中任意一个像素点p(x,y)处的值为原图像左上角到任意点p相应矩形区域的灰度值的总和,其数学公式为

I(p)=∑i=1i≤x∑j=1j≤yI(i,j),(1)

式中:i为行变量;j为列变量。在尺度σ上,点p的Hessian矩阵H(x,σ)定义为

H(x,σ)=Lxx(x,σ)Lxy(x,σ)Lxy(x,σ)Lyy(x,σ),(2)

式中:x为特征点坐标;σ为尺度;L_xx(x,σ)、L_xy(x,σ)和L_yy(x,σ)为图像I在点p与高斯二阶偏导数∂²g(σ)/∂x²的卷积,g(σ)为

g(σ)=12πσ2exp[-(x2+y2)2σ2]。(3)

为简化计算,SURF算法以盒子滤波器近似替代二阶高斯滤波器。在保证性能的同时,又提高了速度。以9 pixel×9 pixel盒子滤波器为例,取尺度σ=1.2,则高斯二阶偏导可近似如图1所示。

假设图像与上述盒子滤波器卷积后得到的Hessian矩阵中参数分别为L_xx、L_xy和L_yy,则(2)式中的H_essian矩阵行列式Det(H_essian)可近似表示为

Det(Hessian)=LxxLyy-(wLxy)2,(4)

式中:w为盒子滤波器的权重系数,当σ=1.2时,w可近似为0.9。建立尺度空间,在尺度空间中一个3 pixel×3 pixel×3 pixel的邻域内取每个点的Det(H_essian)值与当前尺度相邻位置和相邻尺度周围26个邻域的Det(H_essian)值进行比较,得到局部极大值点,通过插值计算可以求得在连续空间中最终的特征点。

2) 生成特征点描述算子。为了保证特征点的旋转不变性,需要为特征点赋予一个主方向,这是生成特征点描述算子的第一步。方法是在以该特征点为中心、6倍尺度为半径的圆形区域内,计算各点的x和y方向上Haar小波响应。然后对这些响应赋予高斯权重,将每60°的扇形区域w内的水平响应d_x和垂直响应d_y相加,得到局部方向矢量N(m_w,θ_w),

mw=∑wdx+∑wdyθw=arctan∑wdx/∑wdy,(5)

图 1. 盒子滤波器。(a) x方向;(b) y方向;(c) xy方向

Fig. 1. Box filter. (a) x-scale; (b) y-scale; (c) xy-scale

下载图片 查看所有图片

式中:m_w为w扇形区域内所有点的水平、垂直Haar小波特征总和;θ_w为w扇形区域内Haar小波特征局部方向角。比较所有局部方向矢量N,将最长矢量θ作为该特征点的主方向,

θ=β|β=max(N)。(6)

生成特征点描述算子的第二步是为特征点建立特征描述符。在获得特征点主方向后,以该主方向为x轴,在特征点周围的邻域内选取20 pixel×20 pixel的矩形区域,将该矩形区域等分为4 pixel×4 pixel的16个子域,再利用Haar小波模板获得沿主方向和垂直于主方向的高斯加权d_x和d_y,分别统计出 ∑dx、∑dx、∑dy和∑dy,形成特征矢量V=∑dx,∑dy,∑dx,∑dy,每个特征矢量是4维,共16个子域,所以最后得到一个64维的特征点描述算子。

2.2 特征点匹配

特征描述算子生成后,利用原图像与目标图像的特征描述算子进行相似度匹配。由于对比度相同的特征点对应的Hessian矩阵迹为同号,不同对比度对应的为异号,因此首先使用Hessian矩阵迹进行特征点的初步匹配,选择同号特征对,然后采用欧式距离匹配法对选择的同号特征对进行相似度判断,找到目标图像与源图像中待匹配特征点的一一对应关系。假设(x₁,x₂,…,x_N),(x'₁,x'₂,…,x'_N)为一对待匹配特征向量,则两向量之间的欧式距离为

D=∑k=1N(xk-x'k)2,(7)

式中:k为待匹配特征向量的元素序号;N为总的匹配特征点对数量。

计算目标图像与源图像上所有待匹配特征点的欧式距离,选出最小欧式距离D_near和次最小欧式距离D_{sub_near},求出D_near/D_{sub_near}的比值η,将η与设定阈值D_th进行比较,若η<D_th,则证明目标图像上该特征点与源图像上D_near对应点匹配;否则不匹配。实验表明,当D_th在0.6~0.9之间时,可以取得较好的匹配效果^[11]。

2.3 特征点筛选

通过初步匹配得到的特征点中还存在大量错误的匹配点,为了保证透视变换模型的计算精度,应尽量消除错误匹配点,保证图像配准的质量和效果,选用随机采样一致性算法(RANSAC)对特征点进行筛选。RANSAC是目前广泛采用的一种剔除错误匹配点的算法,由Fischler和Bolles于1981年最先提出^[12]。使用RANSAC算法的基本流程如下:

1) 随机从初步匹配集P中抽出4组特征点对(此4个特征点对不共线),计算变换矩阵H,记为模型M;

2) 计算P中所有数据与模型M的投影误差,若小于阈值t,加入内点集Z,并记录该模型下的统计误差;

3) 重复上述步骤,当计算得到新的模型时,比较其统计误差e与e_min大小,如果e更小,则更新模型M和e_min;

4) 输出最大内点集Z的最优模型M。

取4组特征点对是由于本文选用几何变换模型为透视变换模型,模型中有8个未知数,至少需要8组线性方程求解,一组特征点对可以列出两个方程,因此选择4组特征点对。将阈值设置为t=0.7时,可以获得较好的匹配正确率^[13-14]。

2.4 透视变换

根据匹配点转换源图像与目标图像间的坐标关系,即两幅图像之间的透视变换矩阵H可表示为

H=h0h1h2h3h4h5h6h71,(8)

式中:h₀,h₁,…,h₇为矩阵H的元素,若p(u,v),q(x,y)为匹配的特征点对,则投影变换公式为

(x,y,w)T=H×(u,v,1)T,(9)

式中:(u,v)为原始图像坐标;(x/w,y/w)为变换之后的目标图像坐标,由1.3节得到的最优模型M可得H中的各参数h₀,h₁,…,h₇,再利用透视变换完成原图像到目标图像的配准。

3 多通道光谱图像配准方法

3.1 基于SURF算法的光谱图像配准

使用SURF提取与透视变换对多通道光谱图像进行配准,先任取一个通道的光谱图像作为目标图像,依次计算其他通道光谱图像到目标图像的变换矩阵H,经过透视变换,可以得到配准后的多通道光谱图像。由(9)式可知,透视变换配准公式为

x'=xw=h0u+h1v+h2h6u+h7v+1y'=yw=h3u+h4v+h5h6u+h7v+1,(10)

式中:(x',y')为原图像上的像素点(u,v)通过(10)式转换到目标图像上的坐标。

图2(a)中,Target为目标图像,Channel 1、Channel 2和Channel 3为待配准图像,使用SURF对3个通道的光谱图像配准后的效果如图2(b)所示。

图 2. 多通道光谱原图和SURF配准后的图像。(a)目标图像与待配准多通道光谱图像;(b)目标图像与SURF配准后的多通道光谱图像

Fig. 2. Multi-channel spectral original image and image after SURF registration. (a) Target image and multi-channel spectral image to be registered; (b) target image and multi-channel spectral image after SURF registration

下载图片 查看所有图片

从图2可以看出,由于透视变换,SURF配准后各通道图像的边缘都存在不同程度的黑边,对所有通道进行叠加,得到的光谱立方体图像若具有不规则边缘,则会对后期的像元分析、光谱重建以及颜色再现造成误差,故需要对不规则边缘进行截取处理。

3.2 最大内接矩形检测

处理不规则边界主要是为了解决目标对象的最大内接矩形^[15]问题。在获取目标物体最小外接矩形后,在最小外接矩形内查找该物体的最大内接矩形,获取最大内接矩形面积。目前,常用的获取最大内接矩形的方法有遍历法和中心扩散法^[15]。

1) 遍历法。给出一个坐标点,计算以该点为右上角并且所有点都位于目标物体内的所有矩形面积,比较矩形面积,得到以该点为右上角的最大矩形面积。对目标物体内每一点均调用上述方法,即可求得目标物体最大内接矩形。

2) 中心扩散法。计算最小外接矩形中心点O( x^, y^),即

x^=(pl+pr)2y^=(pu+pd)2,(11)

式中:p_l、p_u分别为最小外接矩形的左上角顶点横坐标与纵坐标;p_r、p_d分别为最小外接矩形的右下角顶点横坐标与纵坐标。设定4个变量p'_u、p'_d、p'_l、p'_r,并赋初值

p'l=x^-1,p'l≤plp'r=x^+1,p'r≤prp'u=y^-1,p'u≤pup'd=y^+1,p'd≤pd。(12)

在区间(p'_l,p'_r)上遍历p'_u行,判断像素是否为0,若当前行像素均不为0,将p'_u行向上移动,即p'_u-1,否则停止移动;同样在区间(p'_l,p'_r)上遍历p'_d,若当前行像素值均不为0,则p'_d行下移,即p'_d+1,否则停止移动;在区间(p'_u,p'_d)上遍历p'_l,若当前列像素值均不为0,将p'_l列左移,即p'_l-1,否则停止移动;同样在区间(p'_u,p'_d)上遍历p'_r,若当前列像素值均不为0,则p'_r列右移,即p'_r+1,否则停止移动。当p'_l、p'_r、p'_u、p'_d均停止移动时,计算其围成的面积。

在实际使用中,遍历法具有很强的稳健性,但是算法复杂度较高,不适合大面积目标物体检测,中心扩散法算法复杂度较低,但是对不规则的物体边缘十分敏感,且检测准确性较低。

3.3 基于MSM-SURF的多通道光谱图像配准

针对SURF初步配准和现有最大内接矩形检测方法存在的问题,提出了将SURF提取与MSM相融合的MSM-SURF算法,用于实现多通道光谱图像配准。MSM-SURF算法框图如图3所示。

图 3. MSM -SURF算法框图

Fig. 3. Block diagram of MSM-SURF algorithm

下载图片 查看所有图片

将SURF初步配准后的n个通道光谱图像M_g(g=1,2,…,n)叠加,得到矩阵I,叠加过程为

Iab=0,min({k|mk,ab})=0max({k|mk,ab}),min({k|mk,ab})≠0,(13)

式中:a、b分别为矩阵M_g的行和列;I_ab为矩阵I的第a行b列元素;m_k,ab为初步配准后第k个通道的光谱图像矩阵M_k在a、b处的像素值。580,620,715 nm通道光谱图像矩阵M_g经过(13)式叠加后得到的I矩阵效果图如图4所示。

图 4. 叠加效果图

Fig. 4. Effect image of superposition

下载图片 查看所有图片

从图4可知,SURF初步配准后得到的叠加矩阵将各个通道图像边缘的无效区域都进行了累积叠加,图像立方体的边缘是不规则的。对叠加矩阵I进行二值化,完成最大内接矩形检测。

针对遍历方法和中心扩散方法求解中的问题,本文提出了MSM方法进行最大内接矩形检测,采用最大子矩阵方法的思路如下。

在给定矩阵中,假设矩阵元素只包含0与1,找出其中不包含任何0元素的所有子矩阵,计算所有子矩阵面积,其中面积最大的子矩阵即为最大子矩阵。算法框图如图5所示。

主要包含以下步骤:

1) 将矩阵I中值大于阈值Q的点的值设定为1,其余点值均为0,对目标矩阵I进行二值化。由于多通道光谱图像是灰度图像,经透视变换后边缘无效区域像素值均为0,因此,本文中将阈值Q设置为0。

Iij=0,Iij=QIij,Iij>Q。(14)

2) 创建矩阵N,大小与矩阵I相同;对N初始化,令N_1j=I_1j,当i>1时,若I_ij≠0,则N_ij=I_ij+I_(i-1)j,否则N_ij=0,即

Nij=I1j,i=10,i>1,Iij=0Iij+I(i-1)j,i>1,Iij≠0。(15)

3) 创建数组A_arr[r_ow],r_ow为矩阵N的列数,A_arr存放数组原始数据,遍历数组;

(1) 当i=0时,此时栈为空,将数组索引i=0入栈;当i≠0时,先获取栈顶元素p_op,若A_arr[i]≥A_arr[p_op],则将i入栈;

(2) 若A_arr[i]<A_arr[p_op],弹出栈顶元素p_op,令h=A_arr[p_op],w=p_op-p'_next-1,其中p'_next为p_op出栈后新的栈顶元素,令S_k=h×w,重复步骤(1),直到数组遍历完成;

(3) 比较S_k大小,求出其中最大值S_max={k=1,2,…,n|max(S_k)};

4) 按行遍历矩阵N,令A_arr[r_ow]=N_i,重复执行3),得到每一行S_max,第i行中S_max围成区域即为矩阵I第1~i行中最大子矩阵S'_max,令P={i=1,2,…,n|max(S'_max)},则P即为矩阵I最大子矩阵。

计算得到的最大子矩阵如图6所示。

使用最大子矩阵对SURF初步配准后的每个通道的光谱图像M_g进行截取,得到每个通道截取后的配准结果,MSM-SURF配准后的多通道光谱图像如图7所示。

图 5. 最大子矩阵算法流程图

Fig. 5. Flow chart of maximum submatrix algorithm

下载图片 查看所有图片

图 6. 最大子矩阵

Fig. 6. Maximum submatrix

下载图片 查看所有图片

图 7. MSM-SURF配准后的多通道光谱图像。(a) 580 nm;(b) 620 nm;(c) 715 nm

Fig. 7. Multichannel spectral images after MSM-SURF registration. (a) 580 nm; (b) 620 nm; (c) 715 nm

下载图片 查看所有图片

综上所述,MSM-SURF算法首先使用SURF提取与投影变换对多通道光谱图像进行初配准,再针对初步配准后的图像边缘像素值为0的无效区域,使用最大子矩阵MSM进行分割,在截取每个通道光谱图像的无效区域的同时,能够最大化地保留有效区域的信息。

4 实验结果及分析

4.1 实验数据

本实验方案采用两组壁画共6个通道的光谱图像作为实验数据。采集系统包括海洋光学公司的SpectroCam VIS型号电荷耦合器件(CCD)多光谱相机(500万像素)、6个窄带干涉滤光片和CIE(Commission Internationale de L'Eclairage)标准A照明光源。实验采用Matlab仿真,平台为Intel i5-2400 3.1 GHz 四核处理器、4G内存。采集得到每组壁画在6个不同波段通道下的光谱图像如图8所示,图中每组均以580 nm通道的壁画为配准目标图像。图8中,每组壁画中的各通道光谱图像都存在不同程度的位移几何畸变,需要进行配准校正。

图 8. 壁画的6通道多光谱图像。(a)壁画A;(b)壁画B

Fig. 8. 6-channel multi-spectral images of murals. (a) Mural A; (b) mural B

下载图片 查看所有图片

4.2 评价标准

为了测试算法的配准性能,从配准精度、特征点匹配准确率、有效区域像素点占比率以及运算时间效率4个方面对两种算法进行客观评价^[15-16]。

1) 配准精度R_MSE。由源图像与目标图像所有匹配特征点对之间坐标的方均根误差R_MSE确定,即

RMSE=1n∑τ=1n(x'τ,y'τ)-f(xτ,yτ)2,(16)

式中:(x'_τ,y'_τ)为目标图像上第τ个特征点坐标;(x_τ,y_τ)为待配准图像上第τ个特征点坐标;f表示坐标转换关系;n为经筛选后的匹配特征点数目。

2) 特征点匹配准确率K。正确匹配特征点对数目在所有匹配特征点对数目中的比例 ,即

K=RN×100%,(17)

式中:R为最终正确匹配特征点对数量。

3) 有效区域像素点占比率R_ate。最大内接矩形提取的区域所包含像素点数目与整幅图像像素点数目的比例,即

Rate=cs×100%,(18)

式中:c为最大内接矩形包含像素点数目;s为整幅图像像素点数目。

4.3 实验结果分析

将SURF与不同的内接矩形检测方法融合。对壁画A和B的数据分别采用遍历SURF(T-SURF)、中心扩散SURF(CD-SURF)和MSM-SURF算法进行配准对比。3种方法对初步配准后叠加矩阵的检测结果如图9所示。

图 9. 3种内接矩形检测方法的检测结果。(a) T-SURF;(b) CD-SURF;(c) MSM-SURF

Fig. 9. Test results with three inscribed rectangle detection methods. (a) T-SURF; (b) CD-SURF; (c) MSM-SURF

下载图片 查看所有图片

从图9可以看出,相比于T-SURF和CD-SURF算法,MSM-SURF算法对叠加矩阵的检测结果的信息保留最大。

对上述检测结果进行量化分析,得到最终配准后的配准精度R_MSE、特征点匹配准确率K、运行时间Time和有效区域像素点占比率R_ate的结果,如表1所示。

从表1可以看出,由于3种算法均使用SURF,在配准精度与特征点匹配正确率上几乎相同;MSM-SURF的有效占比率略高于T-SURF,CD-SURF的有效占比率只有T-SURF与MSM-SURF算法的60%;CD-SURF和MSM-SURF在运行时间上相近,但T-SURF算法的复杂度远高于CD-SURF和MSM-SURF算法,T-SURF的时间效率是CD-SURF与MSM-SURF的3倍左右。综合上述分析,MSM-SURF算法不仅具有较高的检测有效占比率,并且算法时间复杂度较低,更适合于多通道光谱图像的配准。

5 结论

针对多通道光谱图像的特点,提出了一种SURF与最大子矩阵融合的MSM-SURF算法来改进现有匹配和内接矩形检测方法的缺陷。并对采集的实际壁画的多通道光谱图像数据进行了实验。实验结果表明,MSM-SURF算法能有效地将SURF与最大子矩阵检测相结合,解决了SURF配准出现无效边缘区域的问题,同时能够将有效区域信息最大化地保留,对于多通道光谱图像的配准具有较好的实际应用意义。