基于图像分割的视网膜血管图像配准研究

刘志东, 林江莉, 陈 科*

(1. 四川大学 材料科学与工程学院, 四川 成都 610065; 2. 四川城市职业学院 汽车与信息工程学院, 四川 成都 610010)

基于图像分割的视网膜血管图像配准研究

刘志东1,2, 林江莉1, 陈 科1*

(1. 四川大学 材料科学与工程学院, 四川 成都 610065; 2. 四川城市职业学院 汽车与信息工程学院, 四川 成都 610010)

利用不同波长的视网膜图像可以对视网膜血管血氧饱和度进行计算,但需进行配准处理.提出一种基于视网膜图像血管分割的互信息图像配准方法.为了充分利用血管的轮廓信息和灰度信息,提高配准精度,首先对配准图像进行图像分割,提取视网膜图像中的血管轮廓信息；然后对分割后图像中的血管进行相似度计算,并采用Powell优化算法中的黄金分割法一维搜索算法来提升运算速度；最后根据计算的相似度值来完成不同波长图像的配准.研究中算法配准获得变换参数(角度、X方向、Y方向)的误差的均值分别为2.00%、2.53%和2.52%,误差的方差分别为0.57、2.09和0.34,均优于直接互信息配准方法.实验表明：算法可以自动、有效地对不同波长的视网膜血管图像进行配准,并具有良好的可重复性和稳定性.

视网膜图像; 图像分割; 配准; 互信息

人眼的很多疾病,如青光眼、视网膜病变、血管阻塞等,都会引起视网膜血管血氧饱和度值的变化,因此对视网膜血管中血氧饱和度值的无损测量是目前研究的热点问题.A. Harris等[1]根据组织对光谱的吸收特性的研究,提出了基于光谱法的无损视网膜血氧测量方法,该方法根据含氧血红蛋白和脱氧血红蛋白对2个波长的光吸收差异,同时获取2个不同波长的视网膜图像,但由于2幅图像是非共光路引入,因此图像位置会出现偏差,在计算血氧值之前,必须对2幅图像进行配准.

到目前为止,图像配准的方法很多,按照配准过程基本上可以分为两类：基于图像灰度信息和基于图像特征[2-3].基于图像特征的配准方法有一定的不足之处,在特征提取过程中很多情况需要人工干预,同时图像特征的提取本身就是一个比较困难的问题,而且特征提取的精确度也无法保证,会直接影响配准的效果[4].基于图像灰度信息的方法主要利用的是图像中像素的灰度信息,根据概率密度计算公式对2幅图像之间的相似度进行计算,通过优化搜索算法获得配准图像间的变换参数.

基于灰度信息配准的方法,目前主要有联合熵法、相关熵法、最大互信息法、条件熵法等.与基于图像特征配准方法相比,基于灰度信息的配准方法是一种自动的方法,无需人工干预,依赖于图像本身信息进行配准,同时可以适用于不同模态图像的配准[5-7],避免了主观因素的影响.如洪明坚等[8]结合图像二维信息与互信息对图像进行配准；翟利志等[9]利用邻域信息与互信息相结合处理红外与可见光图像配准的问题.越来越多的研究将灰度信息与图像目标结构信息相结合,充分利用图像中存在的纹理、形状、结构等信息,提高配准的准确性是目前研究的热点.

本文针对2幅不同波长视网膜图像进行配准,充分利用图像中血管轮廓信息,采用基于图像分割的互信息配准方法进行配准研究.为了充分利用血管的轮廓信息和减少图像配准中无关信息对配准结果的影响,首先对配准图像进行图像分割,提取视网膜图像中的血管信息；对分割后图像中的血管进行相似度计算,并采用Powell优化算法中的黄金分割法一维搜索算法来提升运算速度；根据计算的相似度值来完成不同波长图像的配准.本文的算法流程如图1所示.

1 图像分割

1.1 图像预处理 根据视网膜血管不同组织结构对光谱的吸收特性不同,血管内部和血管壁在图像中呈现不同的灰度级别,尤其在一些大尺寸的血管中,存在中央反射的问题,为血管分割带来了一定的困难.本文采用数学形态学中的开运算来对血管进行填充,将血管内部和血管壁的灰度值调整到一个级别上去,可以有效去除中央反射问题.同时,开运算还会对图像起到滤波作用,将图像中存在的一些椒盐噪声滤除.视网膜血管在图像中的像素值低于周围眼部其他组织,为了在血管分割过程中提高视网膜血管的可见性,将待处理图像进行图像反转[10].在灰度图像中像素灰度值为0～255,而所谓反转就是指将原灰度值为g的像素值转变为255-g,即原来像素值为0,反转后就变为255,原来为255,反转后就变为0.通过反转,可以将视网膜血管的亮度进行很大程度的提升,为后续处理奠定基础.

1.2 血管增强 视网膜图像中的视盘和其他相关组织,在血管分割过程中都会对分割结果产生一定的影响.,为了避免这些影响,本文采用top-hat变换对血管对比度进行增强[11].

top-hat变换定义为

(1)

其中,IA为反转后图像,O(IA,B)为图像IA数学形态学开运算后结果,运算的结构元素为B,Itop-hat为top-hat变换后图像.top-hat变换的实现原理为反转后图像的血管的灰度值要大于视盘及周边组织的灰度值.通过选择一定尺度的结构元素B对图像进行开运算,图像中的血管信息将被周边组织的信息所覆盖,得到血管被模糊掉的图像,保留了非血管图像信息,再在原图像的基础上减去非血管图像信息,使得图像血管信息与周边其组织信息的对比度进一步提高,从而增强了血管信息.在变换中,结构元素B的尺寸选择就很重要,要保证结构元素B在图像中做开运算时,能够覆盖直径最大的血管,这样才能将图像中的血管信息全部去除,提高变换效果.通过测量本文研究中采用的图像血管最大直径一般在14像素左右,根据结构元素直径要大于血管直径的原则,本文选择的结构元素直径为16像素.血管增强效果如图2(b)所示.

1.3 多窗口检测 单窗口检测首先需要定义一个宽和高都为H像素的检测窗口,检测窗口在被检测图像上滑动并遍历整幅图像,在遍历过程中通过获取最有意义检测线来获取图像血管的特征.单窗口检测的过程如下：

1) 定义检测窗口,宽和高同为H像素;

(2)

其中,P(i,j)为检测窗口覆盖下图像像素灰度值,NH为检测窗口覆盖的像素个数总数;

(3)

其中,P(i,j)为检测窗口检测线上的像素灰度值;

(4)

为了使血管特征函数更好地反映出血管特征,检测窗口的大小选择就比较重要,检测窗口的中心处于血管中心线上时,窗口内的血管像素和背景像素数据近似,因此窗口H一般为血管宽度的2倍,而检测线长度一般与血管宽度相当.但是该方法某些情况下效果不好,如血管紧邻时,容易将二者合并为一条血管.

针对单窗口检测对相邻血管处理效果不好等存在的问题,本文采用多窗口检测的方法,通过改变窗口的宽度,控制检测线的长度,构建多窗口检测器.特定检测窗口长度h下的血管特征响应函数[12]为

(5)

其中,H为最大检测窗口的宽度,h的范围为1≤h≤H,表示检测窗口的宽度.不同宽度的检测窗口在血管特征提取过程中具有不同的作用,长宽度的检测窗口可以有效的对大尺寸的血管进行特征提取,而短宽度的检测窗口则可以有效的处理紧邻的血管,避免了单窗口检测的相邻血管合并的问题.采用线性组合的方式,对不同宽度窗口检测的优点进行合并,获得多窗口检测的最终响应函数

(6)

其中,Nh为不同宽度的检测窗口的数量.

1.4 后处理 图2展示了图像血管分割过程.

通过多窗口检测操作,将图像血管特征进行提取,这里所谓提取并不是已经对图像中的血管进行了分割,而是通过检测操作,获取了图像每一个像素点属于血管或背景的可能性,即检测后的图像像素值集中在2个区域,即图像的直方图集中呈现为2个分明的峰值,有了这个基础,就可以很容易利用阈值分割方法对图像进行分割处理.但是分割结果中存在很多的小的分割区域,这些小的区域在后续配准中存在一定的影响,必须将其去除,如图2(c)所示.本文根据血管的连通性,将连通区域的面积小于100的连通区域去除,如图2(d)所示.

2 图像配准

在血管分割的基础上,获取分割后的血管图像,利用血管轮廓信息与灰度信息,采用互信息配准方法对不同波长下的视网膜图像进行配准[13].互信息用来描述2个系统间的信息相关性,用熵进行定量的数学描述.

图像A的熵定义为

(7)

图像B的熵定义为

(8)

图像A、B的联合熵定义为

3）可以通过搜索引擎，在百度、搜狐上搜索一些熟知的英文新闻网站和英文学习网站，了解最新国内外大事和与四、六级考试相关的资讯。

(9)

其中,PA(i)、PB(j)分别表示图像A、B的灰度概率密度分布,PAB(i,j)表示图像A、B的联合灰度概率分布.

图像A、B的最大互信息为

(10)

根据统计学分析,如果图像A和图像B是完全独立的,则

(11)

同时I(A,B)=0；如果图像A和图像B完全依赖,则PAB(i,j)=PA(i)=PB(j),即H(A)=H(B)=H(A,B),则I(A,B)最大.对于图像A和图像B,当二者的空间位置一样时,2幅图像间的联合概率密度分布最集中.从2幅图像之间关联性角度分析,则指2幅图像间的相关性最大.而从熵的角度来看,即最终计算出来的2幅图像之间的联合熵的值最小,表明2幅图像之间的互信息最大,这个时候表明2幅图像已经配准了.

本研究中互信息配准流程如图3所示.

图像配准目的就是寻找合适的变换参数,所以其过程就是寻找最优变换参数的过程.所谓最优变换参数就是使得2幅图像之间的互信息最大.但是,互信息最大值的计算和搜索是比较慢的,目前应用于互信息计算的优化算法很多,Powell优化算法就是其中一种[14].本文采用Powell优化算法中的一维黄金分割搜索算法.一维黄金分割搜索算法的思路就是,先利用黄金分割法确定一个较小的包含极小点的不确定区间,然后利用抛物线法获得一个极小点,若此极小点落在此不确定区间,则利用该极小点继续进行二次插值；否则放弃该点,改用黄金分割法继续搜索,经过多次迭代后便可求得互信息的最大值[15-16].

3 实验结果

表1中共列举了5组实验数据.将原图与旋转后的图像直接采用直接互信息配准方法和本文中的配准方法进行配准,并计算每一组的配准误差.将配准误差进行均值和方差计算,其中采用直接互信息配准方法的角度、X方向和Y方向参数误差的均值分别为9.68%、9.14%和10.40%,误差的方差分别为2.56、4.30和15.03.采用本文算法的误差均值分别为2.00%、2.53%和2.52%,误差的方差分别为0.57、2.09和0.34.通过对比可以看出,本文算法无论准确性和算法的鲁棒性都优于直接互信息配准方法.

表 1 方式-图像配准参数表

数据实际变换参数角度/(°)XY直接互信息参数角度/(°)XY本文算法参数角度/(°)XY直接方法误差/%角度/(°)XY本文方法误差/%角度/(°)XY15．0010．0013．004．5811．0314．235．079．8513．288．4010．309．461．401．502．1029．00-12．00-8．009．72-13．43-9．899．11-11．69-8．158．0011．9111．131．202．581．8836．00-14．007．005．38-13．018．165．83-13．786．8910．307．0716．572．831．571．574-7．0011．0013．00-6．329．9811．87-7．1311．5513．369．719．278．691．865．002．775-8．00-12．00-16．00-8．96-11．14-16．99-8．22-11．76-16．4712．007．176．192．752．002．94

第二种方式采用设备采集的570nm图像和600nm图像作为配准数据,但是该组数据的旋转和平移的变换参数无法未知,无法采用方式一中的方法进行配准结果的定量判断.基于此,本文采用直观观察的方式进行效果验证,即通过将570nm图像分割得到的血管轮廓按照配准获得的参数加载到600nm图像上,如图5中的轮廓线,该轮廓线就是570nm图像中分割获得的血管轮廓,通过加载到600nm图像中,观察600nm图像中的血管轮廓与轮廓线的重合程度,来判断配准结果的准确性.

如图5显示了直接互信息配准算法与本文算法在配准后570 nm图像中的血管轮廓线叠加到600 nm图像下的叠加效果.

通过图5中2幅血管轮廓的效果叠加图来看,图5(a)和(b)中分别画出了3组不同的区域,通过同种颜色框区域内容的对比,可以看出图5(b)中各个方框中的轮廓线和血管轮廓线的重合度要优于图5(a)中的各区域.通过直观观察的方式,可以证明本文配准方法优于直接互信息的配准方法.

本文方法充分考虑到不同波长视网膜图像背景的复杂性和灰度分布的不均匀性,充分利用图像中血管的轮廓信息和灰度信息,采用多窗口检测的方法对视网膜图像血管进行分割处理,再利用Powell优化算法的黄金分割搜索方法完成图像的配准.通过实验验证,本文方法相对于直接互信息的配准方法,可以很好地提高配准精确度,证明本文方法在不同波长视网膜图像的配准方法是可行的.但是本文还存在一定的问题,如配准过程中,分割的耗时相对要长一点,因此,在后续的研究中,分割方法的运行效率方面是研究的重点之一.

[1] HARRIS A, DINN R B, KAGEMANN L. A review of methods for human retinal oximetry[J]. Ophthalmic Surgery Lasers & Imaging,2003,34(2):152-164.

[2] 李新胜,刘宇. 基于边缘特征的虚实图像精确配准[J]. 四川大学学报(工程科学版),2014,46(5):96-103.

[3] 何敬,李永树,李歆,等. 基于点特征和边缘特征的无人机影像配准方法[J]. 西南交通大学学报(自然科学版),2012,46(6):955-961.

[4] 王红玉,冯筠,崔磊,等. 应用显著纹理特征的医学图像配准[J]. 光学精密工程,2015,23(9):2656-2665.

[5] PLUIM J, MAINTZ J, VIERGEVER M A. Image registration by maximization mutual information and gradient information[J]. IEEE Trans on Medical Imaging,2001,19(8):809-819.

[6] DEEPA M, SARAVANAN T. Automatic image registration using 2D-discrete wavelet transform[J]. Indian J Science & Technology,2016,9(5):1-3.

[7] ARAKI T, IKEDA N, DEY N, et al. A comparative approach of four different image registration techniques for quantitative assessment of coronary artery calcium lesions using intravascular ultrasound[J]. Computer Methods & Programs in Biomedicine,2015,118(2):158-172.

[8] 洪明坚,吕建斌,杨丹,等. 一种新的基于互信息的图像配准方法[J]. 重庆大学学报(自然科学版),2009,32(6):697-700.

[9] 翟利志,王敬东,李鹏. 基于邻域信息的红外与可见光图像互信息配准[J]. 计算机技术与发展,2008,18(10):151-154.

[10] 华梓铮,华泽玺. 结合NSCT高低频特征的图像边缘检测算法[J]. 四川师范大学学报(自然科学版),2014,37(4):612-617.

[11] 高向军. 基于多尺度线性检测的视网膜血管分割[J]. 科学技术与工程,2013,13(23):6820-6824.

[12] NGNYEN U T V, BHUIYAN A, PARK L A F. An effective retinal blood vessel segmentation method using multi-scale line detection[J]. Pattern Recognition,2013,13(3):703-715.

[13] YAMADA M, SIGAL L, RAPTIS M, et al. Cross-domain matching with squared-loss mutual information[J]. IEEE Transactions on Pattern Analysis & Machine Intelligence,2015,37(9):1764-76.

[14] 潘婷婷,陆丽婷,顾绮芳. QPSO算法和Powell法结合的多分辨率医学图像配准[J]. 计算机应用与软件,2014,31(7):237-240.

[15] 乔玉龙,赵源萌,张存林,等. 基于混合Powell法的太赫兹/可见光双波段图像配准[J]. 光学学报,2015,35(2):151-157.

[16] SHAMES R , SADEGHI P, KENNEDY R. Parallel computation of mutual information on the GPU with application to real-time registration of 3D medical images[J]. Computer Methods & Programs in Biomedicine,2010,99(2):133-46.

(编辑 李德华)

Registration of Retinal Vessel Blood Image Based on Image Segmentation

LIU Zhidong1,2, LIN Jiangli1, CHEN Ke1

( 1.CollegeofMaterialsScienceandEngineering,SichuanUniversity,Chengdu610065,Sichuan; 2.DepartmentofAutomobileandInformationEngineering,UrbanVocationalCollegeofSichuan,Chengdu610010,Sichuan)

In order to calculate the blood oxygen saturation of retinal images, different wavelength images should be registered. This paper presents an image registration method based on the segmentation of blood vessel image and mutual information. In the study, in order to reduce the impact of the information on the registration result, the registration of image segmentation, extraction of retinal vessels information in the image; calculating the vascular similarity in the segmented image, and using the Powell optimization algorithm 0.618 one-dimensional search algorithm to improve the speed of operation; the different wavelength of the image registration based on the calculated similarity value. In the study, the error average of parameters (angle,Xdirection,Ydirection) calculated from the registration algorithm is 2.00%, 2.53% and 2.52%, and the variance of the error is 0.57, 2.09 and 0.34, were better than the direct mutual information registration method. Experiments show that the algorithm can automatically and effectively register retinal images with different wavelengths, and has good repeatability and stability.

retinal image; image segmentation; registration; mutual information

2017-02-21

国家自然科学基金(81301286)和四川省科技支撑项目(2014GZ0005)

TP391.4

A

1001-8395(2017)04-0554-07

10.3969/j.issn.1001-8395.2017.04.020

*通信作者简介：陈 科(1982—),男,博士,主要从事医学图像处理的研究,E-mail:chenke@scu.edu.cn