杨 琴，姜自波，刘承超

(曲阜师范大学 信息科学与工程学院,山东 日照 276826)

基于频域快速解卷的血管提取算法

杨 琴，姜自波，刘承超

(曲阜师范大学 信息科学与工程学院,山东 日照 276826)

为了方便识别新生血管，从而有效地对疾病进行诊断，所以需要将血管图像从背景图像中分割出来。文中提出一种在提高血管分辨率的基础上进行血管分割的算法。首先通过特定公式在梯度图像上对血管图像进行纹理增强，这是为了在之后的图像分割中，能识别细微的血管，使其不被忽略掉；然后对模糊图像进行频域上的快速解卷去模糊，消除成像仪器在成像过程中对图像清晰度的影响；由于图像在传输过程中会产生噪声，因此为了去除噪声对血管分割的影响，接着对血管图像进行了中值滤波操作；最后使用最大类间方差法来进行血管的分割操作，因为最大类间方差法可以有效地将图像前景和背景分割开。通过实验对比，直观上证明了该算法在血管分割中的有效性。

频率域解卷；纹理增强；中值滤波；OTSU图像分割

1 概 述

虽然如今用于医学图像成像的仪器已经达到了很高的精度，例如光显微镜，它在横向空间中限制光的衍射，从而能提供微米甚至亚微米的分辨率[1-4]，但是成像仪器也具有其局限性。比如，由于受成像设备的影响、局部体效应(同一个体素中包含多种组织)、患者的体位运动和检查床的匀速直线运动，使得医学图像不可避免地盘现噪声和伪影，边缘模糊和信噪强度不均匀。此外，在图形成像和传输的过程中，图像的像质也会受到一定影响。这些都给医生下达准确的诊断造成了一定的障碍。

为了提高医学图像的可读性，使得医生可以对人体的解剖结构以及病变部位进行更有效的观察和诊断，提高诊断的准确率，医学图像的后期处理成为一个必不可少的过程。现如今许多疾病都伴随着新生血管的产生，例如视网膜表面血管的变化可能暗示着糖尿病、高血压、动脉硬化的发生。但是在疾病发生初期，血管都处于微米甚至纳米层，因此，正确提取和分析血管图像在预防和诊断病患中起到了相当大的作用。最大类间方差法[5]是较早提出的一种图像分割技术，但它对噪音十分敏感，并且仅对类间方差为单峰的图像有较好的分割效果。目前，进行血管提取的方法有很多种。例如，利用Hessian矩阵的特征值来提取血管[6-7]，该方法是一种多分辨率的提取方法。另一种进行血管提取的方法，如文献[8]中的方法，整合了多种图像增强的方法来进行血管提取，虽然得出的效果比较好，但是实验过程繁杂。

去卷积是一种计算密集型图像处理技术，它试图从模糊图像中得到对应的清晰图像。通过该技术，可以提高显微图像的对比度和清晰度[9]。目前，卷积技术可以分为两种：盲解卷积和非盲解卷积。盲解卷积需要事先知道模糊因子，而非盲解卷积则不需要知道模糊因子就可以从模糊图像中恢复出清晰图像。卷积技术的发展是值得关注的问题，目前已经得到了各研究机构的广泛关注[10-12]。

文中提出了基于频率域的非盲解卷积技术，用于提高血管的分辨率。其中的关键技术是在解卷之前使用了纹理增强技术，强化了血管边缘。该技术还能使血管和背景图像较清晰地分隔开，使得Otsu算法能够顺利进行。为了证明文中算法的可行性，分别对脑血管和视网膜血管进行了对比实验，结果显示，该算法能将血管分割出来的同时保留了部分细小血管，并且对于血管边缘的分割也比较清晰。

2 算法原理

2.1 图像解卷

对于一幅模糊图像，若直接对其进行分割操作，有些细节则不能很好地分割出来，所以需要在分割之前对其进行清晰化操作。而解卷技术可以提高显微图像的对比度和清晰度[1]，在频率域上进行解卷操作，可以大大提高操作的速度[13-14]。

假设I是原始未退化的清晰灰度图像，g是模糊或加噪声后的退化图像。文中使用卷积核k对图像I进行卷积模糊并对其添加零均值高斯噪声来获得退化图像g，用矩阵方式可表示为：g=I⊗k+n。假设g和k是已知的，目的就是构建最佳的清晰图像I。从概率方面考虑，就是为了寻找I的最大后验概率：p(I|g,k)。要得到p(I|g,k)的最大值，相当于求-logp(I|g,k)的最小值：

(1)

其中：⊗为二维卷积操作；矩阵ds,s∈{1,2,…,5}为图像I的一阶和二阶导数:{dx,dy,dxx,dxy,dyy}；δs为正权值，用来控制正规化的强度；ωs为理想图像I的过滤操作。

令τ为阈值，对于像素i，当|dsg|i<τ时，ωs=0，否则，ωs=dsg[13]。其中

(2)

(3)

(4)

其中：K=F(k),G=F(g),Ds=F(ds),Ws=F(ωs)；F(•)表示二维离散傅里叶变换；∘表示矩阵元素相乘。

令：

(5)

则式(3)可写为：

(6)

(7)

其中，./表示矩阵元素相除。

2.2 纹理增强

由于直接求图像梯度，梯度值较小的点在梯度图像上看不清，所以有必要对图像纹理做相应处理，使其在梯度图像上能较清楚地显现出来。

图像中梯度值较大的地方，光照强度值的改变也较大，需要对梯度值较大的地方进行限制，这样才能保证整幅图像具有较高的对比度。目的是减小梯度值的大小，并且还能保持它们的方向不变。首先对原始图像g(x,y)取对数，再进行相应的纹理处理。这是由于对数图像的亮度和人眼所观察到的图像亮度很接近，并且对数域里的梯度图像和人眼所观察到的局部对比度很接近[15]。定义纹理增强后的特征图为：

F(x,y)=L'(x,y)Θ(x,y)

(8)

其中：L'(x,y)为原始图像g(x,y)的对数形式；Θ(x,y)为尺度函数，用于改变梯度的大小，但不会改变梯度方向。

这里，令

(9)

然而，真实世界的图像边界是多尺度的，使用一个多分辨率的边缘检测机制可以准确、有效地检测到有意义的图像强度。但是，又不能在每次检测到图像强度的时候都能对其梯度值进行相应地增强或减弱，因为这会在边界强度较强的地方出现光环效应[15]。因此，可将一幅图像分为多层，在每一层中使用一个梯度减弱系数。这样，所有的梯度增强或减弱操作都在单个层中进行，而每一层是针对整幅图像进行的，这样就不会出现光环效应。

首先，构建一个k层的高斯金字塔，在每一层中，计算其梯度值：

在第k层中，尺度函数定义为：

(10)

正如文献[15]中所说的，参数α决定哪个梯度不被改变，参数β决定哪个梯度被减小(β<1)，即梯度大的被减小，比α小的梯度被稍微放大。尺度函数Θ(x,y)使用递归的方法，将每层得到的尺度函数进行逐个相乘并相加：

Θ1(x,y)=ϑ1(x,y)

Θk+1(x,y)=Θk(x,y)+Θk(x,y)ϑk+1(x,y)

(11)

2.3 去 噪

在进行血管分割之前，首先需要进行一定程度的降噪，因为噪声会对血管分割产生一定程度的影响。目前用于去噪的方法有很多种。例如，均值滤波器，它在滤波过程中增加了图像的平滑度，但也丢失了一部分细节；自适应维纳滤波，它对保留图像边缘和高频部分有较大的作用，但是计算量较大；形态学滤波，它适用于图像中的对象尺寸都比较大，且没有细小的细节的图像。

文中采用中值滤波器。中值滤波是基于排序统计理论的一种能有效抑制噪声的非线性数字滤波技术，经常用于去除图像或者其他信号中的噪声。它可以在滤除噪声的同时，保护信号的边缘，使之不被模糊。中值滤波的基本原理是把数字图像或数字序列中一点的值用该点的一个邻域中各点值的中值代替，让周围的像素值接近真实值，从而消除孤立的噪声点。其步骤为：

(1)从图像中的某个采样窗口取出奇数个数据并从小到大排序；

(2)用排序后的中值取代要处理的数据。

2.4 血管分割

基于去模糊和去噪的血管图像，将血管和背景图像分割开，需要选择合适的分割算法。现如今，用于图像分割的方法有很多种。如阈值分割，这种方法计算简单、运算效率高、速度快，它的关键技术在于阈值的选取；基于区域的分割，区域生长技术需要选择一组能正确代表同一区域的种子像素，确定在生长过程中的相似性准则，制定让生长停止的条件或准则。这种技术对噪声比较敏感，可能导致区域内有空洞；区域分裂合并技术，它是区域生长的逆过程。这种方法对复杂图像的分割效果较好，但算法较复杂，计算量大，分裂还可能破坏区域的边界；边缘检测技术，它常用微分算子进行边缘检测，但是这些算子对噪声敏感，因此这种技术只适合于噪声较小、不太复杂的图像。

文中选用Otsu法(最大类间方差法)。它使用的是聚类的思想，把图像的灰度数按灰度级分成两个部分，使得两个部分之间的灰度值差异最大，每个部分内部的灰度差异最小。Otsu法通过方差的计算来寻找一个合适的灰度级别进行分类，并且此算法不受图像亮度和对比度的影响。前景和背景图像的方差f由下面的过程得到：

记t为前景与背景的分割阈值，前景点数占图像总点数的比例为w0，平均灰度为u0；背景点数占图像总点数的比例为w1，平均灰度为u1。

图像的总平均灰度为：

u=w0u0+w1u1

前景和背景图像的方差为：

f=w0*(u0-u)*(u0-u)+w1*(u1-u)*(u1-u)

从N个灰度级遍历t，使得t为某个值的时候，前景和背景的方差f(t)最大，则这个t值便是要求得的阈值。

基于频率域解卷的血管提取方法描述如下：

数据：g为所获得的退化图像；k为卷积核。

步骤1：对g求对数L=log(g+1);

步骤2：求对数图像L的梯度，利用式(9)得到梯度图像;

步骤3：构建高斯金字塔，利用式(10)计算每层的尺度函数ϑ，利用式(11)计算最终的尺度函数Θ；

步骤4：利用式(8)得到最终的纹理增强图像F;

步骤5：对图像F求一阶和二阶导数；

步骤6：利用式(2)和计算得到的F的各阶导数计算w；

步骤7：利用式(5)、(7)得到最终的解卷图像I；

步骤8：对I进行平滑滤波，得到去噪图像I';

步骤9：对去噪后的图像进行血管分割。

3 实 验

文中算法流程如图1所示。

图1 文中算法流程

图中，(a)为模糊的原始图像;(b)为纹理增强的梯度图像;(c)为基于频域解卷的图像;(d)为滤波后的图像;(e)为血管分割图像。

首先，对退化图像进行纹理增强，这有助于细微血管的提取；然后，对模糊图像进行清晰化处理，由于噪音会对血管分割产生一定的影响，所以有必要在提取血管前进行滤波操作。

3.1 实验结果

将Horacio等提出的图像恢复算法直接应用到血管提取中，会发现有些血管被忽略了，如图2所示。

图2 脑血管分割

可以发现，在图(c)中，只分割出血管主干，某些细微血管则被认为是背景忽略掉。在图(d)中，某些灰度值较小的血管没有能分割出来，出现了血管的不连续性。对于文中算法，它在将血管主干提取出来的同时，仍然保留了某些细微血管，这对血管疾病的诊断是非常有必要的。

为了进一步说明文中算法的优越性，使用了另一幅血管图像进行了对比实验，如图3所示。

可以发现，使用Otsu算法可以将主干血管分割出来，但是细微血管却被忽略甚至多根血管被重叠为一根血管。而图(d)的血管边缘模糊，这对判断血管粗细程度造成了影响。通过对比发现，文中算法分割出来的血管纹理清晰。

3.2 参数选择

利用文献[8]中所使用的19×19卷积核进行实验。针对脑血管图像，设置α等于0.000 05倍的平均梯度，β设为0.9时纹理增强效果较好。解卷过程中，对于所有的s，将λ均设为0.000 05，对于一阶导数{dx,dy}，设置τ(1)=0.065，对于二阶导数{dxx,dxy,dyy}，设置τ(2)=0.0325时，实验结果较好。针对视网膜毛细血管，将λ设为0.000 01，其他参数不变。

图3 视网膜毛细血管分割

4 结束语

文中从造成图像不清晰的因素考虑，提出了基于频率域解卷的血管提取方法。同时考虑到Otsu算法进行图像分割的局限性，所以在使用Otsu算法进行血管分割之前，先将血管纹理进行了加强，将背景图像进行了一定程度的灰度值变化，便于后面运用Otsu算法对血管进行分割。纹理增强的关键在于α和β的选取，它们控制背景和前景血管的亮度对比。实验结果表明，在进行血管分割之前，对图像进行解卷以提高图像的分辨率是有必要的。

该算法的局限性在于：只能对背景色相对血管较暗的血管图片进行血管提取。因此，适用性更广、精确度更高的血管提取算法是进一步研究的目标。

[1]MaR,SöntgesS,ShohamS,etal.Fastscanningcoaxialoptoacousticmicroscopy[J].BiomedicalOpticsExpress,2012,3(7):1724-1731.

[2]MaslovK,ZhangHF,HuS,etal.Optical-resolutionphotoacousticmicroscopyforinvivoimagingofsinglecapillaries[J].OpticsLetters,2008,33(9):929-931.

[3]SongL,MaslovK,WangLV.Multifocaloptical-resolutionphotoacousticmicroscopyinvivo[J].OpticsLetters,2011,36(7):1236-1238.

[4]HajirezaP,ShiW,ZempR.Label-freeinvivoGRIN-lensopticalresolutionphotoacousticmicro-endoscopy[J].LaserPhysicsLetters,2013,10(5):055603.

[5]OtsuN.Athresholdselectionmethodfromgray-levelhistograms[J].IEEETransonSystems,ManandCybernetics,1979,9(1):62-66.

[6]SatoY,ShiragaN,ShiragaN,etal.Three-dimensionalmulti-scalelinefilterforsegmentationandvisualizationofcurvilinearstructuresinmedicalimages[J].MedicalImageAnalysis,1998,2(2):143-168.

[7]FrangiAF,NiessenWJ,VinckenKL,etal.Multiscalevesselenhancementfiltering[C]//ProcofMICCAI.[s.l.]:[s.n.],1998:130-137.

[8]YangZY.Multi-parametricquantitativemicrovascularimagingwithoptical-resolutionphotoacousticmicroscopyinvivo[J].OpticsExpress,2014,22(2):1500-1511.

[9]ChenJH,LinRQ,WangHN.Blind-deconvolutionoptical-resolutionphotoacousticmicroscopyinvivo[J].OpticsExpress,2013,21(6):7316-7327.

[10]BanhamM,KatsaggelosA.Digitalimagerestoration[J].IEEESignalProcessingMagzine,1997,14(2):24-41.

[11]CampisiP,EgiazarianK.Blindimagedeconvolution:theoryandapplications[M].BocaRaton:CRCPress,2007.

[12]LevinA,WeissY,DurandFW.Efficientmarginallikelihoodoptimizationinblinddeconvolution[C]//ProcofCVPR.[s.l.]:[s.n.],2011:2657-2664.

[13]FortunatoHE,OliveiraMM.Fasthigh-qualitynon-blinddeconvolutionusingsparseadaptivepriors[J].VisualComputer,2014,30(6):661-671.

[14]DilipK,RobF.Fastimagedeconvolutionusinghyper-laplacianpriors[C]//ProcofNIPS.[s.l.]:[s.n.],2009:1033-1041.

[15]FattalR,LischinskiD,WermanM.Gradientdomainhighdynamicrangecompression[J].ACMTransonGraphics,2002,21(3):249-256.

Vessel Extraction Algorithm Based on Fast Frequency-domain Deconvolution

YANG Qin,JIANG Zi-bo,LIU Cheng-chao

(College of Information Science and Engineering,Qufu Normal University,Rizhao 276826,China)

In order to identify new blood vessels,thus effectively diagnosing disease,it is necessary to segment the vascular image from the background image.A segmentation algorithm based on improved resolution of the vessel is proposed.First the texture is enhanced for blood vessel image through specific formula in the gradient image,in order to identify small blood vessels and not ignore it in the image segmentation later.Secondly,the image deconvolution in frequency domain is made to obtain a sharp image,eliminating the influence on image clarity of imaging instrument in the imaging process.Then so as to remove the influence of noise on the vessel segmentation,a filtering operation is carried on vascular image.Finally,the maximum class variance method is used to segment vessels,because it can effectively distinguish the image foreground and background.The effectiveness of the algorithm in vessel segmentation is verified intuitively by experimental comparisons.

frequency domain deconvolution;texture enhancement;median filtering;OTSU image segmentation

2015-06-16

2015-09-18

时间：2016-02-18

国家自然科学基金资助项目(61572283)；山东省自然科学基金(ZR2011FQ026)

杨 琴(1990-)，女，硕士研究生，CCF会员，研究方向为图像处理与模式识别。

http://www.cnki.net/kcms/detail/61.1450.TP.20160218.1634.050.html

TP391

A

1673-629X(2016)03-0113-04

10.3969/j.issn.1673-629X.2016.03.027