CN103514605B - 基于hd-oct视网膜图像的脉络膜层自动分割方法 - Google Patents

Abstract

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

Description

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

技术领域

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

背景技术

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

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

（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){\∫}_{-\∞}^{\∞}{\∫}_{-\∞}^{\∞}f\left(\mathrm{\ξ}\right)c(\mathrm{\ξ},x)s(f\left(\mathrm{\ξ}\right),f\left(x\right))\mathrm{d\ξ}$

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

步骤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)\&GreaterEqual;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,\mathrm{\&alpha;})=(\mathrm{\&Delta;}{x}^2+\mathrm{\&Delta;}{y}^2)+\frac{w_{M}H(\mathrm{\&Delta;y}-T_p)|\mathrm{\&Delta;y}-T_p|}{1+\exp (-\mathrm{\&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){\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}f\left(\mathrm{\&xi;}\right)c(\mathrm{\&xi;},x)s(f\left(\mathrm{\&xi;}\right),f\left(x\right))\mathrm{d\&xi;}$

式中f和h分别为输入和输出图像，函数c(ξ,x)用于测量邻域中心点x和邻域点ξ之间的空间距离，函数s用于测量两点间的灰度相似性，函数c和函数s都是高斯函数， $k\left(x\right)={\&Integral;}_{-\&infin;}^{\&infin;}{\&Integral;}_{-\&infin;}^{\&infin;}c(\mathrm{\&xi;},x)s(f\left(\mathrm{\&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)\&GreaterEqual;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,\mathrm{\&alpha;})=(\mathrm{\&Delta;}{x}^2+\mathrm{\&Delta;}{y}^2)+\frac{w_{M}H(\mathrm{\&Delta;y}-T_p)|\mathrm{\&Delta;y}-T_p|}{1+\exp (-\mathrm{\&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)\&GreaterEqual;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,\mathrm{\&alpha;})=(\mathrm{\&Delta;}{x}^2+\mathrm{\&Delta;}{y}^2)+\frac{w_{M}H(\mathrm{\&Delta;y}-T_p)|\mathrm{\&Delta;y}-T_p|}{1+\exp (-\mathrm{\&alpha;}(\mathrm{\&Delta;y}-T_p))}$ 中的T_p取为5，α为2，w_M为2000。为了进一步使CSI边界更光滑，采用四阶多项式拟合图搜索得到最终的CSI，如图13中的黑色曲线所示，从图13可知：本发明得到的CSI满足临床医生的判断依据，即CSI下方的灰度一致，且CSI边界光滑连续。

Claims (7)

1.一种基于HD-OCT视网膜图像的脉络膜层自动分割方法，其特征在于，包括以下步骤：

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

步骤2、采用改进的双边滤波算法对输入图像进行去噪处理；

步骤3、根据玻璃体和视网膜的反射率差异定位ILM层；

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

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

步骤6、在RPE层的下边界附近根据垂直梯度算法得到候选的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& if& I_a(x,y)\&GreaterEqual;I_a(x+1,y)\\ 0& 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。

2.根据权利要求1所述的基于HD-OCT视网膜图像的脉络膜层自动分割方法，其特征在于，步骤2采用改进的双边滤波算法对输入图像进行去噪处理具体是将传统的双边滤波的各项同性高斯邻域窗口改为各项异性的高斯邻域窗口，其中，传统的双边滤波算法的公式为：

$h\left(x\right)=k^{-1}\left(x\right){\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}{\&Integral;}_{-\mathrm{\&infin;}}^{\mathrm{\&infin;}}f\left(\mathrm{\&xi;}\right)c(\mathrm{\&xi;},x)s\left(f\right(\mathrm{\&xi;}),f(x\left)\right)d\mathrm{\&xi;}$

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

3.根据权利要求1所述的基于HD-OCT视网膜图像的脉络膜层自动分割方法，其特征在于，步骤3中根据玻璃体和视网膜的反射率差异定位ILM层具体为：

通过阈值得到玻璃体区域，然后寻找玻璃体的最下方边界点从而得到ILM层，所述阈值为41。

4.根据权利要求1所述的基于HD-OCT视网膜图像的脉络膜层自动分割方法，其特征在于，步骤4去除与ILM层相近的高反射率RNFL层，具体为：

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

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

Ms＝M⊕S

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

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

5.根据权利要求1所述的基于HD-OCT视网膜图像的脉络膜层自动分割方法，其特征在于，步骤5利用高反射率和视网膜的结构特性估计RPE层的具体步骤为：

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

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

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

6.根据权利要求1所述的基于HD-OCT视网膜图像的脉络膜层自动分割方法，其特征在于，步骤6中RPE层的下边界附近为RPE层中轴线下方的20个像素窄带区域内；候选的BM边界点的选取方法为：

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

7.根据权利要求1所述的基于HD-OCT视网膜图像的脉络膜层自动分割方法，其特征在于，步骤9采用图搜索和多项式拟合得到最终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,\mathrm{\&alpha;})=({\mathrm{\&Delta;x}}^2+{\mathrm{\&Delta;y}}^2)+\frac{w_{M}H(\mathrm{\&Delta;}y-T_p)|\mathrm{\&Delta;}y-T_p|}{1+\exp (-\mathrm{\&alpha;}(\mathrm{\&Delta;}y-T_p\left)\right)}$

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

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

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