光学相干层析成像的视网膜层状结构自动分割

高用贺，李跃杰,，王立伟，张明蓉

1 北京协和医学院 中国医学科学院生物医学工程研究所，天津市，300192

2 天津市眼科医学设备技术工程中心，天津市，300384

光学相干层析成像的视网膜层状结构自动分割

【作 者】高用贺1，李跃杰1,2，王立伟1，张明蓉1

1 北京协和医学院 中国医学科学院生物医学工程研究所，天津市，300192

2 天津市眼科医学设备技术工程中心，天津市，300384

目的 利用算法实现对视网膜层状结构的自动分割及定量分析是光学相干层析成像技术应用于青光眼及视网膜病变早期诊断的关键。现存处理方法对图像质量要求较高且可靠性不高。该文拟利用改进的非线性复合扩散滤波等方法解决这个问题。方法 首先对自主搭建的OCT系统获得的20幅视网膜图像，通过自动阈值、改进的非线性复合扩散滤波、形态学操作、峰值探测等综合算法，进行分割，比较准确的分割出内界膜（ILM）、外核层（ONL）、内节层和外节层（IS/OS）以及视网膜色素上皮与脉络膜层（RPE_ChCap）边界，最后测量得到视网膜的厚度。结果 本算法对视网膜的分割与专家手动测量有较好的一致性，视网膜中心凹测量结果与Zeiss Stratus OCT视网膜中心厚度212±20 μm数据一致。结论该文提出的算法有希望应用于临床视网膜疾病的诊断。

视网膜；光学相干层析；层状结构；非线性复扩散滤波；形态学操作

0　引言

光学相干层析成像技术(Optical Coherence Tomography，OCT) 是近年来迅速发展起来的一种成像技术，最早由美国麻省理工学院的 Huang D等在1991年成功应用于眼组织[1-2]。研究证实，黄斑中心凹光感受器层厚度随青光眼病程进展呈现出动态、非线性变化的趋势：在青光眼早期明显增厚，而在青光眼中晚期又逐渐下降，厚度接近正常水平[3]。同时，黄斑区视网膜厚度也是诊断及监测糖尿病黄斑水肿遗传性黄斑病变等眼科疾病的重要依据[4]。对这些有重要意义的数据的测量则依赖于对光学相干层析图像中层状结构的精确分割。因此，视网膜层状结构的定量分析成为眼科光学相干层析成像的一项核心技术。在此领域国外机构起步早，2001年，Koozekanani D等[5]发表了视网膜分割的文献，此后每年都有很多关于视网膜层状结构分割的文献发表，并在视网膜层状结构的分割层数与精确度上逐步提升[6-8]。国内OCT成像研究领域进展相对较慢，在对OCT视网膜层状结构的分割领域研究比较少，可检索到的中文文献很少，水平相对国外仍有一定差距[9]。这也严重制约了OCT眼科光学相干层析成像仪国产化的步伐。现有的视网膜层状结构分割算法，大多基于多次帧平均(即对眼底某个位置截面扫描多次)后的图像进行分割。图像噪声大幅减少，虽然分割结果比较理想，但在OCT仪器在对人眼扫描过程中，由于很多原因例如人眼总会有不自主抖动，多次扫描经常会有伪差，大大影响算法分割的准确性。本文提出的视网膜层状结构的自动算法，是基于一帧OCT视网膜A扫描图像基础上，先用改进的自适应非线性复合扩散滤波器[10]减少噪声影响，再使用形态学操作等方法进行分割。用自主研发的OCT平台获得的视网膜图像进行算法测试，对视网膜厚度的测量及层状结构的定量分析取得了良好的结果。

1　视网膜层状结构分割

利用天津市眼科医学设备技术工程中心所搭建的高分辨率SD-OCT样机进行眼底视网膜扫描。医学上，视网膜结构可分为七层，分别为：视网膜神经纤维层(retinal never fiber layer, RNFL)、神经节细胞层和内网状层(inner plexiform layer, IPL)、内核层(inner nuclear layer, INL)、外网状层(outer plexiform layer, OPL)、光感受器内节/外节(inner and outer photoreceptor segment, IS/OS)、视网膜色素上皮与脉络膜层(retinal pigment epithelium and choroid, RPE_ ChCap)，如图1所示。

图1　OCT视网膜图分层结构Fig.1　Layers structures of retina picture.

1.1内界膜(Internal limiting membrane, ILM)定位

首先，对原OCT视网膜图像使用最大类间方差法(简称OTSU)找到一个合适的阈值(threshold)，这个阈值在[0, 1]范围内。再利用这个阈值，用Matlab自带的im2bw函数将灰度图像转换为二值图像。此方法比人为手动设定的阈值效果更好且方便。对二值图像水平方向上(即0o)实施线形膨胀运算，结构元素参数设为4，这样可以将ILM边界填充的更平滑。遍历寻找图像每一列值为1的纵坐标最小的位置，即ILM初始边界。对初始边界第6个点开始，每一个点与其前5个点坐标均值作比较，若差的绝对值大于5，则用此坐标均值代替此点的坐标，以达到平滑边界的目的，最终处理结果为ILM边界，如图2所示。

图2　ILM（内界膜）层边界Fig.2　Inner Limiting Membrane

1.2RPE_Ch和IS/OS边界定位

1.2.1图像预处理

由于OCT系统中，目标的光学性能和移动、光源的光斑的大小和时间相干性、传输光的多次散射和相位畸变等都是造成OCT图像散斑噪声的主要因素。为了对OCT视网膜图像更准确的层分割，显然在保持图像层状信息的基础上，去除影响图像质量的散斑噪声，显得尤为重要。

目前，已有一些算法，如中值滤波，高斯滤波，小波非线性阈值滤波，零幅度算法(Zero-adjustment procedure, ZAP)[11]，清洁算法(CLEAN)[12-13]，增强的Lee滤波器(ELEE)等方法。本文采用一种改进的自适应非线性复扩散滤波器，是由Rui Bernardes于2010年提出的，它基于普通的非线性复合扩散滤波器。相对于其他算法，此算法对OCT视网膜图像具有更好的噪声消除和层状结构边缘保持能力。

一般的非线性各向异性复合扩散是基于式(1)的解：

式中，△●是散度，△是梯度，D是扩散系数。

将式（1）时域前向离散并化为空域有限差分式，即

式中，△和△分别为离散拉普拉斯和梯度算

hh子，是迭代步数n的时间，i, j, m是体数据块的坐标向量，在2D数据中

成立，在（4）式中我们令△h=△y=1且对于2D图像数据α=4。Salinas将扩散系数定义为

式中，i =√-1，k是阈值参数，θ是一个接近于0的相位角。此算法中，扩散使得图像平滑但衰减了层状结构边缘信息。

扩散系数可以简化为

基于OCT视网膜图像的背景（主要是玻璃体、房水）区域强度较低，而视网膜区域强度相对较高的基础上，在低强度区域促进扩散，而在高强度区域降低扩散，这样能更好的保留层状结构边缘信息。因此，要改进参数k，使其在图像低强度区域增强，在图像高强度区域减小，如下式

式中，g=GN,σ*Re(I)，GN,σ是高斯平均内核，*是卷积运算，尺寸为N×N或N×N×N，σ为标准差。

相对于其他非线性复合扩散算法，大多设置的是固定的时间步长(△t)，由于扩散系数依赖于图像的梯度，并且在扩散迭代初始时因为噪声的缘故使梯度值更大，故这里设置为自适应时间步长，即

式中，α=4，且a+b≤1。

将扩散时间最大值设为0.8 s，对OCT视网膜图像进行处理，之后截取ILM层上部一块儿矩形背景区域，求取平均像素强度Ibg，用下式对OCT视网膜图像进行局部增强，即

式中，I为整体OCT视网膜图像，Ibg是所截取部分的平均强度，这里β=2，重复连续进行3次。以达到增强RPE部分的目的，结果如图4所示。

1.2.2形态学分割RPE_ChCap

二值形态学中的运算对象是集合，通常给出一个图像集合和一个结构元素集合，利用结构元素对图像进行操作。其中，结构元素是一个用来定义形态操作中所用到的形状和大小的矩阵，该矩阵仅由0和1组成，可以具有任意的大小和维数，数值1代表邻域内的像数，形态学运算都是对数值为1的区域进行运算。

图4　增强RPE结果Fig.4　Results of enhanced RPE

形态学开运算，即对其先腐蚀后膨胀。可消除细小物体、在纤细处分离物体、平滑较大物体的边界时不明显改变其面积。

首先，同ILM层边界定位中一样，将图4二值化，用Matlab中bwareaopen函数去掉面积70以下区域，结果如图5所示。

图5　二值化结果Fig.5　Results of binarization

接下来开运算操作，为了使OCT视网膜层状结构边缘更加准确，腐蚀操作结构元素选用半径为2的圆盘状结构；膨胀操作结构元素选用长度为20，方向为水平（即0o）的线形结构，用Matlab中bwareaopen函数去掉面积5 000以下区域（即RNFL区域）。遍历寻找图像每一列值为1的纵坐标最小、最大的位置，即分别为IS/OS边界、RPE_Ch边界。为了减小IS/OS层边界左右两端由于线型膨胀操作所引起的误差，以RPE_Ch边界最低点为基准，将图像逐列下移，使RPE_Ch边界拉平为水平直线。

对拉平的图像以相同设置重复如上操作，即对RPE部分局部增强、形态学开运算分割，并边界定位，结果如图6所示。

图6　拉平原图后RPE层分割结果Fig.6　Results of segmentation of RPE after alignment

1.2.3分割ONL

对经过预处理过的拉平的OCT视网膜图像每11列像素宽的A扫描叠加求均值赋于第1列，以第200列像素为例，画出其灰度强度图。最大的峰值出现在RPE_ChCap层内部，从最大峰值处向上使用峰值探测法寻找，第一个峰值即是ONL层内边界。之后，遍历图像所有列，用上述方法找到所有峰值点。将所有峰值点连成光滑曲线，即为ONL层边界，如图7所示，由上往下第二条分割线即为ONL层边界。

图7　ONL层位置Fig.7　Location of ONL

2　视网膜厚度测量

将ILM边界距离RPE_ChCap边界的两个极大值分别定义为旁中心凹（parafovea），在旁中心凹中间寻找ILM边界的最低点并定义为中心凹（fovea）如图8所示。提取ILM边界与IS/OS边界的像素差，再经过定标后生成视网膜层厚度分布图，如图9所示。对所有列厚度加权求平均值，即为此A扫描的视网膜平均厚度。视网膜光感受器层厚度是以中心凹处ILM边界到IS/OS边界处的距离表示，如图8箭头长度所示。

图8　中心凹及旁中心凹位置Fig.8　Location of fovea and parafove

本文在检验算法的准确性及可靠性的过程中，用自主搭建的设备采集20个正常人眼底视网膜OCT图像。运用本文所提出的算法进行自动测量黄斑中心凹光感受器厚度和视网膜平均厚度。并与专家手动分割测量结果进行比较，结果如表1所示。

由表1可知，本文提出的自动算法与专家手动分割测量结果无显著性差异，具有较好的一致性。其中，视网膜中心凹厚度与Zeiss Stratus OCT视网膜中心厚度212±20 μm数据一致[14]。之所以自动测量的结果较手动测量偏低，是由于形态学操作的误差所致，但差异范围在3%以内。

图9　视网膜厚度分布曲线Fig.9　Retinal thickness distribution curve

表1　黄斑中心凹光感受器厚度和视网膜平均厚度(n=20) (μm)Tab.1　Photoreceptor thickness of macular fovea and mean retinal thickness(n=20)

3　结论

光学相干层析系统具有非接触、无损伤、实时、高分辨力等特性，对光学相干层析图像中的层状结构进行定量提取与其他检测方法相比具有很多优点。本文提出的算法，首先通过使用最大类间方差法设定阈值结合线型膨胀运算对ILM层边界定位。其次，通过改进的自适应非线性复合扩散滤波器对OCT视网膜图像进行预处理，既滤除图像中的噪声，同时又保留了更多的图像层状特征，减少了层状结构位置失真。对OCT视网膜图像的预处理，使得背景区域灰度更加均匀集中，对RPE_ChCap区域进行线性增强奠定了基础，进而对IS/OS边界和RPE_ChCap边界进行定位。最后，利用峰值探测法对ONL层进行边界定位。结果显示，此方法具有较好的自适应性和重复性，可以应用在临床OCT系统中。

当然，此方法也有需要改进的地方，还不能对视网膜全部层状结构进行分割。对较严重眼底病变的视网膜图像自动分层的准确度有待提高。今后，可在此算法基础上，利用马尔可夫随机场模型对视网膜OCT图像的全部层状结构进行分割，并添加模式识别，以实现视网膜病变的全自动识别诊断。

[1] Huang D, Swanson EA, Lin CP, et al. Optical coherence tomography[J]. Science, 1991, 254(5035): 1178-1181.

[2] Drexler W, Fujimoto JG. State-of-the-art retinal optical coherence tomography[J]. Prog Retin Eye Res, 2008, 27(1): 45-48.

[3] Quigley HA, Nickells RW, Kerrigan LA, et a1．Retinal ganglion cell death in experimental glaucoma and after axotomy occurs by apoptosis[J]. Invest Ophthalmol Vis Sci, 1995, 36: 774-786．

[4] Ishikawa H, Stein DM, Wollstein G, et al. Macular Segmentation with Optical Coherence Tomography[J]. Invest Ophthalmol Vis Sci, 2005, 46(6): 2012-2017.

[5] Koozekanani D, Boyer K, Roberts C. Retinal thickness measurements from optical coherence tomography using a Markov boundary model[J]. IEEE Trans Med Imag, 2001, 20(9): 900-916.

[6] Tan O, Chopra V, Lu ATH, et al. Detection of macular ganglion cell loss in glaucoma by Fourier-domain optical coherence tomography[J]. Ophthalmology, 2009, 116(12), 2305-2314.

[7] Yazdanpanah A, Hamarneh G, Smith B, et al. Intra-retinal layer segmentation in optical coherence tomography using an active contour approach[C]. Proc MICCAI, 2009: 649-656.

[8] Kajic V, Povazay B, Hermann B, et al. Robust segmentation of intraretinal layers in the normal human fovea using a novel statistical model based on texture and shape analysis[J]. Opt Express, 2010, 18(14): 14730-14744.

[9] 江源源, 周传清, 任秋实. 基于光学相干层析的视网膜图像分割[J]. 北京生物医学工程, 2011, 30(5): 452-456.

[10] Bernardes R, Maduro C, Serranho P, et al. Improved adaptive complex diffusion despeckling fi lter[J]. Opt Express, 2010, 18(23): 24048-24059.

[11] Yung KM, Lee SL, Schmitt JM. Phase-domain processing of optical coherence tomography images[J]. J Biomed Opt, 1999, 4(1): 125-136.

[12] Fried DL. Analysis of the clean algorithm and implications for super resolution[J]. J Opt Soc Am A,1995, 12(5): 853-860.

[13] Schmitt JM. Restoration of optical coherence images of living tissue using the clean algorithm[J]. J Biomed Opt, 1998, 3(1): 66-75.

[14] Legarreta JE, Gregori G, Punjani OS, et al. Macular thickness measurements in normal eyes using spectral domain optical coherence tomography[J]. Ophthalmic Surg Las Imag, 2008, 39(4 Suppl): S43-49.

Automated Segmentation of Retina Layer Structures on Optical Coherence Tomography

【 Writers 】Gao Yonghe1, Li Yuejie1,2, Wang Liwei1, Zhang Mingrong1
1 Peking Union Medical College&Chinese Academy of Medical Sciences Institute of Biomedical Engineering, Tianjin, 300192
2 Tianjin Ophthalmic Medical Device Technology Engineering Center, Tianjin, 300384

retina, optical coherence tomography, layer structure, nonlinear complex diffusion fi ltering, morphological operations

TP391.41

A

10.3969/j.issn.1671-7104.2014.02.004

1671-7104(2014)02-0094-04

2014-01-02

中国医学科学院生物医学工程研究所院所基金（1312）；天津科技创新体系及平台建设计划项目（13ZXGCCX06300）

高用贺，研究生，E-mail: gyhklgyhkl@126.com

李跃杰，研究员，E-mail: liyj_1@sina.com

【 Abstract 】Objective Using the algorithm on the layered structure of the retina and quantitative analysis of the automatic segmentation technique is the key to the early diagnosis of glaucoma and other retinopathy on optical coherence tomography. Existing methods require high qulity image and have low reliability. This paper used the improved complex nonlinear diffuse fi ltering and other methods to solve this problem. Methods This paper includes algorithm such as automatic threshold, improved complex nonlinear diffusion filtering, morphological operations and peak detection. Use the method for the segmentation of 20 retinal layers images which acquired on the self-builded OCT system, the boundary of inner limiting membrane(ILM), outer nuclear layer(ONL), the photoreceptor segments(IS/ OS) and the RPE_ChCap layer are detected accurately. At last, the photoreceptor layer thickness is measured. Results The results of segmentation and measurement are good corresponded with expert manual segmentation and measurements, retinal foveal measurements data is consistent with Zeiss Stratus OCT central retinal thickness 212±20 μm. Conclusion The algorithm proposed is prospective applied to clinical diagnosis of retinal diseases.