CN106558030B - 三维大视野扫频光学相干断层成像中脉络膜的分割方法 - Google Patents
三维大视野扫频光学相干断层成像中脉络膜的分割方法 Download PDFInfo
- Publication number
- CN106558030B CN106558030B CN201611005068.4A CN201611005068A CN106558030B CN 106558030 B CN106558030 B CN 106558030B CN 201611005068 A CN201611005068 A CN 201611005068A CN 106558030 B CN106558030 B CN 106558030B
- Authority
- CN
- China
- Prior art keywords
- image
- segmentation
- dimensional
- cost function
- interface
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Active
Links
- 230000011218 segmentation Effects 0.000 title claims abstract description 95
- 210000003161 choroid Anatomy 0.000 title claims abstract description 39
- 238000012014 optical coherence tomography Methods 0.000 title claims abstract description 35
- 238000000034 method Methods 0.000 title claims abstract description 33
- 230000002207 retinal effect Effects 0.000 claims abstract description 31
- 238000001514 detection method Methods 0.000 claims abstract description 25
- 210000001525 retina Anatomy 0.000 claims abstract description 15
- 238000010845 search algorithm Methods 0.000 claims abstract description 15
- 238000012545 processing Methods 0.000 claims abstract description 9
- 238000001914 filtration Methods 0.000 claims description 15
- 230000002146 bilateral effect Effects 0.000 claims description 14
- 238000004364 calculation method Methods 0.000 claims description 9
- 238000003384 imaging method Methods 0.000 claims description 7
- 238000000354 decomposition reaction Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 5
- 230000009466 transformation Effects 0.000 claims description 5
- 230000000750 progressive effect Effects 0.000 claims description 4
- 230000004256 retinal image Effects 0.000 claims description 4
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 claims description 3
- 238000006467 substitution reaction Methods 0.000 claims description 3
- 238000013519 translation Methods 0.000 claims description 3
- 230000002708 enhancing effect Effects 0.000 claims description 2
- 210000001519 tissue Anatomy 0.000 description 6
- 238000005516 engineering process Methods 0.000 description 5
- 238000000605 extraction Methods 0.000 description 4
- 238000010586 diagram Methods 0.000 description 3
- 201000010099 disease Diseases 0.000 description 3
- 208000037265 diseases, disorders, signs and symptoms Diseases 0.000 description 3
- 238000009499 grossing Methods 0.000 description 3
- 238000010408 sweeping Methods 0.000 description 3
- 210000004204 blood vessel Anatomy 0.000 description 2
- 210000001775 bruch membrane Anatomy 0.000 description 2
- 230000008859 change Effects 0.000 description 2
- 238000004891 communication Methods 0.000 description 2
- 230000003247 decreasing effect Effects 0.000 description 2
- 230000000694 effects Effects 0.000 description 2
- 238000002474 experimental method Methods 0.000 description 2
- PCHJSUWPFVWCPO-UHFFFAOYSA-N gold Chemical compound [Au] PCHJSUWPFVWCPO-UHFFFAOYSA-N 0.000 description 2
- 210000002189 macula lutea Anatomy 0.000 description 2
- 238000003672 processing method Methods 0.000 description 2
- 230000000007 visual effect Effects 0.000 description 2
- 206010003694 Atrophy Diseases 0.000 description 1
- 208000003569 Central serous chorioretinopathy Diseases 0.000 description 1
- 208000010412 Glaucoma Diseases 0.000 description 1
- 206010025421 Macule Diseases 0.000 description 1
- OAICVXFJPJFONN-UHFFFAOYSA-N Phosphorus Chemical compound [P] OAICVXFJPJFONN-UHFFFAOYSA-N 0.000 description 1
- 241001599101 Salix microphylla Species 0.000 description 1
- 206010064930 age-related macular degeneration Diseases 0.000 description 1
- 238000004458 analytical method Methods 0.000 description 1
- 230000037444 atrophy Effects 0.000 description 1
- 230000007547 defect Effects 0.000 description 1
- 238000003745 diagnosis Methods 0.000 description 1
- 230000004402 high myopia Effects 0.000 description 1
- 238000003709 image segmentation Methods 0.000 description 1
- 230000003902 lesion Effects 0.000 description 1
- 208000002780 macular degeneration Diseases 0.000 description 1
- 210000005036 nerve Anatomy 0.000 description 1
- 230000007170 pathology Effects 0.000 description 1
- 238000004445 quantitative analysis Methods 0.000 description 1
- 230000009467 reduction Effects 0.000 description 1
- 230000001105 regulatory effect Effects 0.000 description 1
- 210000003786 sclera Anatomy 0.000 description 1
- 208000011580 syndromic disease Diseases 0.000 description 1
- 238000012360 testing method Methods 0.000 description 1
- 230000002792 vascular Effects 0.000 description 1
- 210000004127 vitreous body Anatomy 0.000 description 1
Images
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T5/00—Image enhancement or restoration
- G06T5/70—Denoising; Smoothing
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/10—Image acquisition modality
- G06T2207/10072—Tomographic images
- G06T2207/10101—Optical tomography; Optical coherence tomography [OCT]
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/20—Special algorithmic details
- G06T2207/20172—Image enhancement details
- G06T2207/20182—Noise reduction or smoothing in the temporal domain; Spatio-temporal filtering
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06T—IMAGE DATA PROCESSING OR GENERATION, IN GENERAL
- G06T2207/00—Indexing scheme for image analysis or image enhancement
- G06T2207/30—Subject of image; Context of image processing
- G06T2207/30004—Biomedical image processing
- G06T2207/30041—Eye; Retina; Ophthalmic
Landscapes
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Engineering & Computer Science (AREA)
- Theoretical Computer Science (AREA)
- Investigating Or Analysing Materials By Optical Means (AREA)
- Eye Examination Apparatus (AREA)
Abstract
本发明公开了一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，进行以下处理：图像增强处理；视盘区域检测：用改进的Canny算子检测视网膜上表面和内‑外视网膜分界面，通过检测内‑外视网膜分界面的缺失，得到视盘区域的二维投影；脉络膜上表面分割：在内‑外视网膜分界面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面；脉络膜下表面分割：在脉络膜下表面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面的初步分割结果；根据初步分割结果，计算梯度信息与区域信息相结合的代价函数，用三维图搜算法获得精确分割结果。本发明对SS‑OCT图像获得三维的脉络膜区域，可计算12×9毫米范围内任意位置的脉络膜厚度和脉络膜总体积。
Description
技术领域
本发明涉及一种计算机图像处理方法，具体涉及一种眼底图像分割方法，尤其是对扫频光学相干断层成像得到的眼底图像中组织层次的分割方法。
背景技术
脉络膜位于视网膜与巩膜之间，具有复杂的血管结构，其主要功能是营养视网膜外层及玻璃体，并有遮光作用，使反射的物象清楚。同时对人的视觉系统起保护作用，对整个视觉神经有调节作用。
很多疾病会引起脉络膜形态的变化，如青光眼、高度近视、年龄相关性黄斑变性、中心性浆液性脉络膜视网膜病变、小柳原田综合症等。对脉络膜厚度和体积的精确定量分析对于发现早期病变、观测病程和研究病理都有着重要作用。
中心频率为1微米左右的扫频光学相干断层成像(SS-OCT)技术是最新的眼底三维扫描技术，与传统的光学相干断层成像技术相比，该频率的光能更好地穿透玻璃体和视网膜组织，成像深度更深，因此能检测到完整的脉络膜上下界面。并且成像速度更快，可以在几秒钟内形成高分辨率的三维扫描图像，成像范围从原来的6×6毫米提高到12×9毫米，图像中可同时包含黄斑和视盘区域。
目前针对OCT图像中的脉络膜自动分割算法大部分都是针对二维高清成像的，其缺陷和不足在于：(1)这类方法没有充分利用三维数据中相邻位置的相关信息，较容易受到图像噪声或伪影的影响，导致分割错误；(2)通过分割只能得到一条扫描线上的脉络膜厚度信息，而无法得到其他位置的厚度信息及脉络膜总体积信息；(3)这些算法通常只适用于以黄斑为中心的图像，而不适用于大视野扫描成像。
而对于三维成像的处理大都是针对SD-OCT图像的，例如，中国发明专利104050672A公开了一种基于SD-OCT视网膜图像的地图状萎缩投影图像生成方法，将三维的SD-OCT图像投影成二维的眼底图像，其后续处理仍然是针对二维图像的。而SD-OCT图像是用分光光度仪来分离不同波长的光波的，其信噪比与SS-OCT图像不同。因而，这些算法只适用于信噪比相对较高的图像，通常不适用于噪声强度较大的SS-OCT三维扫描图像。
因此，有必要提供一种用于SS-OCT三维扫描图像的分割处理方法，以实现三维大视野成像中的脉络膜分割。
发明内容
本发明的发明目的是提供一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，实现对SS-OCT图像中脉络膜的自动分割。
为达到上述发明目的，本发明采用的技术方案是：一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，先获取扫频光学相干断层成像三维图像，再进行以下处理：
(1)图像增强处理，所述图像增强处理包括图像对比度增强和散斑噪声去除；
(2)视盘区域检测：用改进的Canny算子检测视网膜上表面和内-外视网膜分界面，通过检测内-外视网膜分界面的缺失，得到视盘区域的二维投影；
所述的改进的Canny算子的改进方法是，设图像梯度为G＝(Gx,Gz)，其中Gx,Gz分别为水平方向和垂直方向的梯度，|G|1＝|Gx|+|min(0,Gz)|，|G|2＝|min(0,Gz)|，以|G|1和|G|2代替Canny算子中的梯度的模进行边界点判断；
(3)脉络膜上表面分割：在内-外视网膜分界面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面；
(4)脉络膜下表面分割：首先进行初步分割，在脉络膜下表面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面的初步分割结果；然后根据初步分割结果，计算梯度信息与区域信息相结合的代价函数，用三维图搜算法获得精确分割结果。
上述技术方案中，步骤(1)中，所述图像对比度增强的方法为，
采用线性对比度拉伸方法，计算公式如下：
其中f(x,y,z)为扫频光学相干断层成像三维图像，fe(x,y,z)为对比度拉伸后的三维图像，x,y,z分别代表中B扫描图像的宽度方向、B扫描图像的帧数和B扫描图像的高度方向，B扫描图像是指三维图像中的一帧二维图像，b为整幅图像的灰度最大值，a为图像下方背景区域的灰度均值。
背景区域的灰度均值计算如下：
其中H取图像总高度的1/5～1/4。
通过对比度拉伸，将灰度值小于等于a的噪声点置0，同时将图像灰度值拉伸并归一化到[0,1]区间，可以去除一部分噪声，同时提高前景和背景、脉络膜组织和脉络膜血管间的对比度。
步骤(1)中，所述散斑噪声去除的方法是，
用交叉双边滤波器对每一张B扫描图像进行散斑噪声去除，设一张对比度增强后的B扫描图像为
I(x,z)＝fe(x,y0,z) (3)
交叉双边滤波器计算公式如下：
ω(x,z,k,l)＝d(x,z,k,l)r(x,z,k,l) (5)
其中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。
交叉双边滤波结合了中值滤波与双边滤波的优点，解决了传统双边滤波处理脉冲型噪声时值域权重较小，导致平滑不足的问题，能够较好地去除三维扫频OCT图像中的散斑噪声。
步骤(2)中，视盘区域检测包括以下两个步骤：辅助分界面检测和视盘区域提取。
(a)辅助分界面检测
需要检测两个辅助分界面，即视网膜上边界和内-外视网膜分界面。因为这两个界面是视网膜中最明显的两个分界面，因此直接在每张二维B扫描图像Id(x,z)上进行检测。设图像梯度为G＝(Gx,Gz)，其中Gx,Gz分别为水平方向和垂直方向的梯度。视网膜上边界是自上向下暗到亮的边界，Gz为负值，但在视盘区域附近接近垂直边缘，Gz接近0，但Gx绝对值较大。因此在视网膜上边界上，式(8)取值较大。
|G|1＝|Gx|+|min(0,Gz)| (8)
内-外视网膜分界面也是自上向下暗到亮的边界，但在视盘区域内部缺失，不存在垂直边缘。在这个分界面上，式(9)取值较大。
|G|2＝|min(0,Gz)| (9)
传统的Canny算子通过计算梯度的模和角度来判断某个像素点是否为边界点，本方法采用的改进型Canny算子中，用式(8)或(9)代替梯度的模进行边界点判断，可获得两张边界点分布图E1和E2，分别用于提取视网膜上边界和内-外视网膜分界面，通过限制边界z方向梯度为负，可排除较多的干扰性边缘。在边界点分布图E1中，选择每一列的第一个边界点，形成视网膜上边界。在边界点分布图E2中，选择每一列中在已检测到的视网膜上边界下方的第一个边界点，形成内-外视网膜分界面。
(b)视盘区域提取
在视盘区域内部，视网膜组织是断开的，因此对该区域内部的每一列，在视网膜上边界下方检测不到内-外视网膜分界面。因此遍历三维图像的每一列，可以得到视盘区域在x-y平面的二值投影图，即当这一列检测不到内-外视网膜分界面时，令该列在投影图上对应的点等于1，反之为0。在初步检测结果上进一步处理，只保留面积最大的连通区域，对内部孔洞进行填充，并用凸包进行边界平滑，得到最后的视盘二维投影区域，如图3所示。在这个区域内部，不分割脉络膜的上下边界。提取到视盘区域后，在每张B扫描图像上，在视盘内部区域对内-外视网膜分界面进行线性插值，以得到一个连续的分界面，作为对脉络膜分割的约束。
步骤(3)中，用多分辨率三维图搜算法进行分割，获得代价最小的分界面，并且该分界面符合由Δx和Δz规定的平滑性约束，即在x方向相邻两列上分界面的高度变化不超过Δx，在y方向相邻两列上分界面的高度变化不超过Δy，对脉络膜上表面，取代价函数
为三维索贝尔算子得到的z方向梯度的相反数，并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，设Δx＝Δy＝1或2，脉络膜上表面的分割在线性插值后的内-外视网膜分界面下方区域中进行，在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果。
步骤(4)中，在得到脉络膜上表面后，基于上表面位置将视网膜图像平坦化，即将图像中各列上下移动，使得上表面变为一个平面，所述脉络膜下表面分割包括以下步骤：
(a)初步分割
用和上表面分割类似的多分辨率三维图搜算法进行分割，因为下表面为暗到亮的分界面，取梯度代价函数
并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，设Δx＝Δy＝1或2，脉络膜下表面的分割在脉络膜上表面下方区域中进行，在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果；
(b)精确分割
基于初步分割结果，计算基于小波变换和渐进灰度下降距离的区域性代价函数，
首先，设图像中每一列为A(z)，计算其一维小波变换：
其中φJ,k(z)及ψj,k(z)分别为小波变换的尺度函数和小波函数，j和k是平移和尺度因子，aJ,k和dj,k为小波分解所得的近似系数和细节系数，J是分解层数，设所有dj,k＝0，可得到该列图像的平滑近似A0(z)，
本方法采用11-9双正交小波进行4层小波变换，
设z坐标值从下到上逐渐增大，WGID函数定义如下：
取初步分割结果下方，所有WGID的局部最大值计算区域信息代价函数，定义如下：
取(10)式中梯度信息代价函数与区域信息代价函数的加权和，形成用于精确分割的新的代价函数
其中w＝0.4～0.6，基于该代价函数，用三维图搜算法在拉平后的最高分辨率图像的脉络膜上表面下方区域中搜索代价最小的分界面，得到精确分割结果，视盘内部所有列的代价函数设为0，Δx＝Δy＝1或2。
由于上述技术方案运用，本发明与现有技术相比具有下列优点：
1、本发明提供了一种针对中心频率为1微米左右的扫频光学相干断层成像大视野三维图像中脉络膜的自动分割方法，通过检测脉络膜的上下表面，即布鲁赫膜(BM)和脉络膜-巩膜分界面(CSI)，获得三维的脉络膜区域，根据分割可计算12×9毫米范围内任意位置的脉络膜厚度和脉络膜总体积。
2、本发明通过对比度增强和去噪，采用改进的Canny算子检测视网膜上表面和内-外视网膜分界面，结合梯度和区域信息，先粗分割再精确分割实现脉络膜下边界分割，解决了SS-OCT图像中噪声强度较大导致的自动分割困难的问题，获得了显著的技术效果。
附图说明
图1是实施例中交叉双边滤波器去噪效果，其中，(a)原B扫描图像(部分)(b)对比度拉伸结果图像(c)去噪结果图像。
图2是实施例中辅助界面检测结果图，其中，(a)用于检测视网膜上边界的边界点分布图E1，(b)用于检测内-外视网膜分界面的边界点分布图E2，(c)辅助界面及视盘检测结果：上方黑色曲线为视网膜上边界，下方白色曲线为内-外视网膜分界面，白色竖线表示视盘范围。
图3是实施例中视盘区域检测结果，其中，(a)通过内-外视网膜分界面的空缺得到的初步检测结果，(b)去除小区域和孔洞填充后的结果(c)边界平滑后的结果。
图4是实施例中基于小波变换和渐进灰度下降距离(WGID)的区域性代价函数计算示意图，其中，(a)拉平后的脉络膜上边界下方区域，白色虚线表示初始分割结果，白色实线表示最终分割结果，黑色实线表示手动分割结果(b)梯度代价函数，(c)WGID图像，(d)区域性代价函数。
图5是实施例中脉络膜分割结果，其中，(a)含视盘和黄斑中心的B扫描图像中分割结果，黑色曲线表示上边界，白色曲线表示下边界，两条白色竖线表示视盘分割边界(b)周边区域的B扫描图像中分割结果(c)脉络膜厚度图。
具体实施方式
下面结合附图及实施例对本发明作进一步描述：
实施例一：
一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，先获取扫频光学相干断层成像三维图像，再进行以下处理：图像增强、视盘区域检测、脉络膜上表面分割和脉络膜下表面分割。
具体描述如下：
1)图像增强：
图像增强包括以下两个步骤：对比度增强和散斑噪声去除
(a)对比度增强
采用线性对比度拉伸方法，计算公式如下：
其中f(x,y,z)为原三维图像，fe(x,y,z)为对比度拉伸后三维图像，x,y,z分别代表中B扫描图像的宽度方向、B扫描图像的帧数和B扫描图像的高度方向。b取整幅图像的灰度最大值，a取图像下方背景区域的均值，计算如下：
其中H取图像总高度的1/5～1/4。
通过对比度拉伸，将灰度值小于等于a的噪声点置0，同时将图像灰度值拉伸并归一化到[0,1]区间，可以去除一部分噪声，同时提高前景和背景、脉络膜组织和脉络膜血管间的对比度。
(b)散斑噪声去除
用交叉双边滤波器对每一张B扫描图像进行散斑噪声去除，设一张对比度增强后的B扫描图像为
I(x,z)＝fe(x,y0,z) (3)
交叉双边滤波器计算公式如下：
ω(x,z,k,l)＝d(x,z,k,l)r(x,z,k,l) (5)
其中Id(x,z)为去噪结果，IM(x,z)为对I(x,z)进行中值滤波的结果图像。σd为计算空间域权重的高斯标准差，取值为图像宽度的1/40～1/60，σr为计算值域权重的高斯标准差，取值为0.02～0.04，中值滤波窗口大小为3或5。
交叉双边滤波结合了中值滤波与双边滤波的优点，解决了传统双边滤波处理脉冲型噪声时值域权重较小，导致平滑不足的问题，能够较好的去除三维扫频OCT图像中的散斑噪声。
图像对比度增强及去噪结果如图1所示。
2)视盘区域检测：
视盘区域检测包括以下两个步骤：辅助分界面检测和视盘区域提取。
(a)辅助分界面检测
需要检测两个辅助分界面，即视网膜上边界和内-外视网膜分界面。因为这两个界面是视网膜中最明显的两个分界面，因此直接在每张二维B扫描图像Id(x,z)上进行检测。设图像梯度为G＝(Gx,Gz)，其中Gx,Gz分别为水平方向和垂直方向的梯度。视网膜上边界是自上向下暗到亮的边界，Gz为负值，但在视盘区域附近接近垂直边缘，Gz接近0，但Gx绝对值较大。因此在视网膜上边界上，式(8)取值较大。
|G|1＝|Gx|+|min(0,Gz)| (8)
内-外视网膜分界面也是自上向下暗到亮的边界，但在视盘区域内部缺失，不存在垂直边缘。在这个分界面上，式(9)取值较大。
|G|2＝|min(0,Gz)| (9)
传统的Canny算子通过计算梯度的模和角度来判断某个像素点是否为边界点，本方法采用的改进型Canny算子中，用式(8)或(9)代替梯度的模进行边界点判断，可获得两张边界点分布图E1和E2，分别用于提取视网膜上边界和内-外视网膜分界面，通过限制边界z方向梯度为负，可排除较多的干扰性边缘，如图2(a)(b)所示。在边界点分布图E1中，选择每一列的第一个边界点，形成视网膜上边界。在边界点分布图E2中，选择每一列中在已检测到的视网膜上边界下方的第一个边界点，形成内-外视网膜分界面。
(b)视盘区域提取
在视盘区域内部，视网膜组织是断开的，因此对该区域内部的每一列，在视网膜上边界下方检测不到内-外视网膜分界面。因此遍历三维图像的每一列，可以得到视盘区域在x-y平面的二值投影图，即当这一列检测不到内-外视网膜分界面时，令该列在投影图上对应的点等于1，反之为0。在初步检测结果上进一步处理，只保留面积最大的连通区域，对内部孔洞进行填充，并用凸包进行边界平滑，得到最后的视盘二维投影区域，如图3所示。在这个区域内部，不分割脉络膜的上下边界。提取到视盘区域后，在每张B扫描图像上，在视盘内部区域对内-外视网膜分界面进行线性插值，以得到一个连续的分界面，如图2(b)所示，作为对脉络膜分割的约束。
3)脉络膜上表面分割
脉络膜上表面是亮到暗的分界面，本方法中用多分辨率三维图搜算法进行分割。三维图搜算法能够获得代价最小的分界面，并且该分界面符合由Δx和Δz规定的平滑性约束，即在x方向相邻两列上分界面的高度变化不超过Δx，在y方向相邻两列上分界面的高度变化不超过Δy。对脉络膜上表面，取代价函数
为三维索贝尔算子得到的z方向梯度的相反数，并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，以消除视盘内部灰度值变化对检测结果的影响。设Δx＝Δy＝1或2。脉络膜上表面的分割在线性插值后的内-外视网膜分界面下方区域中进行。在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果。
4)脉络膜下表面分割：
在得到脉络膜上表面后，基于上表面位置将视网膜图像平坦化，即将图像中各列上下移动，使得上表面变为一个平面，这样可使各B扫描图像在三维空间中对齐，并且减小检测下表面时的搜索范围。
脉络膜下表面分割包括初步分割和精确分割两个步骤。
(a)初步分割
用和上表面分割类似的多分辨率三维图搜算法进行分割，因为下表面为暗到亮的分界面，取梯度代价函数
并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，以消除视盘内部灰度值变化对检测结果的影响。设Δx＝Δy＝1或2。脉络膜下表面的分割在脉络膜上表面下方区域中进行。在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果。
(b)精确分割
基于初步分割结果，计算基于小波变换和渐进灰度下降距离(WGID)的区域性代价函数。其依据是，在高信噪比OCT图像中，从脉络膜下边界开始，下方巩膜区域的灰度是单调下降的，即可通过灰度下降区间的起始点确定脉络膜下边界。但在低信噪比的三维扫频OCT图像中，需要用小波变换的平滑近似恢复这种单调下降的特性。
首先，设图像中每一列为A(z)，计算其一维小波变换：
其中φJ,k(z)及ψj,k(z)分别为小波变换的尺度函数和小波函数,j和k是平移和尺度因子,aJ,k和dj,k为小波分解所得的近似系数和细节系数,J是分解层数。设所有dj,k＝0，可得到该列图像的平滑近似A0(z)，
本方法采用11-9双正交小波进行4层小波变换。
设z坐标值从下到上逐渐增大，WGID函数定义如下：
取初步分割结果下方，所有WGID的局部最大值计算区域信息代价函数，定义如下：
取(10)式中梯度信息代价函数与区域信息代价函数的加权和，形成用于精确分割的新的代价函数
其中w＝0.4～0.6。基于该代价函数，用三维图搜算法在拉平后的最高分辨率图像的脉络膜上表面下方区域中搜索代价最小的分界面，得到精确分割结果，视盘内部所有列的代价函数仍设为0，设Δx＝Δy＝1或2。下表面分割结果及代价函数如图4所示。
5)实验结果
部分实验结果如图5所示。
两个专家独立手动分割脉络膜上下边界，以结果位置的平均值为金标准。在16人的32只眼的OCT图像上进行测试，结果如表1，其中戴斯系数计算公式为：
其中CHauto及CHgt分别指本方法结果及金标准所得的脉络膜区域像素集合，|.|指集合中像素点个数。
表1脉络膜分割结果误差分析
由表1可知，本方法的误差小于不同专家手动分割结果间的差异，且p-值很小，显示具有统计显著差异。本方法分割结果具有较高的准确性和客观性，能够替代手动分割。
至此，一种适用于中心频率1微米左右的扫频光学相干断层扫描图像的脉络膜自动分割方法已经实现并进行了验证。本发明融合了对比度增强、交叉双边滤波去噪、三维图搜技术、小波变换、渐进灰度下降距离计算等步骤，分割结果具有较高的准确性，能够替代手动分割，对于临床相关眼科疾病的诊断与治疗能起到重要的辅助作用。
Claims (1)
1.一种三维大视野扫频光学相干断层成像中脉络膜的分割方法，先获取扫频光学相干断层成像三维图像，再进行以下处理：
(1) 图像增强处理，所述图像增强处理包括图像对比度增强和散斑噪声去除；
(2) 视盘区域检测：用改进的Canny算子检测视网膜上表面和内-外视网膜分界面，通过检测内-外视网膜分界面的缺失，得到视盘区域的二维投影；
所述的改进的Canny算子的方法是，设图像梯度为
，其中
，
分别为水平方向和垂直方向的梯度，
，
，以
和
代替Canny算子中的梯度的模进行边界点判断；
(3) 脉络膜上表面分割：在内-外视网膜分界面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面；
(4) 脉络膜下表面分割：首先进行初步分割，在脉络膜下表面下方，去除视盘区域后，用三维多分辨率图搜算法，基于梯度代价函数，获得脉络膜上表面的初步分割结果；然后根据初步分割结果，计算梯度信息与区域信息相结合的代价函数，用三维图搜算法获得精确分割结果；
步骤(1)中，所述图像对比度增强的方法为，
线性对比度拉伸方法，计算公式如下：
其中
为扫频光学相干断层成像三维图像，
为对比度拉伸后的三维图像，
分别代表B扫描图像的宽度方向、B扫描图像的帧数和B扫描图像的高度方向，B扫描图像是指三维图像中的一帧二维图像，b为整幅图像的灰度最大值，a为图像下方背景区域的灰度均值；
步骤(1)中，所述散斑噪声去除的方法是，
用交叉双边滤波器对每一张B扫描图像进行散斑噪声去除，设一张对比度增强后的B扫描图像为
交叉双边滤波器计算公式如下：
其中I（x,z）表示图像的灰度值，y 0为该图像所在的帧号，
为去噪结果，ω为（x, z）点与（k,l）点间的权重，k和l分别为图像中任意一点的横、纵坐标，d为ω中的空间域权重分量，r为ω中的值域权重分量，
为对I（x,z）进行中值滤波的结果图像，
为计算空间域权重的高斯标准差，取值为图像宽度的1/40~1/60，
为计算值域权重的高斯标准差，取值为0.02~0.04，中值滤波窗口大小为3或5；
步骤(3)中，用多分辨率三维图搜算法进行分割，获得代价最小的分界面，并且该分界面符合由
和
规定的平滑性约束，即在x方向相邻两列上分界面的高度变化不超过
，在y方向相邻两列上分界面的高度变化不超过
，对脉络膜上表面，取代价函数
为三维索贝尔算子得到的z方向梯度的相反数，并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，设
= 
=1或2，脉络膜上表面的分割在线性插值后的内-外视网膜分界面下方区域中进行，在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果；
步骤(4)中，在得到脉络膜上表面后，基于上表面位置将视网膜图像平坦化，即将图像中各列上下移动，使得上表面变为一个平面，所述脉络膜下表面分割包括以下步骤：
（a）初步分割
用和上表面分割类似的多分辨率三维图搜算法进行分割，因为下表面为暗到亮的分界面，取梯度代价函数
并且根据视盘检测结果，将视盘内部所有列的代价函数设为0，设
=
=1或2，脉络膜下表面的分割在脉络膜上表面下方区域中进行，在多分辨率算法中，将三维图像在z方向上进行四次2阶下采样，得到4个分辨率的图像，先在低分辨率图像上进行分割，然后在高一级分辨率图像上，在已得结果的基础上，在附近进一步精确分割，最终得到原图像上的分割结果；
（b）精确分割
基于初步分割结果，计算基于小波变换和渐进灰度下降距离的区域性代价函数，
首先，设图像中每一列为
，计算其一维小波变换：
其中 
及 
分别为小波变换的尺度函数和小波函数，j和q分别是平移和尺度因子，
和 
分别为小波分解所得的近似系数和细节系数，J是分解层数，设所有
，可得到该列图像的平滑近似
，
本方法采用11-9双正交小波进行4层小波变换，
设z坐标值从下到上逐渐增大，基于小波变换的渐进灰度下降距离函数WGID定义如下：
取初步分割结果下方，所有基于小波变换的渐进灰度下降距离函数的局部最大值计算区域信息代价函数，定义如下：
取（10）式中梯度信息代价函数与区域信息代价函数的加权和，形成用于精确分割的新的代价函数
其中w = 0.4 ~ 0.6，基于该代价函数，用三维图搜算法在拉平后的最高分辨率图像的脉络膜上表面下方区域中搜索代价最小的分界面，得到精确分割结果，视盘内部所有列的代价函数设为0，
=
=1或2。
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611005068.4A CN106558030B (zh) | 2016-11-15 | 2016-11-15 | 三维大视野扫频光学相干断层成像中脉络膜的分割方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201611005068.4A CN106558030B (zh) | 2016-11-15 | 2016-11-15 | 三维大视野扫频光学相干断层成像中脉络膜的分割方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN106558030A CN106558030A (zh) | 2017-04-05 |
| CN106558030B true CN106558030B (zh) | 2020-04-03 |
Family
ID=58444001
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201611005068.4A Active CN106558030B (zh) | 2016-11-15 | 2016-11-15 | 三维大视野扫频光学相干断层成像中脉络膜的分割方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN106558030B (zh) |
Families Citing this family (12)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN107644424B (zh) * | 2017-10-09 | 2020-10-13 | 南京医科大学第一附属医院 | 通过合并edi-oct图像来分割sd-oct图像脉络膜的方法 |
| CN108961261B (zh) * | 2018-03-14 | 2022-02-15 | 中南大学 | 一种基于空间连续性约束的视盘区域oct图像层次分割方法 |
| CN108846827B (zh) * | 2018-04-16 | 2021-10-15 | 江南大学 | 一种基于多圆快速分割眼底视盘的方法 |
| CN108765388B (zh) * | 2018-05-17 | 2020-10-27 | 苏州大学 | 食道内窥oct图像层次结构的自动分割方法和系统 |
| CN109389568B (zh) * | 2018-10-25 | 2022-04-01 | 中国科学院上海光学精密机械研究所 | 自动测量皮肤光学相干层析图像中表皮厚度的方法 |
| CN109615634A (zh) * | 2018-12-13 | 2019-04-12 | 深圳大学 | 光学眼底图像分割方法、装置、计算机设备及存储介质 |
| CN110163838B (zh) * | 2019-04-01 | 2023-04-18 | 江西比格威医疗科技有限公司 | 一种视网膜oct图像分层算法 |
| CN110264472B (zh) * | 2019-05-29 | 2021-05-04 | 浙江科技学院 | 一种基于sd-oct的水果近皮下细胞无损成像方法 |
| CN110473297B (zh) * | 2019-08-20 | 2022-03-29 | 上海联影医疗科技股份有限公司 | 图像处理方法、装置、电子设备和存储介质 |
| CN111369510B (zh) * | 2020-02-28 | 2022-07-01 | 四川大学华西医院 | 一种自动估计脉络膜厚度的方法 |
| CN114332087B (zh) * | 2022-03-15 | 2022-07-12 | 杭州电子科技大学 | Octa图像的三维皮层表面分割方法及系统 |
| CN114677371B (zh) * | 2022-05-06 | 2022-09-02 | 南通市海视光电有限公司 | 一种基于机器视觉的化工界面检测方法 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103514605A (zh) * | 2013-10-11 | 2014-01-15 | 南京理工大学 | 基于hd-oct视网膜图像的脉络膜层自动分割方法 |
| CN105488350A (zh) * | 2015-12-02 | 2016-04-13 | 苏州大学 | 基于反应扩散模型的三维脉络膜新生血管生长预测方法 |
Family Cites Families (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US8914098B2 (en) * | 2009-03-08 | 2014-12-16 | Oprobe, Llc | Medical and veterinary imaging and diagnostic procedures utilizing optical probe systems |
| US9241626B2 (en) * | 2013-03-14 | 2016-01-26 | Carl Zeiss Meditec, Inc. | Systems and methods for improved acquisition of ophthalmic optical coherence tomography data |
-
2016
- 2016-11-15 CN CN201611005068.4A patent/CN106558030B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN103514605A (zh) * | 2013-10-11 | 2014-01-15 | 南京理工大学 | 基于hd-oct视网膜图像的脉络膜层自动分割方法 |
| CN105488350A (zh) * | 2015-12-02 | 2016-04-13 | 苏州大学 | 基于反应扩散模型的三维脉络膜新生血管生长预测方法 |
Non-Patent Citations (2)
| Title |
|---|
| Automated 3-D Retinal Layer Segmentation of Macular Optical Coherence Tomography Images with Serous Pigment Epithelial Detachments;Fei Shi 等;《IEEE Transactions on Medical Imaging》;20140930;第1-12页 * |
| 基于Hessian矩阵和区域生长的脉络膜血管自动检测;邢琳 等;《软件导刊》;20160630;第15卷(第6期);第132-136页 * |
Also Published As
| Publication number | Publication date |
|---|---|
| CN106558030A (zh) | 2017-04-05 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN106558030B (zh) | 三维大视野扫频光学相干断层成像中脉络膜的分割方法 | |
| Yin et al. | User-guided segmentation for volumetric retinal optical coherence tomography images | |
| US9589346B2 (en) | Segmentation and identification of layered structures in images | |
| EP3084717B1 (en) | Systems & methods for ocular anterior segment tracking, alignment, and dewarping using optical coherence tomography | |
| CN108961261B (zh) | 一种基于空间连续性约束的视盘区域oct图像层次分割方法 | |
| Hu et al. | Multimodal retinal vessel segmentation from spectral-domain optical coherence tomography and fundus photography | |
| CN107644424B (zh) | 通过合并edi-oct图像来分割sd-oct图像脉络膜的方法 | |
| US10573007B2 (en) | Image processing apparatus, image processing method, and image processing program | |
| Lang et al. | Segmentation of retinal OCT images using a random forest classifier | |
| CN114092405A (zh) | 一种针对黄斑水肿oct图像的视网膜层自动分割方法 | |
| JP7116763B2 (ja) | 光コヒーレンストモグラフィ（ｏｃｔ）画像の３次元シャドウを低減する信号処理方法 | |
| Shi et al. | Automated choroid segmentation in three-dimensional 1-μ m wide-view OCT images with gradient and regional costs | |
| CN110033496B (zh) | 时间序列三维视网膜sd-oct图像的运动伪差校正方法 | |
| Kayte et al. | Automated Screening of Diabetic Retinopathy Using Image Processing | |
| CN111861977A (zh) | 一种基于机器视觉的眼前节断层图像的特征提取方法 | |
| JP6469838B2 (ja) | 生物組織の３次元ボリュームを表す画像データを分析する方法 | |
| Salarian et al. | Acuurate segmentation of retina nerve fiber layer in OCT images | |
| Kafieh et al. | Optical Coherence Tomography noise reduction over learned dictionaries with introduction of complex wavelet for start dictionary | |
| Kabir | A rule based segmentation approaches to extract retinal blood vessels in fundus image | |
| Jo et al. | Comparison of evaluation parameters in the retinal layer between diabetic cystoid macular edema and postoperative cystoid macular edema after cataract surgery based on a hierarchical approach | |
| Santhi et al. | An efficient approach to locate optic disc center, blood vessels and macula in retinal images | |
| Ding et al. | 3d orthogonal sd-oct volumes registration for the enhancement of pores in lamina cribrosa | |
| CN113627230B (zh) | 一种基于机器视觉的视网膜oct图像自动分割方法 | |
| CN109118499B (zh) | 基于回溯最短路径算法的相干光断层图像的层分割方法 | |
| Sparavigna | GIMP and wavelets for medical image processing: Enhancing images of the fundus of the eye |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |