CN104050664A - 基于oct影像的眼前房角开放程度多特征分类方法 - Google Patents

Abstract

本发明涉及一种在光学相干断层扫描(OpticalCoherenceTomography,OCT)影像中对眼前房角开放程度采用多特征进行分类的方法。本发明首先对眼前房角的OCT图像进行分割处理，从分割图像中获取角膜后边缘与虹膜前边缘，并以两条边缘交点为基准，取前房角中包括完整房角隐窝的部分为感兴趣区域；然后在感兴趣区域内分别求取梯度、共生梯度方向直方图以及Edgelet等三个图像特征对房角隐窝区域从局部到整体进行全面描述；最后采用支持向量机根据获取的三个特征对图像中的房角开放程度做出分类。该发明能实现对前房角闭合状态的准确与快速分类。

Description

基于OCT影像的眼前房角开放程度多特征分类方法

技术领域

本发明属于医学图像处理领域，涉及一种新的基于OCT影像的眼前房角开放程度多特征分类方法。

背景技术

青光眼是世界范围内致盲的主要原因，大致可以分为开角型和闭角型两大类。房角闭合是原发性闭角型青光眼(primary angle closure glaucoma，PACG)的发病机制。据估计，到2020年全世界将有2100万闭角型青光眼患者，其中有1000万为中国患者。其病理为，周边虹膜与小梁网相接触或与小梁网产生永久性粘连，房水外流受阻，引起眼压升高。因此，前房角开放(或闭合)程度的评价对原发性闭角型青光眼的排查至关重要。频域OCT(FD-OCT)能够在短时间内获得清晰的眼前房角图像，且图像受眼动误差影响较小，因此其已逐渐成为眼前房角的主要成像手段。但由于OCT成像质量受噪声影响较大，在帧平均失败的情况下易造成前房角内解剖结构细节的不明显，直接导致已有的基于解剖结构特征参数测量或者基于图像视觉特征分类的方法无法正确地评估房角闭合状态。而已有的基于图像视觉特征分类的方法，对于成像清晰的图像也仅能从大类上简单区分眼前房角是开放或闭合两个状态，无法实现眼前房角的细化分类。根据图像分类的经验，只有在足够高的特征向量空间下才能获取较优的评估性能。因此，该情况下可以考虑利用图像中高维的视觉特征取代解剖特征，并且取多个特征对房角闭合状态进行全面描述，以获得更为准确的房角分类结果。由此可见，实现高效且客观的房角闭合程度计算机自动评估方法，是图像处理分析技术与眼前房角OCT影像相结合的有重要意义的课题。本发明基于多特征的眼前房角开放程度评估方法能够实现对前房角闭合状态的准确与快速分类。

发明内容

本发明考虑到：(1)传统的前房角评估方法，均基于对眼前房角内一些解剖结构相关特征参数的测量，但解剖结构细节在前房角OCT图像中并非全部能够被辨识，且1-2个解剖结构参数无法完全描述房角闭合状态信息；(2)由于OCT成像质量受噪声或眼动影响较大，易造成前房角隐窝区域细节不明显，因此在成像时受到扰动(如眼颤)或帧平均增强算法失效时，已有的基于解剖结构特征参数测量或者基于图像视觉特征的方法无法正确地对房角闭合状态进行分类。(3)已有的基于图像视觉特征信息的分类方法，特征选择较为单一，对房角隐窝区域描述不够全面，即使对于成像清晰的图像也仅能从大类上简单区分青光眼的种类，无法实现眼前房角的细化分类。因此，如何能够在适用范围更广的前房角FD-OCT图像中实现基于多重特征的房角开放程度分类方法具有较高的研究意义与实用价值。本发明基于局部点、局部线、局部区域等涉及亮度、边缘与区域形状等三个方面的图像特征，通过采用支持向量机分类能够实现对前房角闭合状态的客观且准确地分类方法。

本发明基于多特征的眼前房角开放程度分级评估方法，包括以下步骤：

步骤(1)对眼前房角OCT图像进行图像分割操作，分割出完整的眼前房角房角区域。

步骤(2)根据眼前房角结构特性，对分割后的图像做边缘检测，并根据角膜后边缘与虹膜前边缘的交点，找出房角隐窝区域作为感兴趣区域。

步骤(3)在感兴趣区域中分别计算梯度、Edgelet和CoHOG特征。

步骤(4)利用支持向量机对图像进行多级分类。

本发明的有益效果：

1、对于成像质量较低的图像具有较高的自动评级正确率，能够对眼前房角闭合状态进行细致准确的分类。

2、本发明充分考虑了前房角OCT图像亮度、边缘与区域形状等多个方面的图像特征，自动评级结果性能稳定，不易受成像时外界环境变化的影响，能够更快且客观地获取眼前房角的开放/闭合状态信息。

附图说明

图1(a)原图像I；

图1(b)分割后图像I'_BW；

图2感兴趣区域提取示意图；

图3方向量化图；

图4CoHOG矩阵的偏移量坐标，(中间圆圈为当前像素点，其它圆圈为其计算CoHOG时所参考的偏移坐标位置)。

具体实施方式

以下结合附图对本发明作进一步说明。

步骤(1)对眼前房角OCT图像进行图像分割操作，分割出完整的眼前房角区域，其步骤如下：

1-1、对图像I(图1(a))采用7×3的中值滤波器去除前房角OCT图像中的斑点噪声，得到图像P.

1-2、对图像P做形态学灰度重建得到I_R，其计算范围为各像素点的四邻域。

1-3、对I_R进行幂次变换得到图像F，其中将输入灰度值范围为[0.06,0.9]的像素映射至输出灰度值范围[0,1]，幂次参数γ＝0.18。

1-4、采用最大类间方差法对图像F进行阈值分割，得到分割后的图像I'_BW(图1(b))。

步骤(2)根据前房角结构特性，对分割后的图像做边缘检测，并根据角膜后边缘与虹膜前边缘的交点，找出房角隐窝部分作为感兴趣区域，该感兴趣区域用G(i,j)表示，其步骤如下：

2-1、从图像I'_BW顶端到底部逐行扫描，检测得到第一条初始边缘，即为角膜前边缘。

2-2、从I'_BW底部向上至第一条边缘逐行扫描检测第三条边缘，即为虹膜前边缘。

2-3、在I'_BW中从第三条边缘向上至第一条边缘逐行扫描，得到第二条边缘，即为角膜后边缘。

2-4、对检测到的角膜后边缘和虹膜前边缘分别进行直线拟合，并求得两条拟合直线的交点[x,y],取[x-a,y-b,w,h]的区域作为感兴趣区域G(i,j)，其中(x-a,y-b)为感兴趣区域左上角坐标，(w,h)为感兴趣区域宽度和高度，结果如图2所示。

步骤(3)在感兴趣区域中分别计算梯度、Edgelet和CoHOG特征，具体步骤如下：

计算图像中感兴趣区域的梯度值：

A、采用Sobel梯度算子获得感兴趣区域G的梯度图I_g。

B、在I_g中分别计算出各像素点在{0°、45°、90°和135°}四个方向上的梯度强度，并求其均值。

计算图像中感兴趣区域的Edgelet特征值：

为了计算Edgelet特征值，需采用一边长为d(d为奇数)的方形窗口w在I_g内滑动，令w内的像素坐标为(x,y)，x＝1,2,…,d，y＝1,2,…,d；

A、令方形窗口w在I_g中已检测到的角膜后边缘与虹膜前边缘位置上逐点滑动，并计算其中心点c(x,y)处的法线方向角θ_E(x,y)；由于d取值较小，故该计算可近似为计算w内边缘两端所连接得到的直线与纵坐标的交角。并将θ_E(x,y)根据图3中方向量化图所示量化为6个角度，并分别以0-5这六个数字表示，得到角β_E(x,y)。

B、根据I_g中c(x,y)点处的纵向和横向梯度值，计算出该点处的梯度方向角θ_I(x,y)。并将θ_I(x,y)根据图3中方向量化图所示量化为6个角度，并分别以0-5这六个数字表示，得到角β_I(x,y)。

C、由公式(1)可计算出感兴趣区域G中，边缘上点c(x,y)处的edgelet特征值f

$f(E,G,v)=\frac{1}{k}\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{M}^G(u_i+v)|<n^G(u_i+v),n_i^E>|---\left(1\right)$

其中，k为选取的边缘片段E的长度，k＝1,2,…,d；u_i为方形窗口w中的坐标值(x,y),x＝1,2,…,d，y＝1,2,…,d；M^G为I_g中像素点c(x,y)的亮度；v为方形窗口w在原图像中的坐标位移(x',y')；n^G为像素点(u_i+v)的梯度向量；为像素点c的法向量；两个向量的点积可由式(2)近似表示，即可由式(3)得到点乘结果，从而计算得到最终的edgelet特征值f。

$|<n^I(u_i+w),n_i^E>|=L({\mathrm{\&beta;}}_I(u_i+w)-{\mathrm{\&beta;}}_E)---\left(2\right)$

$L\left[n\right]=\left\{\begin{array}{cc}1& n=0\\ \frac{4}{5}& n=\&PlusMinus;1,\&PlusMinus;5\\ \frac{1}{2}& n=\&PlusMinus;2,\&PlusMinus;4\\ 0& n=\&PlusMinus;3\end{array}\right.---\left(3\right)$

计算图像中感兴趣区域的CoHOG特征值：

A、根据I_g中c(x,y)点处的纵向和横向梯度值，计算出该点处的梯度方向角θ_I(x,y)。将G中各像素点的θ_I(x,y)进行量化得到α(x,y)。其中，θ_I(x,y)按每45°进行量化，可以得到8个角度，将其分别表示为α(x,y)，α(x,y)取值为0,1,2,…,7。

B、将G分割为多个不互相重叠的边长为N×M个像素的区域。

C、在每一个子区域内计算共生矩阵C，由步骤a得到梯度方向角已量化为8个角度，那么任何两个梯度方向m和n形成一个角度对，两两组合就有64种组合。同时，计算共生矩阵时还要考虑位置偏移(x_d,y_d)，即一个角度对内的两个梯度方向角的位置关系。于是，在每个子区域内根据式(4)可计算共生矩阵C；

$C(x_d,y_d,m,n)=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{q=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\begin{array}{cc}1,& \mathrm{ifG}(p,q)=\mathrm{mandG}(p+x_d,q+y_d)=n\\ 0,& \mathrm{otherwise}.\end{array}\right.---\left(4\right)$

其中，(p,q)为子区域内的像素坐标，图4中白色圆圈所示；(m,n)为共生矩阵C内的坐标，m和n取值分别为8个角度，即0,1,…,7；(x_d,y_d)为像素点(p,q)的偏移坐标，其取值范围如图4中灰色圆圈所示位置。

D、将每个子区域的共生矩阵C串接构成感兴趣区域的CoHOG特征向量。

步骤(4)采用支持向量机(support vector machine,SVM)根据以上三个特征向量对图像进行分类，从而实现对房角开放/闭合程度的评估。

Claims (1)

1.基于OCT影像的眼前房角开放程度多特征分类方法，其特征在于该方法包括如下步骤：

步骤(1)对眼前房角OCT图像进行图像分割操作，分割出完整的眼前房角房角区域；具体是：

1-1、对眼前房角OCT图像采用7×3的中值滤波器去除前房角OCT图像中的斑点噪声，得到图像P；

1-2、对图像P做形态学灰度重建得到I_R，其计算范围为各像素点的四邻域；

1-3、对I_R进行幂次变换得到图像F，其中将输入灰度值范围为[0.06,0.9]的像素映射至输出灰度值范围[0,1]；

1-4、采用最大类间方差法对图像F进行阈值分割，得到分割后的图像I'_BW；

步骤(2)根据眼前房角结构特性，对分割后的图像做边缘检测，并根据角膜后边缘与虹膜前边缘的交点，找出房角隐窝区域作为感兴趣区域，具体是：

2-1、从图像I'_BW顶端到底部逐行扫描，检测得到第一条初始边缘，即为角膜前边缘；

2-2、从图像I'_BW底部向上至第一条边缘逐行扫描检测第三条边缘，即为虹膜前边缘；

2-3、在图像I'_BW中从第三条边缘向上至第一条边缘逐行扫描，得到第二条边缘，即为角膜后边缘

2-4、对检测到的角膜后边缘和虹膜前边缘分别进行直线拟合，并求得两条拟合直线的交点[x,y],取[x-a,y-b,w,h]的区域作为感兴趣区域G(i,j)，其中(x-a,y-b)为感兴趣区域左上角坐标，(w,h)为感兴趣区域宽度和高度；

步骤(3)在感兴趣区域中分别计算梯度值、Edgelet和CoHOG特征值，具体是：

梯度值的计算方法为：

3-1-1、采用Sobel梯度算子获得感兴趣区域G的梯度图I_g；

3-1-2、在梯度图I_g中分别计算出各像素点在{0°、45°、90°和135°}四个方向上的梯度强度，并求其均值；

Edgelet特征值的计算方法为：

采用一边长为d的方形窗口w在梯度图I_g内滑动，d为奇数，令方形窗口w内的像素坐标为(x,y)，x＝1,2,…,d，y＝1,2,…,d；

3-2-1、令方形窗口w在梯度图I_g中已检测到的角膜后边缘与虹膜前边缘位置上逐点滑动，并计算其中心点c(x,y)处的法线方向角θ_E(x,y)；将θ_E(x,y)量化为6个角度，并分别以0-5这六个数字表示，得到角β_E(x,y)；

3-2-2、根据梯度图I_g中c(x,y)点处的纵向和横向梯度值，计算出该点处的梯度方向角θ_I(x,y)；并将θ_I(x,y)量化为6个角度，并分别以0-5这六个数字表示，得到角β_I(x,y)；

3-2-3、由公式(1)可计算出感兴趣区域G中，边缘上点c(x,y)处的Edgelet特征值f

$f(E,G,v)=\frac{1}{k}\underset{i=1}{\overset{k}{\mathrm{\&Sigma;}}}{M}^G(u_i+v)|<n^G(u_i+v),n_i^E>|---\left(1\right)$

其中，k为选取的边缘片段E的长度，k＝1,2,…,d；u_i为方形窗口w中的坐标值(x,y),x＝1,2,…,d，y＝1,2,…,d；M^G为I_g中像素点c(x,y)的亮度；v为方形窗口w在原图像中的坐标位移(x',y')；n^G为像素点(u_i+v)的梯度向量；为像素点c的法向量；两个向量的点积可由式(2)近似表示，即可由式(3)得到点乘结果，从而计算得到最终的Edgelet特征值f；

$|<n^I(u_i+w),n_i^E>|=L({\mathrm{\&beta;}}_I(u_i+w)-{\mathrm{\&beta;}}_E)---\left(2\right)$

$L\left[n\right]=\left\{\begin{array}{cc}1& n=0\\ \frac{4}{5}& n=\&PlusMinus;1,\&PlusMinus;5\\ \frac{1}{2}& n=\&PlusMinus;2,\&PlusMinus;4\\ 0& n=\&PlusMinus;3\end{array}\right.---\left(3\right)$

CoHOG特征值的计算方法为：

3-3-1、根据梯度图I_g中c(x,y)点处的纵向和横向梯度值，计算出该点处的梯度方向角θ_I(x,y)；将感兴趣区域G中各像素点的θ_I(x,y)进行量化得到α(x,y)；其中，θ_I(x,y)按每45°进行量化，得到8个角度，将其分别表示为α(x,y)，α(x,y)取值为0,1,2,…,7；

3-3-2、将感兴趣区域G分割为多个不互相重叠的边长为N×M个像素的区域；

3-3-3、在每一个子区域内计算共生矩阵C，由得到梯度方向角已量化为8个角度，那么任何两个梯度方向m和n形成一个角度对，两两组合就有64种组合；同时，计算共生矩阵时考虑位置偏移(x_d,y_d)，即一个角度对内的两个梯度方向角的位置关系；于是，在每个子区域内根据式(4)可计算共生矩阵C；

$C(x_d,y_d,m,n)=\underset{p=1}{\overset{N}{\mathrm{\&Sigma;}}}\underset{q=1}{\overset{M}{\mathrm{\&Sigma;}}}\left\{\begin{array}{cc}1,& \mathrm{ifG}(p,q)=\mathrm{mandG}(p+x_d,q+y_d)=n\\ 0,& \mathrm{otherwise}.\end{array}\right.---\left(4\right)$

其中，(p,q)为子区域内的像素坐标；(m,n)为共生矩阵C内的坐标，m和n取值分别为8个角度，即0,1,…,7；(x_d,y_d)为像素点(p,q)的偏移坐标；

3-3-4、将每个子区域的共生矩阵C串接构成感兴趣区域的CoHOG特征向量；

步骤(4)根据步骤(3)得到的三个特征值，再利用支持向量机对图像进行多级分类。