基于HD-OCT视网膜图像的脉络膜层自动分割方法

摘要

本发明公开了一种基于高分辨率频域光学相干断层（HD-OCT）图像的脉络膜层自动分割方法，属于图像处理技术领域。该方法首先对输入的HD-OCT图像进行去噪预处理，并通过定位内界膜去除视网膜神经纤维层附近的高反射率区域，然后根据高反射率信息定位视网膜色素上皮层的下边界，即脉络膜层的上边界，最后利用图搜索方法将利用脉络膜下边界的图像特性得到的候选CSI边界点相连，就得到脉络膜的CSI边界。实验结果表明，本发明中所得到的脉络膜层分割精度较高，与手动分割结果相当，可以代替临床医生手动测量脉络膜层厚度的繁琐耗时工作，对提高医生的工作效率具有重要意义。

说明书

技术领域

本发明涉及一种图像分割处理的方法，特别是一种基于高分辨率频域光学相 干断层（HD-OCT）视网膜图像的脉络膜层自动分割方法。

背景技术

HD-OCT视网膜图像是一种高分辨率的频域光学相干断层成像图像，它可以 有效地呈现脉络膜层的脉络膜巩膜界面（CSI），临床实验表明HD-OCT图像与 增强深度成像光学相干断层成像（EDI-OCT）图像都能够用于测量脉络膜层的厚 度。传统的视网膜层分割方法不适用于脉络膜层的分割，因为脉络膜层的CSI 边界在图像上表现较弱、且由于血管的影响存在边界断裂的现象。最近两年主要 出现了如下三种脉络膜层分割方法：

（1）基于纹理和形状信息的两阶段统计模型。该方法是针对1060nm OCT 系统的分割方法，且需要大量的训练样本。

（2）基于梯度和图论的分割方法。该方法是针对EDI-OCT图像的，在 EDI-OCT图像中，由于成像焦距靠近脉络膜层，所以视网膜色素上皮层（RPE） 的反射率相对来说是最高的，不存在视网膜神经纤维层（RNFL）的干扰问题， 因此通过反射率就可以容易地得到RPE层。该方法中CSI候选点是通过寻找局 部灰度谷底得到的，脉络膜层中的血管影响会导致很多错误的CSI候选点，这对 后续的图论分割造成了一定的困难。

（3）基于图论的多阶段分割方法。该方法是一种半自动的分割方法，分割 的对象是一般的频域光学相干断层成像（SD-OCT）图像。

现有的脉络膜层分割方法不适用于HD-OCT图像，因为不同的成像图像具 有不同的成像特性，如EDI-OCT图像中RNFL层的反射率明显低于RPE层，而 HD-OCT图像中两者的反射率相当。

发明内容

本发明的目的在于提供一种基于HD-OCT视网膜图像的自动脉络膜层分割 方法。

实现本发明的目的的技术解决方案为：一种基于HD-OCT视网膜图像的自 动脉络膜层分割方法，包括以下步骤：

步骤1、采集HD-OCT视网膜图像；

步骤2、采用改进的双边滤波算法对输入图像进行去噪处理；具体是将传统 的双边滤波的各项同性高斯邻域窗口改为各项异性的高斯邻域窗口，其中，传统 的双边滤波算法的公式为：

$h\left(x\right)=k^{-1}\left(x\right){\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f\left(\xi \right)c(\xi ,x)s(f\left(\xi \right),f\left(x\right))\mathrm{d\xi }$

式中f和h分别为输入和输出图像，函数c(ξ,x)用于测量邻域中心点x和邻 域点ξ之间的空间距离，函数s用于测量两点间的灰度相似性，函数c和函数s都 是高斯函数，$k\left(x\right)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }c(\xi ,x)s(f\left(\xi \right),f\left(x\right))\mathrm{d\xi }$是归一化函数。

步骤3、根据玻璃体和视网膜的反射率差异定位内界膜（ILM）层；具体为： 通过阈值得到玻璃体区域，然后寻找玻璃体的最下方边界点从而得到ILM层， 所述阈值为41。

步骤4、去除与ILM层相近的高反射率RNFL层；具体为：

步骤4-1、生成模板图像M，其中ILM层像素置为1，其余像素置为0；

步骤4-2、用数学形态学对模板图像进行膨胀，所用公式为：

Ms=M⊕S

式中‘⊕’表示膨胀算子，S为半径30个像素的圆盘形结构元；

步骤4-3、将与Ms目标区域相连的高反射率RNFL层进行去除。

步骤5、利用高反射率和视网膜的结构特性估计RPE层；具体步骤为：

步骤5-1、将RNFL层以外的高反射率区域置为候选RPE层区域；

步骤5-2、根据RPE层的厚度约束去除厚度大于20个像素的虚假RPE区域；

步骤5-3、将每列的候选RPE区域的中心点作为RPE层的中轴点，然后采 用四阶多项式拟合得到RPE层的中轴线。

步骤6、在RPE层的下边界附近根据垂直梯度算法得到候选的布鲁赫膜 （BM）边界点，然后采用四阶多项式拟合得到最终的BM；RPE层的下边界附 近为RPE层中轴线下方的20个像素窄带区域内；候选的BM边界点的选取方法 为：

将每列中在此窄带区域内的最大垂直梯度点作为候选的BM边界点，然后采 用四阶多项式拟合得到最终的BM。

步骤7、基于步骤6得到的BM拉平输入图像，然后将CSI分割限定在BM 下方的200个像素厚的窄带区域；

步骤8、根据CSI下方像素的反射率渐变特性得到CSI的候选边界点；具体 步骤为：

步骤8-1、生成灰度渐变距离图像D

$D(x,y)=\left(\begin{array}{ccc}D(x+1,y)+1& \mathrm{if}& I_a(x,y)\ge I_a(x+1,y)\\ 0& \mathrm{if}& I_a(x,y)<I_a(x+1,y)\end{array}\right)$

式中(x,y)表示图像像素坐标，$I_a=\left(\begin{array}{c}1\\ -1\end{array}\right)*I$表示垂直方向的灰度差，I为平滑后 的拉平脉络膜区域图像，‘*’表示卷积算子；

步骤8-2、通过灰度渐变距离图像D生成灰度渐变距离图像的垂直差图像D_x

$D_x=\left(\begin{array}{c}-1\\ 1\end{array}\right)*D$

步骤8-3、将每列垂直差图像D_x中的最大值位置作为CSI边界的候选点。

步骤9、通过图搜索在候选CSI边界点中寻找初步的CSI边界，然后采用四 阶多项式拟合得到最终的CSI。具体步骤为：

步骤9-1、构造图G=(V,E)，其中候选CSI边界点作为图G中的顶点V， 每个顶点与后续的k列顶点相连，E为图G中边的集合，图中任意两个顶点a和 b之间的边权重通过下式计算：

$w(a,b)=w(\mathrm{\Delta x},\mathrm{\Delta y},w_M,T_p,\alpha )=(\Delta x^2+\Delta y^2)+\frac{w_{M}H(\mathrm{\Delta y}-T_p)|\mathrm{\Delta y}-T_p|}{1+\exp (-\alpha (\mathrm{\Delta y}-T_p))}$

式中Δx和Δy分别为两顶点的水平和垂直距离，H(·)为Heaviside函数，惩罚因 子w_M为2000，T_p为阈值5，当Δy≥T_p时，边将被赋予额外的很大的惩罚权重， 参数α控制惩罚权重的提升速率，α的取值为2；

步骤9-2、采用Dijkstra算法搜素图G的最短路径，作为初始的CSI边界；

步骤9-3、采用四阶多项式拟合初始的CSI边界得到最终的CSI。

本发明与现有技术相比，其显著优点为：(1)本发明首次给出了一种基于 HD-OCT图像的自动脉络膜层分割方法，该方法能够鲁棒精确地分割脉络膜层的 上下边界；(2)本发明充分利用了脉络膜层下边界的灰度渐变特性，较传统的基 于灰度和梯度的方法具有更好的鲁棒性和定位精度。

下面结合附图对本发明作进一步详细描述。

附图说明

图1是本发明基于HD-OCT视网膜图像的脉络膜层自动分割方法的流程图。

图2是去除与ILM层相近的高反射率RNFL层的流程图。

图3是RPE层中轴线估计的流程图。

图4是CSI候选边界点生成的流程图。

图5是CSI最终边界点生成的流程图。

图6是HD-OCT视网膜图像内部组织示意图。

图7是BM分割及拉平结果图。

图8是BM下方的脉络膜窄带区域图。

图9是灰度渐变距离图像。

图10是灰度渐变距离垂直差图像。

图11是候选CSI边界点示意图。

图12是图搜索结果示意图。

图13是最终的CSI分割结果示意图。

具体实施方式

结合图1，本发明的基于HD-OCT视网膜图像的自动脉络膜层分割方法包括 以下步骤：

步骤1、采集HD-OCT视网膜图像，采用现有的OCT成像设备对视网膜图 像进行采集；

步骤2、采用改进的双边滤波对输入图像进行去噪处理。具体是将传统的双 边滤波的各项同性高斯邻域窗口改为各项异性的高斯邻域窗口，其中，传统的双 边滤波算法的公式为：

$h\left(x\right)=k^{-1}\left(x\right){\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f\left(\xi \right)c(\xi ,x)s(f\left(\xi \right),f\left(x\right))\mathrm{d\xi }$

式中f和h分别为输入和输出图像，函数c(ξ,x)用于测量邻域中心点x和邻域点 ξ之间的空间距离，函数s用于测量两点间的灰度相似性，函数c和函数s都是高 斯函数，$k\left(x\right)={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }c(\xi ,x)s(f\left(\xi \right),f\left(x\right))\mathrm{d\xi }$是归一化函数。

步骤3、根据玻璃体和视网膜的反射率差异定位ILM层。通过阈值得到玻璃 体区域，然后寻找玻璃体的最下方边界点从而得到ILM层，所述阈值为41。

步骤4、去除与ILM层相近的高反射率RNFL层。结合图2，具体为：

步骤4-1、生成模板图像M，其中ILM层像素置为1，其余像素置为0；

步骤4-2、用数学形态学对模板图像进行膨胀，所用公式为：

Ms=M⊕S

式中‘⊕’表示膨胀算子，S为半径30个像素的圆盘形结构元；

步骤4-3、将与Ms目标区域相连的高反射率RNFL层进行去除。

步骤5、利用高反射率和视网膜的结构特性估计RPE层。结合图3，具体为：

步骤5-1、将RNFL层以外的高反射率区域置为候选RPE层区域；

步骤5-2、根据RPE层的厚度约束去除厚度大于20个像素的虚假RPE区域；

步骤5-3、将每列的候选RPE区域的中心点作为RPE层的中轴点，然后采 用四阶多项式拟合得到RPE层的中轴线。

步骤6、在RPE层的下边界附近根据垂直梯度算法得到候选的BM边界点， 然后采用四阶多项式拟合得到最终的BM；RPE层的下边界附近为在RPE层中轴 线下方的20个像素窄带区域内，候选的BM边界点的选取方法为：

将每列中在此窄带区域内的最大垂直梯度点作为候选的BM边界点，然后采 用四阶多项式拟合得到最终的BM；

步骤7、基于分割得到的BM拉平输入图像，然后将CSI分割限定在BM下 方的200个像素厚的窄带区域；

步骤8、根据CSI下方像素的反射率渐变特性得到CSI的候选边界点。结合 图4，具体为：

步骤8-1、生成灰度渐变距离图像D

$D(x,y)=\left(\begin{array}{ccc}D(x+1,y)+1& \mathrm{if}& I_a(x,y)\ge I_a(x+1,y)\\ 0& \mathrm{if}& I_a(x,y)<I_a(x+1,y)\end{array}\right)$

式中(x,y)表示图像像素坐标，$I_a=\left(\begin{array}{c}1\\ -1\end{array}\right)*I$表示垂直方向的灰度差，I为平滑后 的拉平脉络膜区域图像，‘*’表示卷积算子；

步骤8-2、通过灰度渐变距离图像D生成灰度渐变距离图像的垂直差图像Dx

$D_x=\left(\begin{array}{c}-1\\ 1\end{array}\right)*D$

步骤8-3、将每列垂直差图像D_x中的最大值位置作为CSI边界的候选点。

步骤9、通过图搜索在候选CSI边界点中寻找初步的CSI边界，然后采用四 阶多项式拟合得到最终的CSI。结合图5，具体为：

步骤9-1、构造图G=(V,E)，其中候选CSI边界点作为图G中的顶点V， 每个顶点与后续的k列顶点相连，E为图G中边的集合，图中任意两个顶点a和 b之间的边权重通过下式计算：

$w(a,b)=w(\mathrm{\Delta x},\mathrm{\Delta y},w_M,T_p,\alpha )=(\Delta x^2+\Delta y^2)+\frac{w_{M}H(\mathrm{\Delta y}-T_p)|\mathrm{\Delta y}-T_p|}{1+\exp (-\alpha (\mathrm{\Delta y}-T_p))}$

式中Δx和Δy分别为两顶点的水平和垂直距离，H(·)为Heaviside函数，惩罚因 子w_M为2000，T_p为阈值5，当Δy≥T_p时，边将被赋予额外的很大的惩罚权重， 参数α控制惩罚权重的提升速率，α的取值为2；

步骤9-2、采用Dijkstra算法搜素图G的最短路径，作为初始的CSI边界；

步骤9-3、采用四阶多项式拟合初始的CSI边界得到最终的CSI。

下面结合实施例对本发明做进一步详细的说明：

本系统发明以HD-OCT视网膜图像作为输入，采用图像处理手段对输入图 像中的脉络膜层进行自动分割。

本实施例的流程如图1所示，通过OCT成像设备采集到的HD-OCT视网膜 图像大小为1024×1024，图6给出了一幅HD-OCT视网膜图像的感兴趣区域，图 中标注了视网膜的几个主要相关组织结构，如脉络膜的上下边界，BM和CSI。 为了便于后续的处理，首先对输入图像进行去噪处理，然后根据玻璃体和视网膜 的反射率差异定位ILM层，即通过阈值区分玻璃体和视网膜区域，对于固定的 成像设备，玻璃体和视网膜区域的反射率特性基本上是固定的，所以本发明将阈 值取为41。由于ILM即为RNFL层的上边界，所以与ILM相连的高反射率区域 基本上都是RNFL层，通过去除RNFL层的干扰，通过RPE层的高反射率特性 就可以粗略地得到RPE层。最后通过四阶多项式拟合RPE层的中轴线，最后在 RPE层中轴线下方的20个像素窄带区域内，将每列中在此窄带区域内的最大垂 直梯度点作为候选的BM边界点，然后采用四阶多项式拟合得到最终的BM。利 用BM可以得到拉平后的图像，图7给出了与图6对应的拉平图像，其中的白色 线为BM。

图8为与图7相对应的BM下方窄带区域图像，从中可以看出，脉络膜层的 CSI下方的灰度呈逐渐减小趋势，因此本发明采用公式 $D(x,y)=\left(\begin{array}{ccc}D(x+1,y)+1& \mathrm{if}& I_a(x,y)\ge I_a(x+1,y)\\ 0& \mathrm{if}& I_a(x,y)<I_a(x+1,y)\end{array}\right)$生成灰度渐变距离图像，如图 9所示，图中灰度值越高，表示灰度渐变距离越大。每列灰度渐变距离最大的点 基本上位于CSI上。

由于每列灰度渐变距离最大的点的上方点为黑色（即灰度渐变距离为0）， 所以通过公式$D_x=\left(\begin{array}{c}-1\\ 1\end{array}\right)*D$得到灰度渐变距离图像的垂直差图像，如图10所示。 将每列垂直差图像（图10）中的最大值位置作为CSI边界的候选点，如图11中 的黑色圆圈所示。

通过图搜索在候选CSI边界候选点中就可以得到初步的CSI边界，如图12 中的白色曲线所示。本发明中权重构造公式 $w(a,b)=w(\mathrm{\Delta x},\mathrm{\Delta y},w_M,T_p,\alpha )=(\Delta x^2+\Delta y^2)+\frac{w_{M}H(\mathrm{\Delta y}-T_p)|\mathrm{\Delta y}-T_p|}{1+\exp (-\alpha (\mathrm{\Delta y}-T_p))}$中的T_p取为 5，α为2，w_M为2000。为了进一步使CSI边界更光滑，采用四阶多项式拟合 图搜索得到最终的CSI，如图13中的黑色曲线所示，从图13可知：本发明得到 的CSI满足临床医生的判断依据，即CSI下方的灰度一致，且CSI边界光滑连 续。