赵 悟,段永璇,段会川2,张 睿,肖宪翠2,岳 媛
医学图像特征匹配[1-2]是当前数字图像处理方面的研究热点之一,如何从待匹配的医学图像中提取出含有图像特征的稳定特征点是其主要任务之一。研究表明,Harris角点检测算法对各种变换(如场景变换、旋转变化等)均能保持较好的适应性,因此在医学图像特征提取中被广泛应用。
Harris角点检测算法提取医学图像特征点的过程是通过构造自相关矩阵,再利用角点响应函数选择初始特征点,最后利用非极大值抑制的方法筛选出符合条件的特征点完成的,该过程运算简便且能获得一定数量的特征点。然而,在实际应用中发现传统的Harris算法检测到的特征点数量不足且图像配准精度不高,因此通过优化Harris算法过程提高特征点数量与图像配准精度是十分有意义的。本文提出了一种优化算法——GM-Harris算法,即利用群搜索优化算法(GSO算法)与互信息相结合的方式优化传统Harris算法的过程。实验表明,优化后的Harris算法不但可以获得较充足的特征点,而且能提高图像配准的精度。
Harris角点检测算法[3-8]是由CHris.Harris和Mike Stephens在Moravec算法的基础上于1988 年提出的,该算法是对Moravec算法的扩充和完善。Harris 算法提取图像特征点的基本思路如下。
公式(1)
第二步,利用高斯核函数G(x,y,σ)对图像进行高斯滤波,得到新的自相关矩阵M2。
第三步,利用角点响应函数R计算原图像上对应的每个像素点的响应值,即R值。角点响应函数R=Det(M2)-k×Tr2(M2),Det(M2)=λ1λ2,Tr(M2)=λ1+λ2,其中λ1、λ2为自相关矩阵M2的特征值,k为经验值。如果某点的角点响应值R大于设定的阈值,则该点就被选定为特征点。
第四步,选取局部的极值点。根据给定的阈值,采用非极大值抑制的方法对不符合条件的极值进行置零处理,从而确定最终的特征点。
医学图像之间的互信息[9-10]能有效反映两幅医学图像之间的相关性,因此可以作为图像相似性度量的方法。互信息来源于信息论中的信息熵,熵表达的是某一事物的复杂性或者是不确定性。一幅图像的熵反映的是该图像中像素灰度值的分布情况,灰度级别越高,熵就越大,图像的信息量就越丰富。图像的熵通常采用概率分布来描述,具体公式如下:
H(A)=-∑aPA(a)logPA(a)
公式(2)
H(B)=-∑bPB(b)logPB(b)
公式(3)
H(A,B)=-∑aPA,B(a,b)logPA,B(a,b)
公式(4)
公式(2)至公式(4)中,H(A)、H(B)为图像A和B的信息熵,H(A,B)为A和B的联合熵,a∈A、b∈B、PA(a)、PB(b)分别为图像A、B的概率分布,PA,B(a,b)为2幅图像的联合概率分布。互信息可以用信息熵表示,具体如下。
I(A,B)=H(A)+H(B)-H(A,B)
公式(5)
群搜索优化算法(GSO)[11-12]来源于自然界中动物的觅食行为,该算法在特定空间中寻找最优解的过程是根据动物在自然界中寻找食物的过程模拟的。GSO算法主要用来处理连续空间的最优值问题,其中的整体称为“种群”,种群中每个单独的个体称为“成员”。种群成员被划分为发现者、加入者、游荡者3类,其中发现者的任务是迭代查找“资源”;加入者接近于发现者,追随发现者并共享已获得的“资源”。为了保持种群的多样性,避免算法陷入局部最优,GSO 算法定义了新成员——游荡者,即在空间中作为游荡者向任意方向随机搜寻。
公式(6)
发现者方案:发现者从0°开始搜寻,在特定空间中依据公式(7)随机寻找与其相邻的右、前、左3点。
公式(7)
r1∈R1为符合正态分布(标准差为1、均值为0)的随机数,r2∈Rn-1为在(0,1)之间分布的随机数,lmax为最大移动距离。
如果公式(7)获得新位置的适应度值优于先前的位置,则发现者就移动到新位置,否则返回先前位置并调整角度继续搜寻。调整角度的公式如下:
φt+1=φt+r2amax
公式(8)
式中,amax为最大偏转角度。
若n次迭代后,发现者未寻找到更佳的位置,则返回到最初0°角的位置,即:
φt+n=φt
公式(9)
式中,n是一个常数。
加入者方案:在每次迭代中,各加入者按公式(10)以随机步长向发现者靠近,追随发现者参与搜寻。
公式(10)
式中,r3∈Rn为在(0,1)之间分布的随机数。
游荡者方案:在第t次迭代中,作为游荡者的第i个成员依据式(11)采用游荡方案,随机搜寻。
公式(11)
本文提出了GSO算法与互信息相结合的方式——GM-Harris算法,并将公式(5)作为GM-Harris算法的适应度函数,以此选取合适的特征点(角点)。GM-Harris算法的基本思路如下。
首先,提取标准图像的特征点,即种群初始化,记录每个成员(像素点)的位置(水平坐标x,垂直坐标y,角度a)以及初始的特征点数目。
其次,依据每个成员的位置参数计算相对应的浮动图像的有关数值。
第三,计算每个成员对应的标准图像与浮动图像的互信息值(适应度值)。
第四,依据GSO算法以及每个特征点(角点)的互信息值更新成员的局部极值和全局极值,再利用迭代公式,即公式(6)、公式(7)、公式(11)更新每个特征点的位置。
最后,判定是否满足GM-Harris算法的终止条件。如果未满足终止条件,则将当前每个特征点的位置作为新的种群参数,重新计算标准图像与浮动图像相应的特征点的适应度值,继续使用GSO算法进行迭代查找。
GM-Harris算法伪代码如下:
创建医学图像特征点提取事件序列
WHILE(迭代次数未达到最大值)//设定终止条件
{//参数设置
设定特征点(角点)种群大小为P,特征点(角点)数目为FP,最大迭代次数为NI;
用目标函数(角点响应函数)筛选初始的特征点数目,再保存当前最优的特征点数;
//迭代查找过程
记录标准图像成员的初始位置,再根据目标函数来获取最初的特征点数目;
再依据每个成员的位置参数计算相对应浮动图像的有关数值;
然后计算每个成员对应的标准图像与浮动图像的互信息值;
依据GSO算法以及每个特征点(角点)的互信息值,逐步更新特征点信息,以获取最新的特征点数;
对于依次获得的特征点数,需与其邻近的前一特征点数进行作差运算:
IF(|差值|>0)
{
记录当前的FP,并将其当作局部最优值进行记录;
}
ELSE
{
寻找到最优的FP,并保存下来,迭代查找结束。
}//END IF
}//END WHILE
根据上述结果,提取出最优的FP。
为定量分析传统Harris算法与GM-Harris算法,本文提出了匹配有效率与算法效率2种评价指标。匹配有效率由特征点(角点)的匹配对数与2幅待匹配图像(标准图像与浮动图像)中提取的特征点总数的比率表示,具体如公式(12)所示。
公式(12)
式中,PM代表2幅待匹配图像的特征点匹配对数,P1与P2分别表示从标准图像与浮动图像中提取的特征点数目,r代表比率值。
在相同实验对象的条件下,该比值越大,说明某一算法检测到的特征点越有效,反之其检测到的特征点匹配有效率较低。
算法效率指算法执行时间。算法执行时间需通过依据该算法编制的程序在计算机上运行时所消耗的时间来度量,采用特征点检测的花费时间(用t表示)评价算法的运算效率。若某一算法检测特征点耗时较少,说明该算法可以高效提取特征点,反之则说明该算法提取特征点的性能较低。
将GM-Harris算法与传统Harris算法进行比较实验,采用颅脑CT图像作为实验图像(图像来源于医学影像图库),在Windows 10系统下,以Matlab 2016A作为开发环境。具体的实验内容如下。
图1、图2分别为颅脑的标准图像与浮动图像,是从实验中选出的Harris算法与GM-Harris算法的特征点效果图。
图1 标准图像
图2 浮动图像
如表1所示,在相同实验图像的条件下,Harris算法获取的特征点(角点)数目比传统的GM-Harris算法要少。相比之下,本文提出的GM-Harris算法可以更好地解决Harris算法漏选特征点的不足,从而获得最佳的特征点提取效果。
表1 特征点数比较
图3显示了在实验图像变换不同角度时,Harris算法与GM-Harris算法获取特征点的情况,这是对两种算法特征点提取效果的进一步说明。据图3可知,在图像变换不同角度时,GM-Harris算法比Harris算法获得的特征点更充足,说明GM-Harris算法具有更好的适应性。
GM-Harris算法与Harris算法的评价指标比较见表2。通过对比表2中的r值,发现 GM-Harris算法的匹配有效率比Harris算法高6%左右。另外,通过对比t值(是进行2 000次运算后得到的均值),发现GM-Harris算法的运行时间比Harris算法节省1min左右,说明本文提出的GM-Harris算法具有较好的特征点提取性能。
图3 旋转变换与特征点数目关系
表2评价因子比较
算法r值/%t值/minHarris算法27.63.19GM-Harris算法33.42.21
综上所述,对于不同的医学图像,GM-Harris算法不但可以获得较充足的特征点(角点),还可以提高图像的配准精度(匹配有效率)以及运算效率,在一定程度上弥补了传统Harris算法的不足。
本文在实际应用中发现Harris算法检测到的特征点(角点)数量不足且图像配准精度不高。为了解决这一弊端,提出了一种优化算法(GM-Harris算法),即采用群搜索优化算法(GSO算法)与互信息相结合的方式优化传统Harris算法,并利用相关评价指标对两种算法进行定量分析。实验表明,GM-Harris算法对于不同的医学图像均具有较好的适应性,基于GM-Harris算法不仅可以获得较充足的特征点(角点),而且还能提高图像的配准精度。