CN104143087A - 基于sd-oct的挫伤性视网膜内外节缺失检测方法 - Google Patents

Abstract

本发明公开了一种基于SD-OCT图像的挫伤性视网膜内节/外节缺失三维自动检测方法，包括以下步骤：（1）图像预处理，采用多尺度三维图割方法将视网膜内部自动分割成11个表面，提取位于第7表面和第8表面之间的内节/外节区域作为感兴趣区域，进行平坦化处理和双线性滤波增强；（2）对感兴趣区域中每个体素提取5类（共计57个）特征；（3）采用主成份分析方法对特征进行优化选择；（4）将特征样本分成训练集和测试集，对训练集样本采用Adaboost算法训练集成分类器；（5）对测试集样本进行缺失/非缺失识别；（6）对识别结果进行血管轮廓影响排除和孤立点排除等后处理，计算相应的缺失体积，缺失体积识别误差小，准确性好。

Description

基于SD-OCT的挫伤性视网膜内外节缺失检测方法

技术领域

本发明属于视网膜图像处理与分析方法，尤其是对SD-OCT(频域光学相干断层成像)视网膜图像中的内节/外节缺失的三维自动定量检测方法。

背景技术

眼外伤是导致视力损伤甚至致盲的主要因素之一，是我国单眼盲的首位致盲原因。当眼球表面的损伤传递到后端的视网膜时，视网膜震荡伤的特点是灰白变色，或是当外伤闭合后，视网膜仍然浑浊。人眼的组织病理学研究和动物研究表明，光感受器的内节/外节连接区损伤是视网膜震荡伤的发病原。而光感受器的内节/外节与光传导紧密相关，因此，内节/外节的完整性与视力密切相关。

SD-OCT(频域光学相干断层成像)是一种无损评估视网膜生理和病理指标的强有力技术，能快速获得高分辨的断层三维图像。在SD-OCT图像中，光感受器内节/外节是介于外界膜和色素上皮层之间的一个高亮反射层。内节/外节完整性的遭破坏程度是光感受器受损程度的一个代表性指标，通常与视网膜震荡伤导致的视力低下以及其他视网膜疾病有关。

在SD-OCT图像中，通过定量检测内节/外节区域缺失范围来定量评估光感受器的损伤程度是一个非常有意义的研究方向。有文献报道在一个B扫描面中进行内节/外节缺失长度的手动检测。但是，这是基于单个二维界面图像。此外，有文献提出采用部分OCT投影图像或者冠状面投影图像等技术更好地可视化光感受器的整体性和缺失区域。但是，这些检测内节/外节缺失面积的技术仍然是基于二维图像的，而且是手动的，这在选择缺失区域边界时可能会引入主观因素。关于内节/外节缺失区域边界的自动检测和内节/外籍缺失面积计算方法也有文献报道，但是，这种方法是基于缺失区域是圆形的假设，而实际上，缺失区域的形状可能是任意形状。

更重要的是，内节/外节是一层有厚度的薄膜，应该以内节/外节缺失的体积来衡量光感受器的损伤程度。这方面的研究还处于起步阶段，现有技术对视网膜内外节缺失检测也仅仅只针对局部较小的范围进行检测，特征提取不全面，并且分类不平衡，例如，CN103679198A公开了一种基于K近邻的视网膜内节外节缺失自动检测，提取以黄斑中心为中心、直径1mm的内节外节区域为感兴趣区域，进行5个特征提取，提取特征少，特征不够典型，且由于K近邻分类器的约束，采用的K近邻分类存在分类不平衡问题，分类性能差；其次，在后处理与缺失体积计算中，对识别分类得到的结果没有进行血管轮廓影响排除和孤立点排除，缺失体积误差大，准确性差。

发明内容

本发明克服了现有技术的不足，提供了一种具有可行性和有效性的视网膜内节/外节缺失三维自动定量检测方法，能够将内节/外节全部区域作为感兴趣区域，全面提取体素特征，充分考虑不平衡分类问题的特点，精确计算视网膜内外节缺失区域。

为解决上述技术问题，本发明采用的技术方案为：

一种基于SD-OCT的挫伤性视网膜内外节缺失检测方法，该方法主要包括6个步骤：图像预处理，体素特征提取，特征选择，集成分类器的训练，体素缺失识别，后处理与缺失体积计算。

一种基于SD-OCT的挫伤性视网膜内外节缺失检测方法，即包括以下步骤：

(1)图像预处理：采用多尺度三维图割方法进行视网膜内部分层，将视网膜SD-OCT图像自动分割成10个分层，产生11个表面，将包含内节/外节区域的第7层和第8层之间的数据提取出来，参考第11层进行平坦化，将平坦化后第7层和第8层之间的内节/外节区域的SD-OCT图像作为感兴趣区域，进行双线性滤波增强；

(2)体素特征提取：对步骤(1)所述平坦化后的内节/外节区域的所有体素的灰度值进行归一化处理，根据手动标记的金标准，对内节/外节区域内的缺失体素和非缺失体素分别提取5类特征，所述5类特征包括归一化的灰度值、分块均值、分块标准差、13个方向的灰度差绝对值、基于灰度共生矩阵的特征；

(3)特征优化选择：对步骤(2)中所提取的特征采用主成份分析法进行特征选择；

(4)集成分类器的训练：将经步骤(3)选择的特征分成训练集和测试集；在训练集中，按1：1比例，从非缺失样本中随机抽取数量与缺失样本数量相同的特征样本，采用Ababoost算法进行集成分类器的训练；

(5)体素缺失/非缺失的分类识别：对测试集样本进行缺失/非缺失识别，采用步骤(4)训练得到的集成分类器对测试样本进行分类，对每个体素的缺失与否进行分类识别；

(6)后处理与缺失体积计算：对步骤(5)分类识别得到的缺失结果进行血管轮廓影响排除和孤立点排除，并计算相应的缺失体积。

本发明技术方案具体描述如下。

(1)图像预处理

图像预处理主要包括以下三个步骤：视网膜内部分层，包含内节/外节区域的感兴趣区域提取和平坦化，以及图像的双线性滤波增强。

(1-1)视网膜内部分层

人类视网膜是一个相当薄的组织，厚度不到0.5mm，是眼睛最重要的组成部分。视网膜本身具有相当复杂的结构，基本结构类似于一个“三层蛋糕”——三层神经细胞被两层类神经连接层隔开，可以更精细地分成10个层。视网膜内部分层对于分析视网膜病变如眼外伤的严重程度、黄斑水肿的形成等具有重要作用。采用多尺度三维图搜索方法将视网膜SD-OCT图像自动分割成10个分层，产生11个表面。所述10个分层自上至下依次为：神经纤维层，神经节细胞层，内丛状层，内核层，外丛状层，外核层，内节层，外节层，维尔赫夫氏膜和视网膜色素上皮层；

三维图搜索方法是基于图论的方法，采用从粗糙到精细的不同分辨率，来检测视网膜的内部各表面。图搜索采用基于边界的代价函数，当代价函数最小时，即找到各个表面；

(1-2)内节/外节区域提取和平坦化

完成了各分层的分割后，可以将包含内节/外节区域的第7层和第8层之间的数据提取出来，并参考视网膜色素上皮层的下界即第11层进行平坦化，因为第11层的分层效果相对鲁棒。本发明主要平坦化后的第7层与第8层之间的SD-OCT图像作为感兴趣区域，研究内节/外节的缺失检测。即将包含内节/外节区域的第7层和第8层之间的数据提取出来，参考第11层进行平坦化，将平坦化后第7层和第8层之间的内节/外节区域的SD-OCT图像作为感兴趣区域，进行双线性滤波增强。

(2)特征提取

本发明是一种基于纹理特征分类的内节/外节缺失检测方法，首先，对感兴趣区域内的所有体素的灰度值进行归一化处理，然后提取了5类特征，5类特征包括归一化的灰度值、分块均值、分块标准差、13个方向的灰度差绝对值、基于灰度共生矩阵的特征；5类特征共计57个特征：归一化的灰度值(特征1)，分块均值(特征2)，分块标准差(特征3)，13个方向的灰度差绝对值(步长分别取1和2，特征4-5)，基于灰度共生矩阵的特征，基于灰度共生矩阵的特征包括对比度(13个方向，特征6-18)、相关性(13个方向，特征19-31)、能量(13个方向，特征32-44)、均匀性(13个方向，特征45-57)。

(2-1)归一化灰度值表示每个体素的灰度级，对感兴趣区域内的所有体素的灰度值采用如式(1)所示的方法进行线性归一化至0～255。

$I_{\mathrm{normalized}}(i,j,k)=\frac{I_{\mathrm{original}}(i,j,k)-I_{\min }}{I_{\max }-I_{\min }}\×255---\left(1\right)$

其中，(i,j,k)表示体素的坐标，I_normalized(i,j,k)表示体素(i,j,k)归一化后的灰度值，I_original(i,j,k)表示体素(i,j,k)原始的灰度值，I_min和I_max分别表示感兴趣区域内体素灰度值的最小值和最大值。内节/外节缺失区域的一个重要特点是灰度值低于非缺失区域，即内节/外节缺失区域的灰度值低于非缺失区域的灰度值。

(2-2)分块均值用于表示以体素(i,j,k)为中心的分块的灰度均值，5×5×5分块的灰度均值M_block(i,j,k)，计算公式如式(2)：

$M_{\mathrm{block}}(i,j,k)=\frac{1}{125}\underset{l=i-2}{\overset{i+2}{\mathrm{\Σ}}}\underset{m=j-2}{\overset{j+2}{\mathrm{\Σ}}}\underset{n=k-2}{\overset{k+2}{\mathrm{\Σ}}}{I}_{\mathrm{normalized}}(l,m,n)---\left(2\right)$

由于内节/外节缺失区域的灰度值整体比非缺失区域的灰度值小，所以内节/外节缺失区域的分块均值M_block(i,j,k)小于非缺失区域的分块均值。

(2-3)分块标准差STD_block(i,j,k)用于表示以体素(i,j,k)为中心的5×5×5分块的灰度值变化情况，表征灰度值在5×5×5分块中的分散程度；可根据分块标准差STD_block(i,j,k)识别分块体素全部位于缺失区域、分块体素全部位于非缺失区域、或者分块中同时包含缺失区域和非缺失区域的体素，计算公式如公式(2)：

${\mathrm{STD}}_{\mathrm{block}}(i,j,k)=\sqrt{\frac{\underset{l=i-2}{\overset{i+2}{\mathrm{\Σ}}}\underset{m=j-2}{\overset{j+2}{\mathrm{\Σ}}}\underset{n=k-2}{\overset{k+2}{\mathrm{\Σ}}}{(I_{\mathrm{normalized}}(l,m,n)-M_{\mathrm{block}}\left(\mathrm{imjmk}\right))}^2}{124}}---\left(3\right)$

依据分块标准差STD_block(i,j,k)的值识别分块体素位于缺失区域或者非缺失区域或者同时包含缺失区域和非缺失区域的体素；如果某分块中的体素灰度值都很接近，即该分块全部位于缺失区域或者全部位于非缺失区域，即STD_block(i,j,k)的值属于一个经验值偏小的范围，STD_block(i,j,k)值的范围与图像本身体素灰度值有关，例如小于所有体素灰度值STD_block(i,j,k)平均值；如果分块中同时包含缺失区域和非缺失区域的体素，则对应的STD_block(i,j,k)的值属于一个经验值偏大的范围，例如，大于所有体素灰度值STD_block(i,j,k)平均值。通过STD_block(i,j,k)的值识别分块体素位于缺失区域或者非缺失区域或者同时包含缺失区域和非缺失区域的体素，STD_block(i,j,k)的值识别范围没有一个明确的范围，基于操作者的经验值识别，一般如果某个分块的STD_block(i,j,k)小于所有体素灰度值STD_block(i,j,k)平均值，分块全部位于缺失区域或者全部位于非缺失区域，否则，分块中同时包含缺失区域和非缺失区域的体素；

(2-4)13个方向上的灰度差绝对值之和定义为：中心点体素(i,j,k)的灰度值与周围体素(13个方向)的灰度差绝对值之和，13个方向上的灰度差绝对值反映中心点体素在缺失区域或非缺失区域的位置，依据所述13个方向上的灰度差绝对值之和AID_step(i,j,k)判断中心点体素(i,j,k)位于缺失与非缺失交界处或者缺失区域或者非缺失区域；若该特征值属于一个比较大的范围，一般大于或者等于0.5，说明该中心点体素位于缺失与非缺失交界处；若该特征值接近于零，一般小于0.5，说明中心点体素位于缺失区域或者非缺失区域，这里选择了两种距离步长，step＝1，2，即step为1或者2，具体计算公式如下：

${\mathrm{AID}}_{\mathrm{step}}(i,j,k)=\underset{13\mathrm{directions}}{\mathrm{\Σ}}|I_{\mathrm{normalized}}(i,j,k)-I_{\mathrm{normalized}}(l,m,n)|,\mathrm{step}=\mathrm{1,2}---\left(4\right)$

其中所述13个方向的选取方法为：假设α₁表示X轴与X-Y平面投影图的夹角，α₂表示Z轴与X-Y平面投影图的夹角，那么所述13个方向分别为(α₁,α₂)＝(0,90°)、(45°,90°)、(90°,90°)、(135°,90°)、(0,45°)、(180°,45°)、(90°,45°)、(-90°,45°)、(0,0)、(45°,45°)、(135°,45°)、(-45°,45°)、(-135°,45°)；13个方向对应的l、m、n取值分别为(i+step,j,k)、(i+step,j+step,k)、(i,j+step,k)、(i-step,j+step,k)、(i+step,j,k+step)、(i-step,j,k+step)、(i,j+step,k+step)、(i,j-step,k+step)、(i,j,k+step)、(i+step,j+step,k+step)、(i-step,j+step,k+step)、(i+step,j-step,k+step)、(i-step,j-step,k+step)；表示13个方向上灰度差绝对值之和。

(2-5)基于灰度共生矩阵的特征描述多个立体切面之间灰度级的空间相关性，基于灰度共生矩阵的特征为三维体数据的灰度共生矩阵，采用了一种三维方法在13个方向上跨切面搜索不同灰度级，可以构造13个灰度共生矩阵；13个灰度共生矩阵构造方法具体包括以下步骤：

对每个5×5×5分块构造所述13个不同方向(步骤(2-4)中选取的13个方向)上的灰度共生矩阵；然后，分别提取4个基于该灰度共生矩阵的特征：(a)对比度，用于度量立体图像的局部对比度，当大灰度差情况发生概率大时，对比度值较大；(b)相关性，用于表征体素对中两个体素的相关性，当体素对的灰度值越相关时，该相关性特征值越大；(c)能量，用于度量体素对的重复出现次数，当体素对重复次数越多时，能量值越大；(d)均匀性，用于表征体素对的局部均匀性，当每个体素对的灰度级越相近时，该均匀性特征值越大。

对比度Contrast_direction m(i,j,k)、相关性Correlation_direction m(i,j,k)、能量Energy_direction m(i,j,k)和均匀性Homogeneity_direction m(i,j,k)分别采用式(5)、式(6)、式(7)和式(8)进行计算：

${\mathrm{Contrast}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}{(x-y)}^{2}p(x,y),m=\mathrm{1,2},...,13---\left(5\right)$

${\mathrm{Correlation}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}\frac{(x-{\mathrm{\μ}}_x)(y-{\mathrm{\μ}}_y)p(x,y)}{{\mathrm{\σ}}_x{\mathrm{\σ}}_y},m=\mathrm{1,2},...,13---\left(6\right)$

${\mathrm{Energy}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}{p}^2(x,y),m=\mathrm{1,2},...,13---\left(7\right)$

${\mathrm{Homogeneity}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}\frac{p(x,y)}{1+|x-y|},m=\mathrm{1,2},...,13---\left(8\right)$

其中，p(x,y)表示灰度共生矩阵中第(x,y)个元素；N_g表示量化后的灰度级(一般基于经验值选取量化后的主要灰度级)，取N_g＝8；μ_x和σ_x分别表示第x行的均值和标准差；μ_y和σ_y分别表示第y列的均值和标准差。

(3)特征优化选择

特征选择目的是降低特征空间的维数，从大的特征集中选择主要的特征，对于训练样本，本发明采用主成份分析方法进行特征选择。对于步骤(2)所述提取的5类特征采用主成份分析方法进行特征选择，将经过特征选择后的结果作为新的特征，根据新特征所代表的数据总方差的比例大于等于90％的原则，选择前12个主成分作为新特征；

完成特征优化选择后，感兴趣区域内的每个体素用一个12维的特征向量来表征，步骤(4)的集成分类器将在特征空间中进行。

(4)集成分类器的训练

由于缺失体素的数量远少于非缺失体素的数量，两类样本的类间分布不平衡，属于典型的非平衡分类问题。传统的分类器，大都建立在训练样本集分布基本均衡的前提下，且以总体分类精度作为目标，在本问题中，由于缺失类样本相对稀缺以及不可避免的噪声问题，传统的分类器如支持向量机、贝叶斯分类器、分类与回归树方法、决策树、线性判别分析、双线性判别分析等，都会出现分类器对数量上占优势的多数类(非缺失类)存在明显的偏好，使得多数类的分类精度提高，而少数类(缺失类)的分类精度降低。本发明通过下采样方式舍弃大部分多数类样本，以缺失样本数量与非缺失样本数量的比例为1：1比例构造新的训练集，并采用分类器集成的方法来进一步克服不平衡分类问题，以决策树为弱分类器，采用Adaboost算法进行集成分类器的训练。

(5)体素缺失/非缺失的分类识别

采用步骤(4)训练得到的集成分类器，对内节/外节区域内的体素进行分类识别，分成缺失/非缺失两类。

(6)后处理与缺失体积计算

由于视网膜血管轮廓的影响，SD-OCT图像的B扫描图的相应位置上会出现上下贯穿的阴影，在IS/OS连接区内也会出现血管阴影，这些阴影往往会被误判为缺失，因此需要进行排除。本设计首先将三维OCT数据投影至X-Y平面，并对投影图进行血管轮廓检测，若步骤(5)得到的缺失识别结果的X坐标值与血管轮廓的X坐标值相等，则认为是缺失误识别，将对应的识别结果改为非缺失，可有效地排除血管轮廓对内节/外节缺失检测的影响。

由于内节/外节缺失的生理特点，孤立的缺失点或者孤立的非缺失点应该不存在，因此，对检测结果采用形态学方法进行了孤立点的排除。

经后处理的检测结果是以体素个数为单位的，最后检测结果乘以体素的分辨率，转化为检测到的内节/外节缺失体积。

步骤(6)具体包括以下步骤：

(6-1)血管轮廓影响的排除：将第1至第10层的全部三维SD-OCT数据投影至X-Y平面，对投影图进行血管轮廓检测，将步骤(5)的缺失识别结果的X坐标值等于血管轮廓的X坐标值，则认为是缺失误识别，将对应的识别结果修改为非缺失区域；

(6-2)将步骤(6-1)的缺失识别结果采用形态学方法进行孤立点的排除；

(6-3)将步骤(6-2)的孤立点的排除后缺失识别结果乘以体素的分辨率，转化为检测到的内节/外节缺失体积。

本发明有益效果包括：

本发明提供了一种基于SD-OCT图像、具有可行性和有效性的视网膜内节/外节缺失三维自动检测方法，将感兴趣区域扩展至整个内节/外节连接区；提取57个体素特征，采用了主成份分析方法进行了特征选择与优化，充分考虑到了不平衡分类问题的存在，采用了Adaboost算法训练集成分类器，从测试结果来看，分类性能比之前申请有较大提高，视网膜内节外节检测准确，同时，本发明对识别分类得到的结果进行血管投影和孤立点排除，缺失体积识别误差小，准确性好；

进一步地，本申请扩大了训练和测试的数据集，包括15个外伤眼和15个正常眼；对于眼外伤等临床常见眼科疾病的诊断与治疗起到了重要的辅助作用。

附图说明

图1为本发明的内节/外节缺失三维自动检测的整体框图；

图2为图2视网膜的分层效果图(SD-OCT图像的X-Z平面图)，分成了10个层(11个表面)，自上至下依次为：神经纤维层，神经节细胞层，内丛状层，内核层，外丛状层，外核层，内节层，外节层，维尔赫夫氏膜和视网膜色素上皮层；

图3为平坦化后的SD-OCT图像和感兴趣区域(X-Z平面图)；

图4为双线性滤波后的SD-OCT图像和感兴趣区域(X-Z平面图)；

图5为血管轮廓的影响，(a)箭头所示为血管轮廓在OCT图像X-Z平面图中的位置；(b)视网膜OCT图像在X-Y平面上的投影图以及血管；(a)和(b)图中箭头处的血管位置一一对应；

图6内节/外节缺失检测结果与金标准。图(a)表示由观察者1手动标记的一个外伤眼数据的内节/外节缺失三维金标准，图(b)表示对该外伤眼数据采用本发明方法检测到的内节/外节缺失区域；图(c)表示采用本发明方法对正常眼数据的检测结果；

图7为外伤数据和正常数据缺失体积检测结果的比较。

具体实施方式

下面结合附图对本发明作更进一步的说明。

本方法的基本框图如附图1所示，主要包括6个步骤：图像预处理，体素特征提取，特征优化选择，集成分类器训练，体素缺失/非缺失分类识别，后处理与缺失体积计算。具体描述如下。

(1)图像预处理

图像预处理主要包括以下三个步骤：视网膜内部分层，包含内节/外节区域的感兴趣区域的提取和平坦化，图像的双线性滤波增强。

(a)视网膜内部分层

人类视网膜是一个相当薄的组织，厚度不到0.5mm，是眼睛最重要的组成部分。视网膜本身具有相当复杂的结构，基本结构类似于一个“三层蛋糕”——三层神经细胞被两层类神经连接层隔开，可以更精细地分成10个层。视网膜内部分层对于分析视网膜病变如眼外伤的严重程度、黄斑水肿的形成等具有重要作用。采用多尺度三维图搜索方法将视网膜SD-OCT图像自动分割成10个分层，产生11个表面，如附图2所示。三维图搜索方法是基于图论的方法，采用从粗糙到精细的不同分辨率，来检测视网膜的内部各表面。图搜索采用基于边界的代价函数，当代价函数最小时，即找到各个表面。

(b)内节/外节区域提取和平坦化

完成了各分层的分割后，可以将包含内节/外节区域的第7层和第8层之间的数据提取出来，并参考视网膜色素上皮层(第10层)的下界即第11层进行平坦化，因为第11层的分层效果相对鲁棒。本发明主要平坦化后的第7层与第8层之间的SD-OCT图像作为感兴趣区域，研究内节/外节的缺失检测。平坦化后的SD-OCT图像和感兴趣区域(X-Z平面图)如附图3所示，双线性滤波后的SD-OCT图像和感兴趣区域(X-Z平面图)如图4所示，滤除图像中的散斑噪音。

(2)体素特征提取

本发明是一种基于纹理特征分类的内节/外节缺失检测方法。首先，对感兴趣区域内的所有体素的灰度值进行归一化处理，然后提取了5类特征，共计57个特征：归一化的灰度值(特征1)，分块均值(特征2)，分块标准差(特征3)，13个方向的灰度差绝对值(步长分别取1和2，特征4-5)，基于灰度共生矩阵的特征，包括对比度(13个方向，特征6-18)，相关性(13个方向，特征19-31)，能量(13个方向，特征32-44)，均匀性(13个方向，特征45-57)。

归一化灰度值表示每个体素的灰度级，对感兴趣区域内的所有体素的灰度值采用如式(1)所示的方法进行线性归一化至0～255。

$I_{\mathrm{normalized}}(i,j,k)=\frac{I_{\mathrm{original}}(i,j,k)-I_{\min }}{I_{\max }-I_{\min }}\×255---\left(1\right)$

其中，(i,j,k)表示体素的坐标，I_normalized(i,j,k)表示体素(i,j,k)归一化后的灰度值，I_original(i,j,k)表示体素(i,j,k)原始的灰度值，I_min和I_max分别表示感兴趣区域内体素灰度值的最小值和最大值。内节/外节缺失区域的一个重要特点是灰度值低于非缺失区域。

分块均值用于表示以体素(i,j,k)为中心的小分块的灰度均值，这里计算了5×5×5小分块的灰度均值M_block(i,j,k)，计算公式如下：

$M_{\mathrm{block}}(i,j,k)=\frac{1}{125}\underset{l=i-2}{\overset{i+2}{\mathrm{\Σ}}}\underset{m=j-2}{\overset{j+2}{\mathrm{\Σ}}}\underset{n=k-2}{\overset{k+2}{\mathrm{\Σ}}}{I}_{\mathrm{normalized}}(l,m,n)---\left(2\right)$

由于内节/外节缺失区域的灰度值整体比非缺失区域的灰度值小，所以内节/外节缺失区域的分块均值M_block(i,j,k)小于非缺失区域的分块均值。

分块标准差STD_block(i,j,k)用于表示以体素(i,j,k)为中心的5×5×5小分块的灰度值变化情况，可以表征灰度值在5×5×5分块中的分散程度。可根据分块标准差识别分块体素全部位于缺失区域、分块体素全部位于非缺失区域、或者分块中同时包含缺失区域和非缺失区域的体素，计算公式如下：

${\mathrm{STD}}_{\mathrm{block}}(i,j,k)=\sqrt{\frac{\underset{l=i-2}{\overset{i+2}{\mathrm{\Σ}}}\underset{m=j-2}{\overset{j+2}{\mathrm{\Σ}}}\underset{n=k-2}{\overset{k+2}{\mathrm{\Σ}}}{(I_{\mathrm{normalized}}(l,m,n)-M_{\mathrm{block}}\left(\mathrm{imjmk}\right))}^2}{124}}---\left(3\right)$

13个方向上的灰度差绝对值之和AID_step(i,j,k)定义为：中心点体素(i,j,k)的灰度值与周围体素(13个方向)的灰度差绝对值之和，该特征可以很好地反映中心点体素在缺失区域或非缺失区域的大概位置，若该特征值大于0.5，说明该中心点体素位于缺失与非缺失交界处；若该特征值小于0.5，说明中心点体素位于较大块缺失区域或者非缺失区域，这里选择了两种距离步长，step＝1，2，具体计算公式如下：

${\mathrm{AID}}_{\mathrm{step}}(i,j,k)=\underset{13\mathrm{directions}}{\mathrm{\Σ}}|I_{\mathrm{normalized}}(i,j,k)-I_{\mathrm{normalized}}(l,m,n)|,\mathrm{step}=\mathrm{1,2}---\left(4\right)$

其中所述13个方向的选取方法为：假设α₁表示X轴与X-Y平面投影图的夹角，α₂表示Z轴与X-Y平面投影图的夹角，那么所述13个方向分别为(α₁,α₂)＝(0,90°)、(45°,90°)、(90°,90°)、(135°,90°)、(0,45°)、(180°,45°)、(90°,45°)、(-90°,45°)、(0,0)、(45°,45°)、(135°,45°)、(-45°,45°)、(-135°,45°)；13个方向对应的l、m、n取值分别为(i+step,j,k)、(i+step,j+step,k)、(i,j+step,k)、(i-step,j+step,k)、(i+step,j,k+step)、(i-step,j,k+step)、(i,j+step,k+step)、(i,j-step,k+step)、(i,j,k+step)、(i+step,j+step,k+step)、(i-step,j+step,k+step)、(i+step,j-step,k+step)、(i-step,j-step,k+step)；表示13个方向上灰度差绝对值之和。

三维体数据的灰度共生矩阵可以描述多个立体切面之间灰度级的空间相关性。本设计采用了一种三维方法在13个方向上跨切面搜索不同灰度级，可以构造13个灰度共生矩阵。假设α₁表示X轴与X-Y平面投影图的夹角，α₂表示Z轴与X-Y平面投影图的夹角，那么13个方向分别为(α₁,α₂)＝(0,90°)、(45°,90°)、(90°,90°)、(135°,90°)、(0,45°)、(180°,45°)、(90°,45°)、(-90°,45°)、(0,0)、(45°,45°)、(135°,45°)、(-45°,45°)、(-135°,45°)。这里，对每个5×5×5小分块构造了以上13个不同方向上的灰度共生矩阵。然后，分别提取4个基于该灰度共生矩阵的特征：(a)对比度，用于度量立体图像的局部对比度，当大灰度差情况发生概率大时，对比度值较大；(b)相关性，用于表征体素对中两个体素的相关性，当体素对的灰度值越相关时，该相关性特征值越大；(c)能量，用于度量体素对的重复出现次数，当体素对重复次数越多时，能量值越大；(d)均匀性，用于表征体素对的局部均匀性，当每个体素对的灰度级越相近时，该均匀性特征值越大。对比度Contrast_direction m(i,j,k)、相关性Correlation_direction m(i,j,k)、能量Energy_direction m(i,j,k)和均匀性Homogeneity_direction m(i,j,k)分别采用式(5)、式(6)、式(7)和式(8)进行计算。

${\mathrm{Contrast}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}{(x-y)}^{2}p(x,y),m=\mathrm{1,2},...,13---\left(5\right)$

${\mathrm{Correlation}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}\frac{(x-{\mathrm{\μ}}_x)(y-{\mathrm{\μ}}_y)p(x,y)}{{\mathrm{\σ}}_x{\mathrm{\σ}}_y},m=\mathrm{1,2},...,13---\left(6\right)$

${\mathrm{Energy}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}{p}^2(x,y),m=\mathrm{1,2},...,13---\left(7\right)$

${\mathrm{Homogeneity}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\Σ}}}\underset{y=1}{\overset{N_g}{\mathrm{\Σ}}}\frac{p(x,y)}{1+|x-y|},m=\mathrm{1,2},...,13---\left(8\right)$

其中，p(x,y)表示灰度共生矩阵中第(x,y)个元素；N_g表示量化后的灰度级，取N_g＝8；μ_x和σ_x分别表示第x行的均值和标准差；μ_y和σ_y分别表示第y列的均值和标准差。

(3)特征优化选择

特征选择是模式识别中常用的一种技术，目的是降低特征空间的维数，从大的特征集中选择主要的特征。本发明采用了主成份分析方法进行特征选择。根据新特征所代表的数据总方差的比例大于等于90％，选择前12个主成分作为新特征。

完成特征优化选择后，感兴趣区域内的每个体素可用一个12维的特征向量来表征，接下来的分类将在特征空间中进行。

(4)集成分类器的训练

由于缺失体素的数量远少于非缺失体素的数量，两类样本的类间分布不平衡，属于典型的非平衡分类问题。传统的分类器，大都建立在训练样本集分布基本均衡的前提下，且以总体分类精度作为目标，在本问题中，由于缺失类样本相对稀缺以及不可避免的噪声问题，传统的分类器如支持向量机、贝叶斯分类器、分类与回归树方法、决策树、线性判别分析、双线性判别分析等，都会出现分类器对数量上占优势的多数类(非缺失类)存在明显的偏好，使得多数类的分类精度提高，而少数类(缺失类)的分类精度降低。本发明通过下采样方式舍弃大部分多数类样本，以缺失样本数量与非缺失样本数量1：1的比例构造训练集，并采用分类器集成的方法来进一步克服不平衡分类问题，以决策树为弱分类器，采用Adaboost算法进行集成分类器的训练。

(5)体素缺失/非缺失的分类识别

采用上述训练得到的集成分类器，对内节/外节区域内的体素进行分类识别，分成缺失/非缺失两类。

(6)后处理与缺失体积计算

如附图5(a)所示，箭头所示为血管轮廓在OCT图像X-Z平面图中的位置，由于视网膜血管轮廓的影响，SD-OCT图像的B扫描图的相应位置上会出现上下贯穿的阴影，在IS/OS连接区内也会出现血管阴影，这些阴影往往会被误判为缺失，因此需要进行排除。首先将第1至第10层的全部三维OCT数据投影至X-Y平面，如图5(b)所示，视网膜OCT图像在X-Y平面上的投影图以及血管，图5(a)和图5(b)在图中箭头处的血管位置一一对应，本申请对投影图进行血管轮廓检测，如步骤(5)中的缺失识别结果的X坐标值与血管轮廓的X坐标值相等，则认为是缺失误检测，将相应的识别结果改为为非缺失，可有效地排除血管轮廓对内节/外节缺失检测的影响。

由于内节/外节缺失的生理特点，孤立的缺失点或者孤立的非缺失点应该不存在，因此，对缺失识别结果进一步采用形态学方法进行了孤立点的排除。

经后处理的缺失识别结果是以体素个数为单位的，最后缺失识别结果乘以体素的分辨率，转化为检测到的内节/外节缺失体积。

(7)实验方法与结果

内节/外节缺失的金标准由2个观察者在眼科医生指导下独立手动标注。对15个存在内节/外节缺失的眼外伤病人和15个正常人数据进行了本发明方法的测试。采用了留一交叉验证方法来检验本发明方法的可行性和有效性。由于在训练阶段，对非缺失样本进行了随机降采样，为了减少随机误差的影响，采用Adaboost算法训练了5个集成分类器模型，并对15个外伤数据和15个正常数据的内节/外节区域分别进行了5次缺失检测。

对外伤病人和正常人的内节/外节缺失检测体积的均值及其95％置信区间进行了比较，对采用t检验方法评估两组数据缺失体积的统计差异性。

为了进一步评估本发明方法的性能，将检测结果与第1个观察者手动标记的金标准进行比较。采用灵敏度(Sensitivity，SEN)、特异度(Specificity，SPE)和均衡正确检测率作为评估的客观指标，分别采用式(9)、式(10)和式(11)计算：

$\mathrm{SEN}=\frac{\mathrm{TP}}{\mathrm{TP}+\mathrm{FN}}\×100\%---\left(9\right)$

$\mathrm{SPE}=\frac{\mathrm{TN}}{\mathrm{TN}+\mathrm{FP}}\×100\%---\left(10\right)$

$\mathrm{BAR}=\frac{\mathrm{SEN}+\mathrm{SPE}}{2}---\left(11\right)$

其中，TP、TN、FP和FN分别表示真阳、真阴、假阳和假阴。

附图6(a)给出了由观察者1手动标记的一个外伤眼数据的内节/外节缺失三维金标准结果(附图中心的阴影部分即为内节/外节缺失区域)，附图6(b)给出了对该外伤眼数据采用本发明方法进行内节/外节缺失检测的结果，附图6(c)给出了本发明方法对一个正常眼数据进行检测的结果。附图6(a)理解为挫伤性视网膜内节/外节缺失区域期望的检测结果，附图6(b)采用本发明技术方案检测的结构，对比附图6(a)的期望的检测结果，本发明技术方案检测的结果基本与期望的检测结果一致，准确性高，附图6(c)给出的本发明方法对一个正常眼数据进行检测的结果，基本无缺失区域，是实际眼部状况一致。

由附图6的检测结果可见，本方法能有效地检测出内节/外节的缺失区域。

如图7所示，图7为外伤眼数据和正常眼数据缺失体积检测结果的比较，图中，外伤眼数据和正常眼数据中宽矩形图标表示均值，宽矩形图标中间的窄实线图标表示95％置信区间。采用本发明方法对15个正常数据进行内节/外节缺失检测的体积均值为0.0037mm³，其95％置信区间为(0.0005,0.0069)mm³；对15个外伤眼数据进行内节/外节缺失检测的体积均值为0.1626mm³，其95％置信区间为(0.0198,0.3054)mm³；外伤眼数据和正常眼数据缺失体积检测结果的比较表明，本发明的检测结果与实际眼部状况一致，正常眼的内节/外节缺失检测的体积以及95％置信区间的数据范围远远小于外伤眼的内节/外节缺失检测的体积以及95％置信区间的数据范围，本发明技术方案检测结果与眼部实际状况一致。对正常数据和外伤数据的缺失体积检测结果进行了t检验，结果表明两者存在显著的统计差异(P＝9.9112×10^-8<0.001)。

表1给出了15个外伤数据的内节/外节缺失检测结果。对于15个外伤数据来说，本发明方法的灵敏度为81.76％±10.52％，特异性为75.29％±7.12％，均衡正确检测率为78.51％±7.06％。对于15个正常数据，本发明方法的特异性为99.03％±0.73％。表明本发明方法是一种较有效的内节/外节缺失检测方法。

表1 15个外伤数据的缺失检测体积、灵敏度、特异度和均衡正确检测率

至此，一种自动检测由眼外伤引起的视网膜内节/外节缺失的方法已经实现并进行了验证。本发明融合了三维图割技术对视网膜SD-OCT进行内部分层、内节/外节区域三维纹理特征的提取、基于主成份分析方法的特征优化选择、基于Adaboost算法的集成分类器的训练与测试以及血管轮廓影响排除和孤立点排除等后处理，实现了视网膜内节/外节缺失的自动三维定量检测，对于眼外伤等临床常见眼科疾病的诊断与治疗起到了重要的辅助作用。

以上所述仅是本发明的优选实施方式，应当指出：对于本技术领域的普通技术人员来说，在不脱离本发明原理的前提下，还可以做出若干改进和润饰，这些改进和润饰也应视为本发明的保护范围。

Claims (6)

1.一种基于SD-OCT的挫伤性视网膜内外节缺失检测方法，其特征在于：包括以下步骤：

(1)图像预处理：采用多尺度三维图割方法进行视网膜内部分层，将视网膜SD-OCT图像自动分割成10个分层，产生11个表面，将包含内节/外节区域的第7层和第8层之间的数据提取出来，参考第11层进行平坦化，将平坦化后第7层和第8层之间的内节/外节区域的SD-OCT图像作为感兴趣区域，进行双线性滤波增强；

(2)体素特征提取：对步骤(1)所述平坦化后的内节/外节区域的所有体素的灰度值进行归一化处理，对内节/外节区域内的缺失体素和非缺失体素分别提取5类特征，所述5类特征包括归一化的灰度值、分块均值、分块标准差、13个方向的灰度差绝对值和基于灰度共生矩阵的特征；

(3)特征优化选择：对步骤(2)中所提取的特征采用主成份分析法进行特征选择；

(4)集成分类器的训练：将经步骤(3)选择的特征分成训练集和测试集，在训练集中，按缺失样本数量与非缺失样本数量的比例为1：1的比例从非缺失样本中随机抽取数量与缺失样本数量相同的特征样本，采用Ababoost算法进行集成分类器的训练；

(5)体素缺失/非缺失的分类识别：对测试集样本进行缺失/非缺失识别，采用步骤(4)训练得到的集成分类器对测试集样本进行分类，对每个体素的缺失与否进行分类识别；

(6)后处理与缺失体积计算：对步骤(5)分类识别得到的结果进行血管轮廓影响排除和孤立点排除，并计算视网膜内节/外节的缺失体积。

2.根据权利要求1所述的基于SD-OCT的挫伤性视网膜内外节缺失检测方法，其特征在于：所述步骤(1)具体包括以下步骤：

(1-1)视网膜内部分层：采用多尺度三维图搜索方法将视网膜SD-OCT图像自动分割成10个分层，产生11个表面；

(1-2)内节/外节区域提取和平坦化：将包含内节/外节区域的第7层和第8层之间的数据提取出来，参考第11层进行平坦化，将平坦化后的第7层与第8层之间的SD-OCT图像作为感兴趣区域，进行内节/外节的缺失检测。

3.根据权利要求1所述的基于SD-OCT的挫伤性视网膜内外节缺失检测方法，其特征在于：所述步骤(2)所述5类特征的提取具体包括以下步骤：

(2-1)归一化灰度值表示每个体素的灰度级，对所述感兴趣区域内的所有体素的灰度值采用如式(1)所示的方法进行线性归一化至0～255；

$I_{\mathrm{normalized}}(i,j,k)=\frac{I_{\mathrm{original}}(i,j,k)-I_{\min }}{I_{\max }-I_{\min }}\&times;255---\left(1\right)$

其中，(i,j,k)表示体素的坐标，I_normalized(i,j,k)表示体素(i,j,k)归一化后的灰度值，I_original(i,j,k)表示体素(i,j,k)原始的灰度值，I_min和I_max分别表示感兴趣区域内体素灰度值的最小值和最大值；所述内节/外节缺失区域的灰度值低于非缺失区域的灰度值；

(2-2)分块均值用于表示以体素(i,j,k)为中心的分块的灰度均值，所述分块为5×5×5分块，所述5×5×5分块的灰度均值M_block(i,j,k)，计算公式如式(2)：

$M_{\mathrm{block}}(i,j,k)=\frac{1}{125}\underset{l=i-2}{\overset{i+2}{\mathrm{\&Sigma;}}}\underset{m=j-2}{\overset{j+2}{\mathrm{\&Sigma;}}}\underset{n=k-2}{\overset{k+2}{\mathrm{\&Sigma;}}}{I}_{\mathrm{normalized}}(l,m,n)---\left(2\right)$

(2-3)分块标准差STD_block(i,j,k)用于表示以体素(i,j,k)为中心的5×5×5分块的灰度值变化情况，表征灰度值在5×5×5分块中的分散程度；根据分块标准差STD_block(i,j,k)识别分块体素全部位于缺失区域或者分块体素全部位于非缺失区域或者分块中同时包含缺失区域和非缺失区域的体素，计算公式如下：

${\mathrm{STD}}_{\mathrm{block}}(i,j,k)=\sqrt{\frac{\underset{l=i-2}{\overset{i+2}{\mathrm{\&Sigma;}}}\underset{m=j-2}{\overset{j+2}{\mathrm{\&Sigma;}}}\underset{n=k-2}{\overset{k+2}{\mathrm{\&Sigma;}}}{(I_{\mathrm{normalized}}(l,m,n)-M_{\mathrm{block}}\left(\mathrm{imjmk}\right))}^2}{124}}---\left(3\right)$

(2-4)13个方向上的灰度差绝对值之和定义为：中心点体素(i,j,k)的灰度值与13个方向周围体素的灰度差绝对值之和，反映中心点体素在缺失区域或非缺失区域的位置，依据所述13个方向上的灰度差绝对值之和AID_step(i,j,k)判断中心点体素(i,j,k)位于缺失与非缺失交界处或者缺失区域或者非缺失区域，选择两种距离步长计算13个方向上的灰度差绝对值之和AID_step(i,j,k)，两种距离步长为step＝1，2，具体计算公式如公式(4)：

${\mathrm{AID}}_{\mathrm{step}}(i,j,k)=\underset{13\mathrm{directions}}{\mathrm{\&Sigma;}}|I_{\mathrm{normalized}}(i,j,k)-I_{\mathrm{normalized}}(l,m,n)|---\left(4\right)$

表示13个方向上的灰度差绝对值之和；

所述13个方向的选取方法为：假设α₁表示X轴与X-Y平面投影图的夹角，α₂表示Z轴与X-Y平面投影图的夹角，那么所述13个方向分别为(α₁,α₂)＝(0,90°)、(45°,90°)、(90°,90°)、(135°,90°)、(0,45°)、(180°,45°)、(90°,45°)、(-90°,45°)、(0,0)、(45°,45°)、(135°,45°)、(-45°,45°)、(-135°,45°)；

(2-5)基于灰度共生矩阵特征描述若干立体切面之间灰度级的空间相关性，所述灰度共生矩阵为三维体数据的灰度共生矩阵，采用三维方法在13个方向上跨切面搜索不同灰度级，构造13个灰度共生矩阵，13个灰度共生矩阵构造方法具体包括以下步骤：

对每个5×5×5分块构造所述13个不同方向上的灰度共生矩阵；

分别提取4个基于所述灰度共生矩阵的特征，所述4个基于灰度共生矩阵的特征为：(a)对比度、(b)相关性、(c)能量和(d)均匀性；

对比度Contrast_direction m(i,j,k)、相关性Correlation_direction m(i,j,k)、能量Energy_direction m(i,j,k)和均匀性Homogeneity_direction m(i,j,k)分别采用式(5)、式(6)、式(7)和式(8)进行计算：

${\mathrm{Contrast}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N_g}{\mathrm{\&Sigma;}}}{(x-y)}^{2}p(x,y),m=\mathrm{1,2},...,13---\left(5\right)$

${\mathrm{Correlation}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N_g}{\mathrm{\&Sigma;}}}\frac{(x-{\mathrm{\&mu;}}_x)(y-{\mathrm{\&mu;}}_y)p(x,y)}{{\mathrm{\&sigma;}}_x{\mathrm{\&sigma;}}_y},m=\mathrm{1,2},...,13---\left(6\right)$

${\mathrm{Energy}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N_g}{\mathrm{\&Sigma;}}}{p}^2(x,y),m=\mathrm{1,2},...,13---\left(7\right)$

${\mathrm{Homogeneity}}_{\mathrm{direction\; m}}(i,j,k)=\underset{x=1}{\overset{N_g}{\mathrm{\&Sigma;}}}\underset{y=1}{\overset{N_g}{\mathrm{\&Sigma;}}}\frac{p(x,y)}{1+|x-y|},m=\mathrm{1,2},...,13---\left(8\right)$

其中，p(x,y)表示灰度共生矩阵中第(x,y)个元素；N_g表示量化后的灰度级，取N_g＝8；μ_x和σ_x分别表示第x行的均值和标准差；μ_y和σ_y分别表示第y列的均值和标准差。

4.根据权利要求1所述的基于SD-OCT的挫伤性视网膜内外节缺失检测方法，其特征在于，所述步骤(3)具体包括以下步骤：对于步骤(2)所述提取的5类特征采用主成份分析方法进行特征选择，将经过特征选择后的结果作为新的特征，根据新特征所代表的数据总方差的比例大于或者等于90％的原则，选择前12个主成分作为新特征；

完成特征优化选择后，感兴趣区域内的每个体素用一个12维的特征向量表征，在特征空间中进行分类。

5.根据权利要求1所述的基于SD-OCT的挫伤性视网膜内外节缺失检测方法，其特征在于，所述步骤(4)具体包括以下步骤：将经步骤(3)选择的特征分成训练集和测试集，在训练集中，通过下采样方式以缺失样本数量与非缺失样本数量的比例为1：1构造新的训练集，采用分类器集成的方法进行分类，以决策树为弱分类器，采用Adaboost算法进行集成分类器的训练。

6.根据权利要求1所述的基于SD-OCT的挫伤性视网膜内外节缺失检测方法，其特征在于，所述步骤(6)具体包括以下步骤：

(6-1)血管轮廓影响的排除：将第1至第10层的全部三维SD-OCT数据投影至X-Y平面，对投影图进行血管轮廓检测，将步骤(5)的缺失识别结果的X坐标值等于血管轮廓的X坐标值，则该X坐标值处的体素被认为是缺失误识别，将对应的识别结果修改为非缺失区域；

(6-2)将步骤(6-1)的缺失识别结果采用形态学方法进行孤立点的排除；

(6-3)将步骤(6-2)的孤立点的排除后的缺失识别结果乘以体素的分辨率，转化为检测到的内节/外节缺失体积。