CN106600614B - 基于凹凸性的sd-oct视网膜图像cnv分割方法 - Google Patents

Abstract

本发明公开了一种基于凹凸性的频域光学相干断层(SD‑OCT)视网膜图像脉络膜新生血管(CNV)分割方法，属于图像处理技术领域。该方法首先对输入的SD‑OCT图像估计视网膜和脉络膜区域，定位内界膜(ILM)和脉络膜‑巩膜分界面(CSJ)，然后根据视网膜图像的反射率渐变特性估计视网膜色素上皮(RPE)层，利用RPE层的凹凸性估计布鲁赫膜(BM)层，最后根据RPE和BM层的厚度差估计得到初步的CNV区域，通过修正CNV的上边界得到最终的CNV分割结果。实验结果表明，本发明中所给出算法能够鲁棒精确地分割CNV，对方便后续的CNV定量分析和提高医生的工作效率具有重要意义。

Description

基于凹凸性的SD-OCT视网膜图像CNV分割方法

技术领域

本发明涉及一种CNV自动分割方法，特别是一种基于频域光学相干断层(SD-OCT)视网膜图像的CNV自动分割方法。

背景技术

脉络膜新生血管(CNV)是老年黄斑变性(AMD)晚期的一种主要表现形式，会导致视网膜下出血、下积液渗出等影响视力的病变出现。传统的CNV病变测量主要基于荧光造影和靛青绿血管造影等二维成像技术，光学相干(OCT)成像技术可以有效地获得眼底的三维图像，基于OCT血管造影图像可以测量CNV的体积等三维参数，从而更加有效地分析CNV。

目前，只有基于OCT血管造影图像的自动CNV分割方法，还没有基于频域OCT(SD-OCT)的CNV自动分割方法，因为CNV在SD-OCT图像中不能像OCT血管造影图像一样，通过颜色差异较容易地分辨。

发明内容

本发明的目的在于提供一种基于凹凸性的SD-OCT视网膜图像CNV分割方法。

实现本发明的目的的技术解决方案为：一种基于凹凸性的SD-OCT视网膜图像CNV分割方法，包括以下步骤：

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

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

式中f和h分别为输入和输出图像，函数c(ξ,x)用于测量邻域中心点x和邻域点ξ之间的空间距离，函数s用于测量两点间的灰度相似性，函数c和函数s均为高斯函数，是归一化函数。

步骤3、根据反射率特性估计视网膜和脉络膜区域；具体为：

采用模板大小为95×95均值滤波平滑步骤2去噪后的图像，得到图像I，通过求图像I每行的平均值得到垂直投影图像PI，然后采用阈值t对平滑后的图像I进行二值化处理，取二值化结果中面积最大的连通区域为视网膜和脉络膜区域：

t＝PI_{masx_i}+4

其中PI_{masx_i}表示投影图像PI中出现次数最多的灰度值。

步骤4、定位ILM和CSJ边界；具体为：

步骤4-1、采用模板大小为61×1均值滤波平滑步骤2去噪后的图像，得到图像I′；

步骤4-2、生成灰度渐变距离图像D

式中(x,y)表示图像像素坐标，表示垂直方向的灰度差，‘*’表示卷积算子；

步骤4-3、生成灰度渐变距离图像的垂直差分图像同时将垂直差分图像中处于视网膜和脉络膜区域以外的区域置为0；

步骤4-4、取垂直差分图像每列中值最大的前两个点中横坐标较小者为初步的内界膜ILM边界点，然后采用三阶多项式拟合去除错误的内界膜ILM边界点，最后采用直线连接剩下的内界膜ILM边界点得到最终的内界膜ILM边界；

步骤4-5、取步骤3所得视网膜和脉络膜区域的下边界近似为CSJ边界。

步骤5、根据视网膜图像的反射率渐变特性估计RPE层；具体为：

步骤5-1、将垂直差分图像中靠近内界膜ILM边界距离19个像素以内的区域置为0；

步骤5-2、归一化步骤5-1得到的垂直差分图像，然后取归一化垂直差分图像每列中的最大值，如果最大值大于0.2，则将该点置为RPE边界点；

步骤5-3、采用9×9中值滤波平滑步骤5-2得到的视网膜色素上皮RPE边界点得到最终的视网膜色素上皮RPE层。

步骤6、利用RPE层的凹凸性估计BM层；具体为：

步骤6-1、首先在x方向通过如下两步对RPE层进行凸包拟合：(1)将RPE上方区域都填充成目标区域，(2)采用Matlab函数regionprops拟合得到目标区域的凸包；然后将拟合得到的凸包区域的下边界作为初步的BM布鲁赫膜边界；

步骤6-2、在y方向对初步的布鲁赫膜BM边界再采用函数regionprops进行一次凸包拟合，然后将拟合得到的凸包区域的下边界作为最终的BM边界。

步骤7、根据RPE和BM的厚度差估计初步的CNV；具体为：

步骤7-1、生成视网膜色素上皮RPE和布鲁赫膜BM的厚度差图像DH，然后采用高阈值thr_high和低阈值thr_low求厚度差图像DH的二值图像BW_low和BW_high，然后将BW_high中的目标像素作为种子点在BW_low图像中进行种子生长得到初步的CNV投影图像，最后对CNV投影图像中每个连通区域根据面积和平均厚度去除虚假的CNV，即连续区域面积小于15个像素，或者连通区域面积小于200个像素且归一化后的平均厚度小于0.2，就被作为虚假CNV去除；

步骤7-2、对去除虚假CNV的投影图像的目标或CNV区域进行空洞填充得到初步的CNV。

步骤8、修正CNV的上边界得到最终的CNV分割结果。具体为：

步骤8-1、对每个CNV区域的上边界进行修正，在步骤2的去噪图像中，以当前CNV的上边界点为起始点，在该点的垂直向上100个像素区域内寻找与起始点灰度值差异在15以内的连通区域的上边界作为更新后的CNV上边界点；

步骤8-2、对更新后的CNV上表面采用窗口大小为5×5，标准差为1.5的高斯滤波平滑后得到最终的CNV分割结果。

本发明与现有技术相比，其显著优点为：本发明考虑了视网膜黄斑区RPE层的形态特征和反射率特性，克服了CNV病变边界模糊甚至缺失的困难，同时克服了水肿等视网膜其它病变对CNV分割的影响，能够全自动、快速、鲁邦、精确地分割CNV。

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

附图说明

图1是本发明基于凹凸性的SD-OCT视网膜图像CNV分割方法的流程图。

图2是定位ILM和CSJ边界的流程图。

图3是根据视网膜图像的反射率渐变特性估计RPE层的流程图。

图4是利用RPE层的凹凸性估计BM层的流程图。

图5是根据RPE和BM的厚度差估计初步的CNV的流程图。

图6是修正CNV的上边界得到最终的CNV分割结果的流程图。

图7是一帧带组织标注信息的原始SD-OCT视网膜图像。

图8是双边滤波平滑结果图像。

图9是估计得到的视网膜和脉络膜区域图像。

图10是视网膜和脉络膜区域内的灰度渐变距离图像。

图11是灰度渐变距离图像的垂直差分图像。

图12是ILM和CSJ边界定位结果。

图13是RPE层检测结果。

图14是BM边界检测结果。

图15是CNV初步分割结果。

图16是最终的CNV分割结果。

具体实施方式

本发明的一种基于凹凸性的频域光学相干断层(SD-OCT)视网膜图像脉络膜新生血管(CNV)分割方法，首先对输入的SD-OCT图像估计视网膜和脉络膜区域，定位内界膜(ILM)和脉络膜-巩膜分界面(CSJ)，然后根据视网膜图像的反射率渐变特性估计视网膜色素上皮(RPE)层，利用RPE层的凹凸性估计布鲁赫膜(BM)层，最后根据RPE和BM层的厚度差估计得到初步的CNV区域，通过修正CNV的上边界得到最终的CNV分割结果。

结合图1，本发明的基于SD-OCT视网膜图像的地图状萎缩投影图像生成方法包括以下步骤：

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

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

式中f和h分别为输入和输出图像，函数c(ξ,x)用于测量邻域中心点x和邻域点ξ之间的空间距离，函数s用于测量两点间的灰度相似性，函数c和函数s都是高斯函数，是归一化函数。

步骤3、根据反射率特性估计视网膜和脉络膜区域。采用模板大小为95×95均值滤波平滑步骤2去噪后的图像，得到图像I，通过求图像I每行的平均值得到垂直投影图像PI，然后采用阈值t对平滑后的图像I进行二值化处理，取二值化结果中面积最大的连通区域为视网膜和脉络膜区域：

t＝PI_{masx_i}+4

其中PI_{masx_i}表示投影图像PI中出现次数最多的灰度值。

步骤4、定位ILM和CSJ边界。结合图2，具体为：

步骤4-1、采用模板大小为61×1均值滤波平滑步骤2去噪后的图像，得到图像I′；

步骤4-2、生成灰度渐变距离图像D

式中(x,y)表示图像像素坐标，表示垂直方向的灰度差，‘*’表示卷积算子；

步骤4-3、生成灰度渐变距离图像的垂直差分图像同时将垂直差分图像中处于视网膜和脉络膜区域以外的区域置为0；

步骤4-4、取垂直差分图像中每列中值最大的前两个中横坐标较小者为初步的ILM边界点，然后采用三阶多项式拟合去除错误的ILM边界点，最后采用直线连接剩下的ILM边界点得到最终的ILM边界；

步骤4-5、取步骤3所得视网膜和脉络膜区域的下边界近似CSJ边界。

步骤5、根据视网膜图像的反射率渐变特性估计RPE层，结合图3，具体为：

步骤5-1、将垂直差分图像中靠近ILM边界(距离19个像素以内)的区域置为0，以消除ILM边界对RPE层检测的影响；

步骤5-2、归一化步骤5-1得到的垂直差分图像，然后取归一化垂直差分图像每列中的最大值，如果最大值大于0.2，则将该点置为RPE边界点；

步骤5-3、采用9×9中值滤波平滑步骤5-2得到的RPE边界点得到最终的RPE层；

步骤6、利用RPE层的凹凸性估计BM层，结合图4，具体为：

步骤6-1、首先在x方向通过如下两步对RPE层进行凸包拟合：(1)将RPE上方区域都填充成目标区域，(2)采用Matlab函数regionprops拟合得到目标区域的凸包；然后将拟合得到的凸包区域的下边界作为初步的BM布鲁赫膜边界；

步骤6-2、在y方向对初步的BM边界再采用函数regionprops进行一次凸包拟合，然后将拟合得到的凸包区域的下边界作为最终的BM边界；

步骤7、根据RPE和BM的厚度差估计初步的CNV，结合图5，具体为：

步骤7-1、生成RPE和BM的厚度差图像DH，然后采用高阈值thr_high和低阈值thr_low求厚度差图像DH的二值图像BW_low和BW_high，然后将BW_high中的目标像素作为种子点在BW_low图像中进行种子生长得到初步的CNV投影图像，最后对CNV投影图像中每个连通区域根据面积和平均厚度去除虚假的CNV，即连续区域面积小于15个像素，或者连通区域面积小于200个像素且归一化后的平均厚度小于0.2，就被作为虚假CNV去除；

步骤7-2、对去除虚假CNV的投影图像的目标(或CNV)区域进行空洞填充得到初步的CNV；

步骤8、修正CNV的上边界得到最终的CNV分割结果，结合图6，具体为：

步骤8-1、对每个CNV区域的上边界进行修正，因为根据反射率渐变特性得到的RPE边界都要低于真实的RPE上边界，所以在CNV上边界往上的一定距离邻域内寻找真实的RPE边界。在步骤2的去噪图像中，以当前CNV的上边界点为起始点，在该点的垂直向上100个像素区域内寻找与起始点灰度值差异在15以内的连通区域的上边界作为更新后的CNV上边界点；

步骤8-2、对更新后的CNV上表面采用窗口大小为5×5，标准差为1.5的高斯滤波平滑后得到最终的CNV分割结果。

本发明考虑了视网膜黄斑区RPE层的形态特征和反射率特性，克服了CNV病变边界模糊甚至缺失的困难，同时克服了水肿等视网膜其它病变对CNV分割的影响，能够全自动、快速、鲁邦、精确地分割CNV。

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

实施例

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

本实施例的流程如图1所示，通过OCT成像设备采集到的三维SD-OCT视网膜图像大小为1024×512×128，对应视网膜2mm×6mm×6mm的区域，图7给出了一帧原始的SD-OCT视网膜图像，图中标注了视网膜的几个主要相关组织结构(ILM：内界膜，RPE：视网膜色素上皮层，BM：布鲁赫膜，CNV：脉络膜新生血管，CSI：脉络膜巩膜边界)。图8为双边滤波平滑结果，图9为基于反射率特性估计得到的视网膜和脉络膜区域，由于视网膜和脉络膜区域的反射率明显高于玻璃体和巩膜等其它区域，所以通过一个全局阈值就可以容易地得到。图10为视网膜和脉络膜区域内的灰度渐变距离图像。对图10进行垂直差分得到图11所示的灰度渐变距离图像的垂直差分图像。对垂直差分图像中上方的强边缘点进行三阶多项式拟合得到近似的ILM边界(图12中的黑线)，取图9中感兴趣(白色)区域的下边界作为CSJ边界(图12中的白线)。然后在ILM和CSJ层之间寻找垂直差分图像梯度较大的点作为RPE，图13中的白色曲线即为得到的RPE层。对RPE表面进行x和y两个方向的凸包拟合得到BM边界(图14中的黑线)。根据RPE和BM的厚度差估计初步的CNV，如图15中白线表示CNV的初步边界，但初步CNV区域的上边界离真实CNV边界有一定误差，因此采用灰度差异更新CNV上边界，最后通过高斯平滑得到最终的CNV分割结果，如图16所示。

从图16可知，本发明得到的CNV分割结果精度较高，能够克服视网膜层形态变化剧烈和CNV下边界缺失的问题，这为临床医生的疾病诊断提供了方便。

Claims (8)

1.一种基于凹凸性的SD-OCT视网膜图像CNV分割方法，其特征在于，包括以下步骤：

步骤1、采集频域光学相干断层SD-OCT视网膜图像；

步骤2、采用双边滤波算法对步骤1采集的SD-OCT视网膜图像进行去噪处理；

步骤3、根据反射率特性估计视网膜和脉络膜区域；

步骤4、定位内界膜ILM和脉络膜-巩膜分界面CSJ的边界；

步骤5、根据视网膜图像的反射率渐变特性估计步骤1采集的SD-OCT视网膜色素上皮RPE层；

步骤6、利用视网膜色素上皮RPE层的凹凸性估计布鲁赫膜BM层；

步骤7、根据视网膜色素上皮RPE和布鲁赫膜BM层的厚度差估计初步的视网膜图像脉络膜新生血管CNV；

步骤8、修正脉络膜新生血管CNV的上边界得到最终的脉络膜新生血管CNV分割结果。

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

式中f和h分别为输入和输出图像，函数c(ξ,x)用于测量邻域中心点x和邻域点ξ之间的空间距离，函数s用于测量两点间的灰度相似性，函数c和函数s均为高斯函数，是归一化函数。

3.根据权利要求1所述的基于凹凸性的SD-OCT视网膜图像CNV分割方法，其特征在于，步骤3根据反射率特性估计视网膜和脉络膜区域，具体为：

采用模板大小为95×95均值滤波平滑步骤2去噪后的图像，得到图像I，通过求图像I每行的平均值得到垂直投影图像PI，然后采用阈值t对平滑后的图像I进行二值化处理，取二值化结果中面积最大的连通区域为视网膜和脉络膜区域：

t＝PI_{masx_i}+4

其中PI_{masx_i}表示投影图像PI中出现次数最多的灰度值。

4.根据权利要求1所述的基于凹凸性的SD-OCT视网膜图像CNV分割方法，其特征在于，步骤4定位内界膜ILM和脉络膜-巩膜分界面CSJ的边界具体为：

步骤4-1、采用模板大小为61×1均值滤波平滑步骤2去噪后的图像，得到图像I′；

步骤4-2、生成灰度渐变距离图像D

式中(x,y)表示图像像素坐标，表示垂直方向的灰度差，‘*’表示卷积算子；

步骤4-3、生成灰度渐变距离图像的垂直差分图像同时将垂直差分图像中处于视网膜和脉络膜区域以外的区域置为0；

步骤4-4、取垂直差分图像每列中值最大的前两个点中横坐标较小者为初步的内界膜ILM边界点，然后采用三阶多项式拟合去除错误的内界膜ILM边界点，最后采用直线连接剩下的内界膜ILM边界点得到最终的内界膜ILM边界；

步骤4-5、取步骤3所得视网膜和脉络膜区域的下边界近似为CSJ边界。

5.根据权利要求4所述的基于凹凸性的SD-OCT视网膜图像CNV分割方法，其特征在于，步骤5根据视网膜图像的反射率渐变特性估计步骤1采集的SD-OCT视网膜色素上皮RPE层，具体为：

步骤5-1、将垂直差分图像中靠近内界膜ILM边界距离19个像素以内的区域置为0；

步骤5-2、归一化步骤5-1得到的垂直差分图像，然后将归一化垂直差分图像每列中的最大值大于0.2的像素点置为RPE边界点；

步骤5-3、采用9×9中值滤波平滑步骤5-2得到的视网膜色素上皮RPE边界点得到最终的视网膜色素上皮RPE层。

6.根据权利要求1所述的基于凹凸性的SD-OCT视网膜图像CNV分割方法，其特征在于，步骤6利用视网膜色素上皮RPE层的凹凸性估计布鲁赫膜BM层，具体为：

步骤6-1、首先在x方向通过如下两步对RPE层进行凸包拟合：(1)将RPE上方区域都填充成目标区域，(2)采用Matlab函数regionprops拟合得到目标区域的凸包；然后将拟合得到的凸包区域的下边界作为初步的BM布鲁赫膜边界；

步骤6-2、在y方向对初步的布鲁赫膜BM边界再采用函数regionprops进行一次凸包拟合，然后将拟合得到的凸包区域的下边界作为最终的BM边界。

7.根据权利要求1所述的基于凹凸性的SD-OCT视网膜图像CNV分割方法，其特征在于，步骤7根据视网膜色素上皮RPE和布鲁赫膜BM的厚度差估计初步的CNV，具体为：

步骤7-1、生成视网膜色素上皮RPE和布鲁赫膜BM的厚度差图像DH，然后采用高阈值thr_high和低阈值thr_low求厚度差图像DH的二值图像BW_low和BW_high，然后将BW_high中的目标像素作为种子点在BW_low图像中进行种子生长得到初步的CNV投影图像，最后对CNV投影图像中每个连通区域根据面积和平均厚度去除虚假的CNV，即连续区域面积小于15个像素，或者连通区域面积小于200个像素且归一化后的平均厚度小于0.2，就被作为虚假CNV去除；

步骤7-2、对去除虚假CNV的投影图像的CNV区域进行空洞填充得到初步的CNV。

8.根据权利要求1所述的基于凹凸性的SD-OCT视网膜图像CNV分割方法，其特征在于，步骤8修正CNV的上边界得到最终的CNV分割结果，具体为：

步骤8-1、对每个CNV区域的上边界进行修正，在步骤2的去噪图像中，以当前CNV的上边界点为起始点，在该点的垂直向上100个像素区域内寻找与起始点灰度值差异在15以内的连通区域的上边界作为更新后的CNV上边界点；

步骤8-2、对更新后的CNV上边界点采用窗口大小为5×5，标准差为1.5的高斯滤波平滑后得到最终的CNV分割结果。