张石, 佘黎煌, 王雅凡, 苏婷
东北大学 计算机科学与工程学院, 辽宁 沈阳 110169
收稿日期：2017-11-14
基金项目：中央高校基本科研业务费专项资金资助项目(N150403002)。
作者简介：张石(1963-), 男, 辽宁抚顺人, 东北大学教授, 博士生导师。

摘要：引入善于描述非稳定图像的多重去趋势波动理论, 提出一种基于多重去趋势分析的视网膜图像分割方法.该方法采用直方图均衡化对图像进行预处理来增强血管影像, 然后采用多重去趋势波动分析计算图像的广义赫斯特指数, 并利用血管指数特性来分割血管, 最后用形态学进行图像后处理, 得到最终的血管图像.基于DIARETDE0和DIARETDE1两个数据库进行实验.结果表明, 该方法在处理视网膜病变图像时有较好的完整性和连通性, 能够较好地提取血管主体, 具有很好的临床应用价值.
关键词：视网膜血管分割多重分形去趋势波动分析病变图像
Retinal Image Segmentation Based on Multifractal Detrended Fluctuation Analysis
ZHANG Shi, SHE Li-huang, WANG Ya-fan, SU Ting
School of Computer Science & Engineering, Northeastern University, Shenyang 110169, China
Corresponding author: SHE Li-huang, E-mail: shelihuang@ise.neu.edu.cn
Abstract: A method based on multifractal detrended fluctuation analysis(DFA), which is adaptable for unsteady images, is introduced to segment retina images. In this method, histogram equalization is first used to precondition images for the enhancement of blood vessels, and then, the generalized Hurst index of images is calculated by multifractal DFA. The blood vessel index is used to segment the blood vessels. Finally, morphological postprocessing is used to get the final blood vessel images. Experiments based on DIARETDE0 and DIARETDE1 databases show that the proposed method has better integrity and connectivity in the treatment of retinal lesions, and it can extract the main body of blood vessels better, which has a good clinical value.
Key words: retinasegmentation of blood vesselmultifractaldetrended fluctuation analysislesion image
眼底视网膜血管是人体血管部分中唯一可无创直接检测到的, 很多眼科疾病、心血管疾病和糖尿病的病变会直接反映在视网膜血管的网络结构上, 引起相应的结构变化.在临床上, 医生可以通过视网膜血管的直径、分支形态、角度等结构变化进行诊断, 也可以通过直接观察眼底图像病变来进行检测; 因此, 彩色视网膜图像血管的检测与提取对于疾病筛查与诊断具有重要意义^[1].
目前眼底血管分割方法有很多, 主要有：基于模式识别的方法, 其中分为监督分类^[2]和无监督分类^[3]; 基于匹配滤波的方法^[4]; 基于数学形态学的方法^[5]; 基于模型的方法^[6]; 基于血管追踪的方法^[7].
分形理论在图像领域中的应用已取得许多成果, 研究表明, 仅由一个分形维数来描述复杂的非线性系统是不够的, 在各个复杂形体形成过程中, 其局域条件十分重要.为了进一步了解在分形体形成过程中局域条件的作用, 文献[8]提出了多重分形方法.文献[9]提出了基于容度的多重分析法, 将图片中不同类型的纹理区分开来.文献[10]采用多重分形频谱分析的方法, 对原木CT图像进行边缘检测, 取得了很好的结果, 为边缘提取提供了新的方法.文献[11]利用多重分形理论分析了油菜病虫害叶片, 能够清晰地分割出叶片边缘轮廓, 准确地定位病变区域, 还能保留较多的细节.文献[12]应用多重分形方法对医学图像进行分割, 证明了方法是可行的, 对比较复杂的图像也有较好的分割效果.大量实验证明, 针对复杂图像, 多重分析分割算法优于传统分割算法, 具有良好分割结果的同时还有很好的抗噪性.
病变的视网膜图像中血管含有病变斑点, 血管局部严重变形, 导致图像的局部灰度发生剧烈变化, 灰度分布不稳定, 形成某种局部趋势, 现有的分割方法不能很好地分割出主体血管.对于非平稳对象, 首先应将某种局部趋势去除, 然后再对剩余序列也就是矩阵进行分析, 本文提出一种针对病变图像的视网膜眼底血管分割方法, 利用基于多重分形维数的去趋势波动分析来处理病变图像, 提取病变图像血管主体, 相比其他方法取得了较好的效果.
1 图像预处理对数据库内的视网膜眼底图像提取绿色通道图像, 对绿色通道图像进行直方图均衡化增强处理, 有效提高图像的对比度, 得到了更多的图像细节, 有利于下一步血管的提取(见图 1).
图 1(Fig. 1)

图 1 视网膜预处理结果Fig.1 Retinal preconditioning results

2 多重分形去趋势波动分析方法一个二维图像可以离散成一个M×N阶矩阵X=(X_{i, j}), (i=1, 2, …, M; j=1, 2, …, N).首先将图像划分为等长度s的非重叠子区域, 大小为M_s×N_s, 每个子区域定义为X_{m, n}:

(1)

其中1≤i, j≤s, r=(m－1)s和t=(n－1)s, 每个子区域表面像素累积和为

(2)

然后用平面来拟合每个子区域G_{m, n}(i, j)的趋势, 并确定残差矩阵中的元素g_{m, n}(i, j), 通过残差矩阵可以求出每个子区域的去趋势函数F²(m, n, s):

(3)

已知每个子区域的去趋势函数, 计算x阶波动函数:

(4)

当x=0时, 根据洛必达法则, 由

(5)

可以求出F₀(s).
最后通过s取值范围的变化(6≤s≤min(M, N)), 得到一组波动函数和尺度s的幂值关系:

(6)

式中h(x)是广义赫斯特指数, 当t=2时是著名的表面赫斯特指数.h(x)描述图像特征, 根据不同的h(x)值, 提取不同区域, 据此分割图像, 并得到每个像素点的图像特征.以每个像素点(i, j)为中心设置大小为a×a的移动窗口, 求其赫斯特指数, 得到L_h(x).图像广义赫斯特指数表示图像局部的非平稳特性.
3 视网膜图像分割算法本算法基于多重分形去趋势分析的思想并加以改进.每个图像可以看作256个灰度级的二维矩阵, 可以用上述算法来计算视网膜图像h(x), 根据血管特征来分割图像.算法步骤如下.
1) 研究表明眼底血管具有自相似性, 因此眼底图像可以作为自相似表面.设定大小为a×a移动窗口, 遍历整幅图片, 对每个像素点I(i, j)按照上述算法求出广义赫斯特指数L_h(x), 找到最大值L_{h(x), max}和最小值L_{h(x), min}, 确定间断范围[L_{h(x), min}, L_{h(x), max}].移动窗口的尺寸会严重影响算法的准确性:窗口过大, 错失细节; 过小, 则没有足够的点来拟合曲线.通过反复试验, 将窗口尺寸设置为7×7.
2) 将间断范围均分为n段, 经过大量试验, 选取n=40.
3) 在每个间断中使用简单的盒子维数计算由L_h(x)组成的图像, 即在[L_{h(x), min}, L_{h(x), max}]内的子图像; 用大小为λ×λ的框覆盖子图像, 当框内具有段内的像素值时, 记录框的数量, 盒子遍历整个图像, 记录框的个数N(λ).
4) 通过尺度λ的变化(λ=2, 4, 6, 8, 16), 根据公式, 可以求出由L_h(x)构成图像的盒子维数D.
5) 重复步骤3)和步骤4), 得到一系列的D(L_h(x))：D₁(L_h(x)), D₂(L_h(x)), …, D_n(L_h(x)).
D(L_h(x))作为L_h(x)的谱函数, 表示图像的全局奇异信息, 可以根据D(L_h(x))的数值确定图像的奇异性并作为分割标准来分割图像.如果D(L_h(x))数值接近1, 说明相应的像素点为平滑边界点; 如果数值接近2, 则相应的像素点位于光滑表面.通过对数值结果的判断, 分析图像的奇异性, 设定D(L_h(x))范围来确定图像的奇异区域和奇异边界, 本文为了提取出主体血管, 经过多次试验, 将D(L_h(x))范围设定为[1.7~1.9].
4 形态学后处理对处理后的图像进行形态学闭运算：

(7)

其中I为处理后的图像, StrucElem选取单位为1的圆形结构元素.去除图像孤立的点, 得到最终结果.
5 实验结果与分析5.1 稳定性测试为了体现本文算法的稳定性, 选取带有手动分割结果的DRIVE数据库进行实验, 定量分析算法的鲁棒性, 并与其他算法进行比较.通常使用三个统计测量指标来评价算法结果, 分别是敏感度(TPR)、特殊性(FPR)和正确性(ACC).TPR表示血管被正确分类为血管的标准化测量, FPR度量非血管点被正确分类的比例, ACC则是每一个像素点被正确分类的情况.从理论上来说, 三个指标的计算结果越高越好, 但是FPR和TPR是一对矛盾的指标, 算法旨在寻求两者之间的平衡.
三个测量标准对应的计算公式如下:

(8)

(9)

(10)

式中各符号含义见表 1.
表 1(Table 1)

表 1 分割符号描述Table 1 Description of segmentation symbols 分割结果 正确分割的像素点 错误分割的像素点 血管点 TP FP 背景点 TN FN 总计 P N

表 1 分割符号描述 Table 1 Description of segmentation symbols

实验在视网膜图像中加入不同程度椒盐噪声和不同程度的乘性噪声, 噪声参数分别为0，0.05, 0.1, 0.15.对图像进行同样的预处理, 用不同的算法进行视网膜血管提取.对比算法分别是COSFIRE算法^[4], 形态学^[13]和Hessian矩阵算法, 实验结果如图 2所示.
图 2(Fig. 2)

图 2 采用不同算法的图像分割性能Fig.2 Performance of image segmentation using different algorithms

综合两组折线图可以看出:初始时本文算法的分割精度稍低于对比算法, 但是分割结果也较为出色, 随着噪声参数值的增加, 算法的准确度、灵敏度和特异度都在降低; 但是在不同类型噪声的干扰下, 本文算法受噪声的影响较小, 误差范围在理想范围内, 波动范围相对于其他算法较为平稳, 说明本文算法具有较强的稳定性, 抗噪能力强, 善于描述非稳定图像.
5.2 病变图像的分割实验为了验证算法对病变图像有着良好的分割效果, 从公开用于糖尿病视网膜病变检查的标准数据库DIARETDE0和DIARETDE1中选取含有糖尿病性病变和微血管瘤病变的图像.DIARETDE0数据库中包含130张眼底图像, 其中20张为健康的眼底图像, 其余110张为不同程度的糖尿病视网膜病变图像.DIARETDE1数据库包含89张眼底图像, 其中84张图像含有微血管瘤病变.
图 3为本文算法分割结果:其中第一行为本文算法在DIARETDE0数据库中的分割结果, 第二行为本文算法在DIARETDE1数据库中的分割结果.从图 3可以看出, 本文提出的基于多重分形去趋势波动分析的算法较好地分割出了病变图像主体血管, 血管网络较为完整.
图 3(Fig. 3)

图 3 病变图像分割结果Fig.3 Segmentation results of lesion images

为了验证本文算法分割病变图像的优势, 选取DIARETDE0和DIARETDE1中的图像进行处理, 并将结果与文献[13](基于形态学算法)和文献[4](COSFIRE算法)的视网膜分割算法的处理结果进行比较.如图 4、图 5所示.
图 4(Fig. 4)

图 4 不同算法对DIARETDE0数据库中图像的分割结果Fig.4 Segmentation results of the images from DIARETDE0 database using different algorithms (a)—原始图像；(b)—基于形态学算法；(c)—COSFIRE算法；(d)—本文算法.

图 5(Fig. 5)

图 5 不同算法对DIARETDE1数据库中图像的分割结果Fig.5 Segmentation results of the images from DIARETDE1 database using different algorithms (a)—原始图像；(b)—基于形态学算法；(c)—COSFIRE算法；(d)—本文算法.

从图 4和图 5中可以清晰地看出:经过相同的预处理步骤, 基于形态学的算法分割结果受病变影响严重, 过多地分割出血管, 将病变区域作为血管分割出来, 血管主体不清晰; COSFIRE算法的图像分割结果受病变区域影响较为严重, 过少地分割出血管, 丢失很多细节和血管; 而本文算法处理病变图像时能很好地保持稳定, 有效提取血管主体, 分割较为清晰, 能显示复杂图像的细节部分, 说明了算法的准确性、通用性和临床实用性.
6 结语本文采用多重分形去趋势分析来处理视网膜病变图像, 利用广义赫斯特指数作为指标将血管从背景中提取出来, 并针对DIARETDE0和DIARETDE1两个数据库中具有不同病变的图像进行了实验, 同时将本文算法与其他方法做了比较.可以看出, 本文算法处理病变图像时, 能够较为准确地分割出血管轮廓, 满足对血管分割完整性的要求; 同时算法的抗噪性比较好, 体现了算法的鲁棒性.但是本文算法的迭代过程较为耗时.今后的研究重点是将本文算法与遗传算法、神经网络等其他优化算法结合起来, 提高算法的运行速度与精度.
参考文献

[1]	Li Q, You J, Zhang D, et al. Vessel segmentation and width estimation in retinal images using multiscale production of matched filter responses[J].Expert Systems with Applications, 2011, 39(9): 7600–7610.
[2]	Tolias Y, Panas S. A fuzzy vessel tracking algorithm for retinal images based on fuzzy clustering[J].IEEE Transactions on Medical Imaging, 1998, 17(2): 263–273.DOI:10.1109/42.700738
[3]	Zhu T. Fourier cross-sectional profile for vessel detection on retinal images[J].Computerized Medical Imaging and Graphics, 2010, 34(3): 203–212.DOI:10.1016/j.compmedimag.2009.09.004
[4]	Azzopardi G, Strisciuglio N, Vento M, et al. Trainable COSFIRE filters for vessel delineation with application to retinal images[J].Medical Image Analysis, 2015, 19(1): 46–57.DOI:10.1016/j.media.2014.08.002
[5]	Djaroudib K, Ahmed A T, Zidani A. Textural approach for mass abnormality segmentation in mammographic images[J].International Journal of Computer Science Issues, 2013, 10(6): 125–131.
[6]	Sun K, Chen Z, Jiang S, et al. Morphological multiscale enhancement, fuzzy filter and watershed for vascular tree extraction in angiogram[J].Journal of Medical Systems, 2011, 35(5): 811–824.DOI:10.1007/s10916-010-9466-3
[7]	Yin Y, Adel M, Bourennane S. Retinal vessel segmentation using a probabilistic tracking method[J].Pattern Recognition, 2012, 45(4): 1235–1244.DOI:10.1016/j.patcog.2011.09.019
[8]	Grassberger P. An optimized box-assisted algorithm for fractal dimensions[J].Physica Letters:A, 1990, 148(1/2): 63–68.
[9]	李会方, 俞卞章. 一种基于多重分形新特征的图像分割算法[J].光学精密工程, 2003, 11(6): 627–631. ( Li Hui-fang, Yu Bian-zhang. Image segmentation approach based on new multifractal feature vectors[J].Optics and Precision Engineering, 2003, 11(6): 627–631.DOI:10.3321/j.issn:1004-924X.2003.06.017)
[10]	韩书霞, 戚大伟, 于雷. 基于多重分形理论的原木CT腐朽图像分析与处理[J].森林工程, 2007, 23(5): 15–18. ( Han Shu-xia, Qi Da-wei, Yu Lei. Analysis and processing of decayed log CT image based on multifractal theory[J].Forest Engineering, 2007, 23(5): 15–18.DOI:10.3969/j.issn.1001-005X.2007.05.005)
[11]	施文, 邹锐标, 王访, 等. 基于多重分形的油菜病虫害叶片图像分割[J].湖南农业大学学报(自然科学版), 2014, 40(5): 556–560. ( Shi Wen, Zou Rui-biao, Wang Fang, et al. Multifractal theory for image segmentation of rapeseed leaf affected by diseases and insects[J].Journal of Hunan Agricultural University(Natural Sciences), 2014, 40(5): 556–560.)
[12]	金春兰, 黄华, 刘圹彬. 基于多重分形的医学图像分割方法[J].中国组织工程研究与临床康复, 2010, 14(9): 1535–1538. ( Jin Chun-lan, Huang Hua, Liu Kuang-bin. Medical image segmentation based on multifractal theory[J].Journal of Clinical Rehabilitative Tissue Engineering Research, 2010, 14(9): 1535–1538.DOI:10.3969/j.issn.1673-8225.2010.09.004)
[13]	Heneghan C, Flynn J, O'Keefe M, et al. Characterization of changes in blood vessel width and tortuosity in retinopathy of prematurity using image analysis[J].Medical Image Analysis, 2002, 6(4): 407–429.DOI:10.1016/S1361-8415(02)00058-0