技术领域
[0001] 本
发明涉及一种计算机
图像处理方法，具体涉及一种
眼底图像分割方法，尤其是对扫频光学相干断层成像得到的眼底图像中组织层次的分割方法。
背景技术
[0002] 脉络膜位于
视网膜与巩膜之间，具有复杂的血管结构，其主要功能是营养视网膜外层及玻璃体，并有遮光作用，使反射的物象清楚。同时对人的视觉系统起保护作用，对整个视觉神经有调节作用。
[0003] 很多
疾病会引起脉络膜形态的变化，如
青光眼、高度近视、年龄相关性黄斑变性、中心性
浆液性脉络膜视网膜病变、小柳原田综合症等。对脉络膜厚度和体积的精确定量分析对于发现早期病变、观测病程和研究病理都有着重要作用。
[0004] 中心
频率为1微米左右的扫频光学相干断层成像(SS-OCT)技术是最新的眼底
三维扫描技术，与传统的光学相干断层成像技术相比，该频率的光能更好地穿透玻璃体和视网膜组织，成像深度更深，因此能检测到完整的脉络膜上下界面。并且成像速度更快，可以在几秒钟内形成高
分辨率的三维扫描图像，成像范围从原来的6×6毫米提高到12×9毫米，图像中可同时包含黄斑和视盘区域。
[0005] 目前针对OCT图像中的脉络膜自动分割
算法大部分都是针对二维高清成像的，其
缺陷和不足在于：(1)这类方法没有充分利用三维数据中相邻
位置的相关信息，较容易受到图像噪声或伪影的影响，导致分割错误；(2)通过分割只能得到一条扫描线上的脉络膜厚度信息，而无法得到其他位置的厚度信息及脉络膜总体积信息；(3)这些算法通常只适用于以黄斑为中心的图像，而不适用于大视野扫描成像。
[0006] 而对于三维成像的处理大都是针对SD-OCT图像的，例如，中国发明
专利104050672A公开了一种基于SD-OCT视网膜图像的地图状萎缩投影图像生成方法，将三维的SD-OCT图像投影成二维的眼底图像，其后续处理仍然是针对二维图像的。而SD-OCT图像是用分光光度仪来分离不同
波长的光波的，其
信噪比与SS-OCT图像不同。因而，这些算法只适用于信噪比相对较高的图像，通常不适用于噪声强度较大的SS-OCT三维扫描图像。
[0007] 因此，有必要提供一种用于SS-OCT三维扫描图像的分割处理方法，以实现三维大视野成像中的脉络膜分割。
发明内容
[0008] 本发明的发明目的是提供一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，实现对SS-OCT图像中脉络膜的自动分割。
[0009] 为达到上述发明目的，本发明采用的技术方案是：一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，先获取扫频光学相干断层成像三维图像，再进行以下处理：
[0010] (1)图像增强处理，所述图像增强处理包括图像
对比度增强和散斑噪声去除；
[0011] (2)视盘区域检测：用改进的Canny算子检测视网膜上表面和内-外视网膜分界面，通过检测内-外视网膜分界面的缺失，得到视盘区域的二维投影；
[0012] 所述的改进的Canny算子的改进方法是，设图像梯度为G＝(Gx,Gz)，其中Gx,Gz分别为
水平方向和垂直方向的梯度，|G|1＝|Gx|+|min(0,Gz)|，|G|2＝|min(0,Gz)|，以|G|1和|G|2代替Canny算子中的梯度的模进行边界点判断；
[0013] (3)脉络膜上表面分割：在内-外视网膜分界面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面；
[0014] (4)脉络膜下表面分割：首先进行初步分割，在脉络膜下表面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面的初步分割结果；然后根据初步分割结果，计算梯度信息与区域信息相结合的代价函数，用三维图搜算法获得精确分割结果。
[0015] 上述技术方案中，步骤(1)中，所述
图像对比度增强的方法为，
[0016] 采用线性对比度拉伸方法，计算公式如下：
[0017]
[0018] 其中f(x,y,z)为扫频光学相干断层成像三维图像，fe(x,y,z)为对比度拉伸后的三维图像，x,y,z分别代表中B扫描图像的宽度方向、B扫描图像的
帧数和B扫描图像的高度方向，B扫描图像是指三维图像中的一帧二维图像，b为整幅图像的灰度最大值，a为图像下方背景区域的灰度均值。
[0019] 背景区域的灰度均值计算如下：
[0020]
[0021] 其中H取图像总高度的1/5～1/4。
[0022] 通过对比度拉伸，将灰度值小于等于a的噪声点置0，同时将图像灰度值拉伸并归一化到[0,1]区间，可以去除一部分噪声，同时提高前景和背景、脉络膜组织和脉络膜血管间的对比度。
[0023] 步骤(1)中，所述散斑噪声去除的方法是，
[0024] 用交叉双边
滤波器对每一张B扫描图像进行散斑噪声去除，设一张对比度增强后的B扫描图像为
[0025] I(x,z)＝fe(x,y0,z) (3)
[0026] 交叉双边滤波器计算公式如下：
[0027]
[0028] ω(x,z,k,l)＝d(x,z,k,l)r(x,z,k,l) (5)
[0029]
[0030]
[0031] 其中I(x,z)表示图像的灰度值，y0为该图像所在的帧号，Id(x,z)为去噪结果，ω为(x,z)点与(k,l)点间的权重，k和l分别为图像中任意一点的横、纵坐标，d为ω中的空间域权重分量，r为ω中的值域权重分量，IM(x,z)为对I(x,z)进行中值滤波的结果图像，σd为计算空间域权重的高斯标准差，取值为图像宽度的1/40～1/60，σr为计算值域权重的高斯标准差，取值为0.02～0.04，中值滤波窗口大小为3或5。
[0032] 交叉双边滤波结合了中值滤波与双边滤波的优点，解决了传统双边滤波处理脉冲型噪声时值域权重较小，导致平滑不足的问题，能够较好地去除三维扫频OCT图像中的散斑噪声。
[0033] 步骤(2)中，视盘区域检测包括以下两个步骤：辅助分界面检测和视盘区域提取。
[0034] (a)辅助分界面检测
[0035] 需要检测两个辅助分界面，即视网膜上边界和内-外视网膜分界面。因为这两个界面是视网膜中最明显的两个分界面，因此直接在每张二维B扫描图像Id(x,z)上进行检测。设图像梯度为G＝(Gx,Gz)，其中Gx,Gz分别为水平方向和垂直方向的梯度。视网膜上边界是自上向下暗到亮的边界，Gz为负值，但在视盘区域附近接近垂直边缘，Gz接近0，但Gx绝对值较大。因此在视网膜上边界上，式(8)取值较大。
[0036] |G|1＝|Gx|+|min(0,Gz)| (8)
[0037] 内-外视网膜分界面也是自上向下暗到亮的边界，但在视盘区域内部缺失，不存在垂直边缘。在这个分界面上，式(9)取值较大。
[0038] |G|2＝|min(0,Gz)| (9)
[0039] 传统的Canny算子通过计算梯度的模和
角度来判断某个
像素点是否为边界点，本方法采用的改进型Canny算子中，用式(8)或(9)代替梯度的模进行边界点判断，可获得两张边界点分布图E1和E2，分别用于提取视网膜上边界和内-外视网膜分界面，通过限制边界z方向梯度为负，可排除较多的干扰性边缘。在边界点分布图E1中，选择每一列的第一个边界点，形成视网膜上边界。在边界点分布图E2中，选择每一列中在已检测到的视网膜上边界下方的第一个边界点，形成内-外视网膜分界面。
[0040] (b)视盘区域提取
[0041] 在视盘区域内部，视网膜组织是断开的，因此对该区域内部的每一列，在视网膜上边界下方检测不到内-外视网膜分界面。因此遍历三维图像的每一列，可以得到视盘区域在x-y平面的二值投影图，即当这一列检测不到内-外视网膜分界面时，令该列在投影图上对应的点等于1，反之为0。在初步检测结果上进一步处理，只保留面积最大的连通区域，对内部孔洞进行填充，并用凸包进行边界平滑，得到最后的视盘二维投影区域，如图3所示。在这个区域内部，不分割脉络膜的上下边界。提取到视盘区域后，在每张B扫描图像上，在视盘内部区域对内-外视网膜分界面进行线性插值，以得到一个连续的分界面，作为对脉络膜分割的约束。
[0042] 步骤(3)中，用多分辨率三维图搜算法进行分割，获得代价最小的分界面，并且该分界面符合由Δx和Δz规定的平滑性约束，即在x方向相邻两列上分界面的高度变化不超过Δx，在y方向相邻两列上分界面的高度变化不超过Δy，对脉络膜上表面，取代价函数为三维索贝尔算子得到的z方向梯度的相反数，并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，设Δx＝Δy＝1或2，脉络膜上表面的分割在线性插值后的内-外视网膜分界面下方区域中进行，在多分辨率算法中，将三维图像在z方向上进行四次2阶下
采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的
基础上，在附近进一步精确分割，最终得到原图像上的分割结果。
[0043] 步骤(4)中，在得到脉络膜上表面后，基于上表面位置将视网膜图像平坦化，即将图像中各列上下移动，使得上表面变为一个平面，所述脉络膜下表面分割包括以下步骤：
[0044] (a)初步分割
[0045] 用和上表面分割类似的多分辨率三维图搜算法进行分割，因为下表面为暗到亮的分界面，取梯度代价函数
[0046]
[0047] 并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，设Δx＝Δy＝1或2，脉络膜下表面的分割在脉络膜上表面下方区域中进行，在多分辨率算法中，将三维图像在z方向上进行四次2阶
下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果；
[0048] (b)精确分割
[0049] 基于初步分割结果，计算基于
小波变换和渐进灰度下降距离的区域性代价函数，[0050] 首先，设图像中每一列为A(z)，计算其一维小波变换：
[0051]
[0052] 其中φJ,k(z)及ψj,k(z)分别为小波变换的尺度函数和小波函数，j和k是平移和尺度因子，aJ,k和dj,k为小波分解所得的近似系数和细节系数，J是分解层数，设所有dj,k＝0，可得到该列图像的平滑近似A0(z)，
[0053]
[0054] 本方法采用11-9双
正交小波进行4层小波变换，
[0055] 设z坐标值从下到上逐渐增大，WGID函数定义如下：
[0056]
[0057] 取初步分割结果下方，所有WGID的局部最大值计算区域信息代价函数，定义如下：
[0058]
[0059] 取(10)式中梯度信息代价函数与区域信息代价函数的加权和，形成用于精确分割的新的代价函数
[0060]
[0061] 其中w＝0.4～0.6，基于该代价函数，用三维图搜算法在拉平后的最高分辨率图像的脉络膜上表面下方区域中搜索代价最小的分界面，得到精确分割结果，视盘内部所有列的代价函数设为0，Δx＝Δy＝1或2。
[0062] 由于上述技术方案运用，本发明与
现有技术相比具有下列优点：
[0063] 1、本发明提供了一种针对中心频率为1微米左右的扫频光学相干断层成像大视野三维图像中脉络膜的自动分割方法，通过检测脉络膜的上下表面，即布鲁赫膜(BM)和脉络膜-巩膜分界面(CSI)，获得三维的脉络膜区域，根据分割可计算12×9毫米范围内任意位置的脉络膜厚度和脉络膜总体积。
[0064] 2、本发明通过对比度增强和去噪，采用改进的Canny算子检测视网膜上表面和内-外视网膜分界面，结合梯度和区域信息，先粗分割再精确分割实现脉络膜下边界分割，解决了SS-OCT图像中噪声强度较大导致的自动分割困难的问题，获得了显著的技术效果。
附图说明
[0065] 图1是
实施例中交叉双边滤波器去噪效果，其中，(a)原B扫描图像(部分)(b)对比度拉伸结果图像(c)去噪结果图像。
[0066] 图2是实施例中辅助界面检测结果图，其中，(a)用于检测视网膜上边界的边界点分布图E1，(b)用于检测内-外视网膜分界面的边界点分布图E2，(c)辅助界面及视盘检测结果：上方黑色曲线为视网膜上边界，下方白色曲线为内-外视网膜分界面，白色竖线表示视盘范围。
[0067] 图3是实施例中视盘区域检测结果，其中，(a)通过内-外视网膜分界面的空缺得到的初步检测结果，(b)去除小区域和孔洞填充后的结果(c)边界平滑后的结果。
[0068] 图4是实施例中基于小波变换和渐进灰度下降距离(WGID)的区域性代价函数计算示意图，其中，(a)拉平后的脉络膜上边界下方区域，白色虚线表示初始分割结果，白色实线表示最终分割结果，黑色实线表示手动分割结果(b)梯度代价函数，(c)WGID图像，(d)区域性代价函数。
[0069] 图5是实施例中脉络膜分割结果，其中，(a)含视盘和黄斑中心的B扫描图像中分割结果，黑色曲线表示上边界，白色曲线表示下边界，两条白色竖线表示视盘分割边界(b)周边区域的B扫描图像中分割结果(c)脉络膜厚度图。
具体实施方式
[0070] 下面结合附图及实施例对本发明作进一步描述：
[0071] 实施例一：
[0072] 一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，先获取扫频光学相干断层成像三维图像，再进行以下处理：图像增强、视盘区域检测、脉络膜上表面分割和脉络膜下表面分割。
[0073] 具体描述如下：
[0074] 1)图像增强：
[0075] 图像增强包括以下两个步骤：对比度增强和散斑噪声去除
[0076] (a)对比度增强
[0077] 采用线性对比度拉伸方法，计算公式如下：
[0078]
[0079] 其中f(x,y,z)为原三维图像，fe(x,y,z)为对比度拉伸后三维图像，x,y,z分别代表中B扫描图像的宽度方向、B扫描图像的帧数和B扫描图像的高度方向。b取整幅图像的灰度最大值，a取图像下方背景区域的均值，计算如下：
[0080]
[0081] 其中H取图像总高度的1/5～1/4。
[0082] 通过对比度拉伸，将灰度值小于等于a的噪声点置0，同时将图像灰度值拉伸并归一化到[0,1]区间，可以去除一部分噪声，同时提高前景和背景、脉络膜组织和脉络膜血管间的对比度。
[0083] (b)散斑噪声去除
[0084] 用交叉双边滤波器对每一张B扫描图像进行散斑噪声去除，设一张对比度增强后的B扫描图像为
[0085] I(x,z)＝fe(x,y0,z) (3)
[0086] 交叉双边滤波器计算公式如下：
[0087]
[0088] ω(x,z,k,l)＝d(x,z,k,l)r(x,z,k,l) (5)
[0089]
[0090]
[0091] 其中Id(x,z)为去噪结果，IM(x,z)为对I(x,z)进行中值滤波的结果图像。σd为计算空间域权重的高斯标准差，取值为图像宽度的1/40～1/60，σr为计算值域权重的高斯标准差，取值为0.02～0.04，中值滤波窗口大小为3或5。
[0092] 交叉双边滤波结合了中值滤波与双边滤波的优点，解决了传统双边滤波处理脉冲型噪声时值域权重较小，导致平滑不足的问题，能够较好的去除三维扫频OCT图像中的散斑噪声。
[0093] 图像对比度增强及去噪结果如图1所示。
[0094] 2)视盘区域检测：
[0095] 视盘区域检测包括以下两个步骤：辅助分界面检测和视盘区域提取。
[0096] (a)辅助分界面检测
[0097] 需要检测两个辅助分界面，即视网膜上边界和内-外视网膜分界面。因为这两个界面是视网膜中最明显的两个分界面，因此直接在每张二维B扫描图像Id(x,z)上进行检测。设图像梯度为G＝(Gx,Gz)，其中Gx,Gz分别为水平方向和垂直方向的梯度。视网膜上边界是自上向下暗到亮的边界，Gz为负值，但在视盘区域附近接近垂直边缘，Gz接近0，但Gx绝对值较大。因此在视网膜上边界上，式(8)取值较大。
[0098] |G|1＝|Gx|+|min(0,Gz)| (8)
[0099] 内-外视网膜分界面也是自上向下暗到亮的边界，但在视盘区域内部缺失，不存在垂直边缘。在这个分界面上，式(9)取值较大。
[0100] |G|2＝|min(0,Gz)| (9)
[0101] 传统的Canny算子通过计算梯度的模和角度来判断某个像素点是否为边界点，本方法采用的改进型Canny算子中，用式(8)或(9)代替梯度的模进行边界点判断，可获得两张边界点分布图E1和E2，分别用于提取视网膜上边界和内-外视网膜分界面，通过限制边界z方向梯度为负，可排除较多的干扰性边缘，如图2(a)(b)所示。在边界点分布图E1中，选择每一列的第一个边界点，形成视网膜上边界。在边界点分布图E2中，选择每一列中在已检测到的视网膜上边界下方的第一个边界点，形成内-外视网膜分界面。
[0102] (b)视盘区域提取
[0103] 在视盘区域内部，视网膜组织是断开的，因此对该区域内部的每一列，在视网膜上边界下方检测不到内-外视网膜分界面。因此遍历三维图像的每一列，可以得到视盘区域在x-y平面的二值投影图，即当这一列检测不到内-外视网膜分界面时，令该列在投影图上对应的点等于1，反之为0。在初步检测结果上进一步处理，只保留面积最大的连通区域，对内部孔洞进行填充，并用凸包进行边界平滑，得到最后的视盘二维投影区域，如图3所示。在这个区域内部，不分割脉络膜的上下边界。提取到视盘区域后，在每张B扫描图像上，在视盘内部区域对内-外视网膜分界面进行线性插值，以得到一个连续的分界面，如图2(b)所示，作为对脉络膜分割的约束。
[0104] 3)脉络膜上表面分割
[0105] 脉络膜上表面是亮到暗的分界面，本方法中用多分辨率三维图搜算法进行分割。三维图搜算法能够获得代价最小的分界面，并且该分界面符合由Δx和Δz规定的平滑性约束，即在x方向相邻两列上分界面的高度变化不超过Δx，在y方向相邻两列上分界面的高度变化不超过Δy。对脉络膜上表面，取代价函数 为三维索贝尔算子得到的z方向梯度的相反数，并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，以消除视盘内部灰度值变化对检测结果的影响。设Δx＝Δy＝1或2。脉络膜上表面的分割在线性插值后的内-外视网膜分界面下方区域中进行。在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果。
[0106] 4)脉络膜下表面分割：
[0107] 在得到脉络膜上表面后，基于上表面位置将视网膜图像平坦化，即将图像中各列上下移动，使得上表面变为一个平面，这样可使各B扫描图像在三维空间中对齐，并且减小检测下表面时的搜索范围。
[0108] 脉络膜下表面分割包括初步分割和精确分割两个步骤。
[0109] (a)初步分割
[0110] 用和上表面分割类似的多分辨率三维图搜算法进行分割，因为下表面为暗到亮的分界面，取梯度代价函数
[0111]
[0112] 并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，以消除视盘内部灰度值变化对检测结果的影响。设Δx＝Δy＝1或2。脉络膜下表面的分割在脉络膜上表面下方区域中进行。在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果。
[0113] (b)精确分割
[0114] 基于初步分割结果，计算基于小波变换和渐进灰度下降距离(WGID)的区域性代价函数。其依据是，在高信噪比OCT图像中，从脉络膜下边界开始，下方巩膜区域的灰度是单调下降的，即可通过灰度下降区间的起始点确定脉络膜下边界。但在低信噪比的三维扫频OCT图像中，需要用小波变换的平滑近似恢复这种单调下降的特性。
[0115] 首先，设图像中每一列为A(z)，计算其一维小波变换：
[0116]
[0117] 其中φJ,k(z)及ψj,k(z)分别为小波变换的尺度函数和小波函数,j和k是平移和尺度因子,aJ,k和dj,k为小波分解所得的近似系数和细节系数,J是分解层数。设所有dj,k＝0，可得到该列图像的平滑近似A0(z)，
[0118]
[0119] 本方法采用11-9双正交小波进行4层小波变换。
[0120] 设z坐标值从下到上逐渐增大，WGID函数定义如下：
[0121]
[0122] 取初步分割结果下方，所有WGID的局部最大值计算区域信息代价函数，定义如下：
[0123]
[0124] 取(10)式中梯度信息代价函数与区域信息代价函数的加权和，形成用于精确分割的新的代价函数
[0125]
[0126] 其中w＝0.4～0.6。基于该代价函数，用三维图搜算法在拉平后的最高分辨率图像的脉络膜上表面下方区域中搜索代价最小的分界面，得到精确分割结果，视盘内部所有列的代价函数仍设为0，设Δx＝Δy＝1或2。下表面分割结果及代价函数如图4所示。
[0127] 5)实验结果
[0128] 部分实验结果如图5所示。
[0129] 两个专家独立手动分割脉络膜上下边界，以结果位置的平均值为金标准。在16人的32只眼的OCT图像上进行测试，结果如表1，其中戴斯系数计算公式为：
[0130]
[0131] 其中CHauto及CHgt分别指本方法结果及金标准所得的脉络膜区域像素集合，|.|指集合中像素点个数。
[0132] 表1脉络膜分割结果误差分析
[0133]
[0134] 由表1可知，本方法的误差小于不同专家手动分割结果间的差异，且p-值很小，显示具有统计显著差异。本方法分割结果具有较高的准确性和客观性，能够替代手动分割。
[0135] 至此，一种适用于中心频率1微米左右的扫频
光学相干断层扫描图像的脉络膜自动分割方法已经实现并进行了验证。本发明融合了对比度增强、交叉双边滤波去噪、三维图搜技术、小波变换、渐进灰度下降距离计算等步骤，分割结果具有较高的准确性，能够替代手动分割，对于临床相关眼科疾病的诊断与
治疗能起到重要的辅助作用。
