1 引言

医学图像对医生进行病症诊断和治疗有至关重要的作用,然而医学图像的成像质量与选择成像的方式、设备以及在成像的过程中的可能遇到的误差密不可分,这关系到医生对病症诊断的准确性^[1]。医学图像在采集、传输和存储的过程中不可避免地会产生大量噪声,降低了图像的质量并使视觉效果变差,不利于后期图像病症的识别。为了便于医生判断和处理,提高医生诊断的正确性,需要对图像的纹理和细节信息进行增强处理,以降低噪声的干扰。

偏微分方程在图像处理和模式识别领域是一种新的研究方向^[2],由于偏微分方程具有对图像的细节进行提升的特点,因此已广泛应用于图像处理中。目前,国内外的研究学者已经利用偏微分方程在图像去噪^[3-6]、图像增强^[7-9]和图像重建^[10-12]等领域取得了一些研究成果。偏微分的图像增强方法,包括整数阶微分、分数阶微分以及有理数阶微分。整数阶的偏微分方程方法已经在图像重建领域中取得了较好的研究成果^[13],但是经过此方法处理后的图像容易丢失部分纹理信息,在保留图像的纹理和细节方面的处理效果不好,于是分数阶微分作为整数阶微分的扩展应运而生。许多学者对分数阶微分进行了大量研究,尤其是在对图像进行处理的方面。Pu等^[14-16]对分数阶微积分进行了分析,得出了分数阶微分对增强图像纹理和细节信息具有明显效果的结论,首次将分数阶微分引入到了数字图像增强处理当中。黄果等^[17-18]提出了一种可以根据图像特征和信息动态地调整分数阶阶次的图像增强算法。汪成亮等^[19-20]提出了一种可以根据图像的信息自动调整分数阶阶次的自适应图像增强算法。李军成^[21]首次将0~1 阶范围的分数阶微分阶次扩展到了1~2阶。蒋伟等^[22-23]结合了分数阶微分和整数阶微分各自的优点,将分数阶微分扩展到了有理数阶微分。

图像增强是一个需要持续研究的问题,上述方法在图像增强方面有一定的效果,但是在细节和纹理的增强效果方面还有待改进。而局部二值模式(LBP)^[24-25]可以很好地描述图像的纹理特性,也已广泛应用到图像识别^[26]和图像匹配^[27-28]中,因此,本文提出了一种新的基于分数阶微分理论和局部二进制模式方差(LBPV)相结合的医学图像增强算法。运用LBPV理论对图像进行特征提取,找出待处理图像的主方向,对其分数阶微分算子模板系数进行加权增强,再将该方法应用到医学图像中,并与现有的分数阶微分方法进行实验对比。实验结果表明,所提方法在图像纹理和细节处理的效果比现有的分数阶算法效果更好,具有良好的视觉效果。

2 理论知识

2.1 分数阶微分的基本定义

分数阶微积分的Grünwald-Letnikov的定义^[14-16]是从研究连续函数整数阶导数的经典定义出发,将微积分由原来的整数阶n∈Z⁺扩展到分数阶v而得到的,∀v∈R,令其整数部分为[v]。若信号f(t)∈[a,t](t,a分别为微分的上下限,a<t,a∈R,t∈R),则存在m+1(m∈Z)阶连续导数;当v>0时,m至少取v的整数部分[v],则分数阶(v阶)导数的定义为

aGDtv=limh→01hv∑m=0t-ah(-1)mΓ(v+1)m!Γ(v-m+1)f(t-mh),(1)

式中:G为G-L的定义;v为分数阶微分阶次;h为步长;Γ(v)为Gamma函数。为使f(t)为非零极限,需当h→0时,n→∞,故令h=(t-a)/n,于是n=[(t-a)/h]。将信号的持续期[a,t]按单位h=1进行划分,可以推导出一元信号f(t)的分数阶微分表达式为

dvf(t)dtv≈f(t)+(-v)f(t-1)+(-v)(-v+1)2!f(t-2)+…+Γ(-v+1)n!Γ(-v+n+1)f(t-n)。(2)

2.2 LBP概述

LBP是一种用来度量和提取局部图像的纹理信息的有效纹理描述算子。初始定义是在某中心像素及其周围大小为3×3的矩形邻域系统中,把中心点像素的灰度值作为阈值,将周围八个邻域的像素点和中心像素点进行比较,如果邻域点的灰度值小于阈值,则该点用0表示,大于或等于阈值则用1表示,这样得到一个二值图,然后从x正轴方向开始,按逆时针方向进行串联得到一个二进制串,即为LBP模式,将该二进制字符串的数字转换为相应的十进制数,并用十进制数来标识像素点,最终得到的值即为LBP算子的值,如图1所示。

图 1. LBP模式

Fig. 1. LBP mode

下载图片 查看所有图片

由于原始的LBP算子存在对尺度较大的纹理结构无法捕获的问题,并且无法实现旋转不变的设计。Ojala等^[24]又提出了一种采用圆形邻域结构的方式计算LBP模式,如图2所示。定义在该圆形邻域系统上的LBP模式计算式为

LpR=∑p=0P-1s(gp-gc)2p,(3)s(x)=1,x≥00,x<0,(4)

图 2. 圆形邻域系统

Fig. 2. Circular neighborhood system

下载图片 查看所有图片

式中:s(x)是一个符号函数;g_c是中心像素点(i,j)的灰度值;g_p是邻域像素点的值;P是邻域像素的个数;R是邻域的半径。

2.3 LBPV算子

在识别图像的纹理中,通常采用空间结构(例如,LBP模式)和对比度[例如,局部图像纹理的方差(VAR)]描述该区域的纹理特征。然而,空间结构随着旋转的变化而变化,对比度则不然。所以Ojala等^[24]提出了联合直方图,使用两个互补特征LBP/VAR用于旋转不变纹理分类。然而在这种方法中,VAR值是连续的,需要通过量化步骤计算直方图,但是在量化过程中不可避免地引起一些误差,这会影响最终的分类效率,所以引入了LBPV^[25],其不是全局计算LBP和VAR的联合直方图,而是使用局部区域方差的权值作为LBP子模式的权重,有效地将图像的空间模式与局部差异性融为一体,无需量化过程,避免了量化过程中可能出现的问题。另外,LBPV提取的特征向量维数相对较小,计算复杂度低,可以有效节省计算时间。

3 所提方法

3.1 特征提取

假设有一幅大小为M×N的纹理图像,在识别每个像素(i,j)的LBPV模式之后,通过构建直方图表示整个纹理图像为

WPR(k)=∑i=1M∑j=1Nf[LPR(i,j),k],k∈[0,K],(5)fLPR(i,j),k=VPR(i,j),LPR(i,j)=k0,others,(6)VPR=1P∑p=0P-1(gp-u)2,(7)u=1P∑p=0P-1gp,(8)

式中:u为P邻域内的平均灰度;k为LBPV模式;K为LBPV模式的最大值。其中当P=8时的统一LBP模式中的分类如图3所示,共有58种,即K=58,其中黑点和白点分别代表1和0。例如,当全为白点时,LBPV模式为00000000,此时k=0。其他模式(跳变值U>2)归为混合模式,共有一种。图3的前七行中每一列代表同一方向包括的模式,每一行包括8个方向,最后一行代表全0和全1模式。

图 3. 当P=8时的统一LBP模式图

Fig. 3. Diagram of unified LBP mode when P=8

下载图片 查看所有图片

3.2 主方向估测

将主方向定义为模式频率最高的方向。图4为根据上述特征提取的公式计算得到的纹理图像的模式与模式频率之间关系图。而由图3、4的对应关系可得出纹理图像的角度与模式频率的关系图,如图5所示。从图中可以清晰看出角度与模式频率的关系,进而找出频率最大的方向。经过实验发现选取一个主方向时的增强效果没有选取两个主方向时好,

图 4. 纹理图像的模式与模式频率图

Fig. 4. Mode and mode frequency map of textured image

下载图片 查看所有图片

图 5. 纹理图像的角度与模式频率图

Fig. 5. Angle and mode frequency map of textured image

下载图片 查看所有图片

而选取3个主方向时,图像会变得立体,不适合所选图像,因此选择两个主方向。主方向即为出现频率最高的模式方向,分别选择270°和90°作为主方向。不同图像增强时,可能选取的主方向角度以及个数都不同,因此,需要先找到图像的所有方向,按图像的具体要求选取主方向的角度与个数。

3.3 改进后的分数阶掩模

根据(2)式可以推导出f(x)在x和y坐标轴上分数阶偏微分的差分近似表达式分别为

dvf(x,y)dxv≈f(x,y)+(-v)f(x+1,y)+(-v)(-v+1)2!f(x+2,y)+(-v)(-v+1)(-v+2)3!f(x+3,y)+…+Γ(-v+1)n!Γ(-v+n+1)f(x+n,y),(9)dvf(x,y)dyv≈f(x,y)+(-v)f(x,y+1)+(-v)(-v+1)2!f(x,y+2)+(-v)(-v+1)(-v+2)3!f(x,y+3)+…Γ(-v+1)n!Γ(-v+n+1)f(x,y+n),(10)

式中:n为趋于无穷的正整数。

由差分表达式可以得到分数阶微分的差分系数为

a0=1a1=-va2=(-v)(-v+1)2!a3=(-v)(-v+1)(-v+2)3! ︙an=Γ(-v+1)n!Γ(-v+n+1)。(11)

因为图像相邻的像素点之间相关性很强,为了使增强后的算法具有旋转不变的性质,分别构造出8个方向(0°、45°、90°、135°、180°、225°、270°、315°)的掩模模板,为了保持方向一致性,掩模算子方向由LBPV的方向决定,如图6所示。同时,构造分数阶模板,如图7所示。

图 6. LBPV的8个方向

Fig. 6. Eight directions of LBPV

下载图片 查看所有图片

图 7. 分数阶的8方向掩模。(a) 0°;(b) 45°;(c) 90°;(d) 135°;(e) 180°;(f) 225°;(g) 270°;(h) 315°

Fig. 7. Fractional masks for eight directions. (a) 0°; (b) 45°; (c) 90°; (d) 135°; (e) 180°; (f) 225°; (g) 270°; (h) 315°

下载图片 查看所有图片

一般来说,对M×N的图像f(x,y)增强,可用m×n大小的模板对其进行线性滤波,令模板中心坐标为w(0,0),则原有的图像s(x,y)表示为

s(x,y)=∑o=-aa∑r=-bbw(o,r)f(x+o,y+r),(12)

式中:a、b分别为分数阶微分算子模板大小在 x、y坐标轴方向上的半径;(o,r)为模板坐标;w(o,r)为滤波系数。

所采用的掩模为利用LBPV通过特征提取得到的图像主方向,将掩模算子对应图像主方向的方向系数加权处理,其他方向的系数不变,最终得到掩模算子。由于图像不同,每幅图像的主方向也不相同,所以采用的掩模也不相同。改进后的增强图像g(x,y)表示为

g(x,y)=∑o=-aa∑r=-bbw'(o,r)f(x+o,y+r),(13)w'(m,r)=w(o,r),maindirectionisinconsistent withthedirectionofthefractionalmasklw(o,r),maindirectionisconsistent withthedirectionofthefractionalmask,(14)

式中:w'为新的滤波系数。选取的加权倍数l=2。即当分数阶微分算子掩模方向与图像主方向相同时,l=2;当算子掩模方向与主方向不同时,则l=1。

4 实验结果

为了验证算法在图像增强方面的有效性,采用所提算法对outex数据集中的图片进行增强处理,并采取主观评价和指标评价两种方法来评价图像增强的效果。又将该算法应用到医学图像中,并与现有的分数阶算法进行了对比分析。从结果可以看出,与其他方法相比,所提算法在图像增强上有良好的视觉效果,图像纹理和细节信息更加突出。选取平均梯度(AG)和信息熵(E)两个参量作为评价指标。平均梯度反映图像的空间变化率,其值越大,表明图像变化越明显。设一幅大小为M×N的图像f(x,y),其计算公式为

g=1(M-1)(N-1)∑x=1M∑y=1N12∂f(x,y)∂x2+∂f(x,y)∂y2。(15)

信息熵反映图像所含信息的丰富程度,其值越大,表明图像所含信息量越大,图像的细节信息越丰富等。信息熵的计算公式为

E=-∑k=0255p(q)lb[p(q)],(16)

式中:p(q)为图像中像素q出现的概率。

4.1 在Outex数据集中的实验结果

为了证明算法的实用性,采用所提算法对Outex数据集中的图片进行处理,并从数据集中选取了部分图像来检验该算法的有效性。

将图8中原图与增强后的图像对比可以看出,增强后的图像在细节和纹理等方面具有良好的增强效果,并包含更多的纹理和细节信息,图像的对比度也有了明显提高,图像的边界更加清晰并具有良好的视觉效果。此外,表1中的数据与原始图像相比,平均梯度和信息熵均得到了较大提升,从而可以证明该算法在图像增强中的有效性。

4.2 在医学图像中的实验结果

为检验所提算法在医学图像上的效果,将其与文献[ 5]、 [6]和[9]算法进行比较,选取大量术后的肱骨头粉碎性骨折图像和腓骨图像进行实验,并从图像库中选取一幅图像进行增强,仍采用平均梯度和信息熵两个参数来评价图像增强的效果。为了更细致地观察骨折区域的细节信息,还对骨折位置的细节信息进行了放大,如图9、10所示,表2、3分别给出了图9、10的评价参数。

由图9、10可以看出,所提算法的视觉效果优于其他方法,在图像纹理和细节方面的增强效果较好,提高了图像边缘的清晰度。从图9中增强效果的局部图像可以看出,所提算法处理后的骨折区域更为明显(椭圆部分为钉子,两侧钉子中心即为骨折区域),并且含有更丰富的纹理和细节。从图10中增强效果的局部图像也可以看出,所提算法处理后的骨折区域边缘更加明显。由表2、3可以看出,所提算法的信息熵和平均梯度高于文献[ 5]和[12]算法,虽然在评价参数中文献[ 9]算法的平均梯度值更高,但是结合图9、10可以看出,经过文献[ 9]算法处理后的图像含有大量噪声点,存在过度增强的现象,对图像处理的效果不好,因此,所提算法更适合用于医学图像的增强处理。

图 8. 原图与增强后的图像。(a)原图;(b)图8(a)增强后的图像;(c)原图;(d)图8(c)增强后的图像

Fig. 8. Original images and enhanced images. (a) Original images; (b) enhanced images of Fig. 8(a); (c) original images; (d) enhanced images of Fig. 8(c)

下载图片 查看所有图片

图 9. 不同方法下肱骨头的增强图像及其局部图像。(a)原图;(b)文献[ 5]方法;(c)文献[ 6]方法;(d)文献[ 9]方法;(e)所提方法

Fig. 9. Enhanced images of humeral head and its local images with different methods. (a) Original images; (b) method in Ref. [5]; (c) method in Ref. [6]; (d) method in Ref. [9]; (e) proposed method

下载图片 查看所有图片

图 10. 不同方法下腓骨的增强图像及其局部图像。(a)原图;(b)文献[ 5]方法;(c)文献[ 6]方法;(d)文献[ 9]方法;(e)所提方法

Fig. 10. Enhanced images of tibia and its local images with different methods. (a) Original images; (b) method in Ref. [5]; (c) method in Ref. [6]; (d) method in Ref. [9]; (e) proposed method

下载图片 查看所有图片

5 结论

为了有效地对医学图像进行增强处理,提出了一种新的基于分数阶微分理论和LBPV相结合的图像增强算法,并将其应用于医学图像增强中,将所提算法与现有的分数阶微分图像增强算法进行了比较分析,实验结果表明,所提算法在增强图像的细节和纹理信息方面具有良好的效果,提高了图像的对比度,增强了图像的清晰度,对图像有较好的增强效果。与现有的分数阶图像增强算法相比,所提算法增强后的图像在质量和视觉效果等方面均有明显提高,较好地达到了增强图像细节处理的目的,是一种有效的医学图像增强方法。