一种眼底OCT图像中视网膜分层方法与流程

本发明属于医学图像处理领域，涉及oct图像的处理方法，尤其是涉及了一种眼底oct图像中视网膜分层方法。
背景技术：
光学相干层析成像(oct)是一种高分辨率横截面成像技术，通过收集样品表面的微结构背反射光的光子，并与参考光路光形成干涉，实现成像。oct设备最早应用于生物组织的活体分析，主要为眼科领域；随着研究深入，不同厂家生产的oct及其图像给眼科疾病的临床诊断带来巨大的参考作用。视网膜层总厚度和视网膜内层厚度是临床眼科医生的重要指标，如视网膜神经纤维层(rnfl)、视网膜色素上皮(rpe)层厚度是青光眼和糖尿病视网膜病变的重要诊断依据，尤其对糖尿病黄斑水肿等高致盲眼科疾病的诊断实现了早起筛查和病情冬天观测。在实际应用中，如果医生使用自己的经验，在做诊时手动对视网膜进行分析，既费时又费力。近些年，各大oct设备厂商都开始提出自动视网膜分层算法，并对分层效率及分层精度极大关注。目前自动分层算法已经得到广泛的研究和商业应用。但是，oct图像的成像质量极易受散斑噪声、设备对焦的影响，都会降低图像的低图像对比度以及图像中的血管和其他结构都给算法实现自动分层制造层层障碍。正是由于上述问题的存在，分层算法仍然是研究的关注点；视网膜分层方法首先是基于灰度变化，并以此为基础，配合各种滤波和强度峰搜索，这类方法通常是不稳定的；基于灰度变化的方法受制于图像的低对比度和复杂信息如血管等的影响，且所能分出的层数有限，难以分辨细胞层(gcl)和内网层(ipl)的分界线以及rpe层和感光层的分界线。随着研究的深入，结合灰度、梯度信息在应用过程中的限制，基于活动轮廓模型的方法和基于图论的方法是目前分层算法的主流；但是，活动轮廓方法的轮廓特征通常不同精确到像素量级。图论方法由于构造了能量函数和限制条件，并计算了梯度和灰度信息建立，与活动轮廓方法相比，图论方法可保证寻找到全局最小值。然而，图论算法对噪声和图像退化非常敏感。技术实现要素：针对于
背景技术：
中存在的问题，本发明的目的在于提供了一种眼底oct图像中视网膜分层方法，是基于多种特征的光学相干层析成像视网膜层状结构分割处理方法。本发明结合了各种方法的优势，分多步对视网膜图像进行分层，综合考虑了图像的灰度、梯度、a扫描信号轮廓、梯度路径等多种特征，从而达到实现视网膜快速分层的目的。该方法具有较好的重复性与自适应性，分割的精确度也达到领先水平，结果令人满意，符合临床检查中应用要求，为国产oct装置在眼科应用的性能提升奠定技术基础。如图1所示，本发明采用的技术方案是包括以下步骤：1)采集眼底视网膜的oct图像作为原始图像；2)使用权系数矩阵模板遍历整个原始图像进行模板滤波；3)采用以下方式进行处理获得视网膜色素上皮层(rpe)和脉络膜(choroid)的分界线、内界膜(ilm)层、感光细胞层内节(is)和感光细胞层内节(os)之间的分界线、神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线、外丛状层(opl)和外核层(ονl)之间的分界线、内核层(inl)与外丛状层(opl)之间的分界线以及内丛状层(ipl)和内核层(inl)之间的分界线，完成眼底oct图像中的视网膜分层。所述步骤3)具体为：3.1)针对oct图像的每一列记录灰度最大值所在像素点，将各列的灰度值最大值所在像素点连接一起，作为视网膜色素上皮层(rpe)和脉络膜(choroid)之间的分界线；3.2)首先使用模板[-1,-1,1,1]对原始图像进行滤波；接着以原始图像从上到下的前10行像素作为背景，计算图像背景处的噪声方差σ，以噪声方差σ的n倍作为阈值，滤波后的原始图像进行二值化处理；然后使用sobel边缘检测算子检测二值化处理后的图像中的边界，使用样条曲线拟合边界作为内界膜(ilm)层(通常认为ilm没有厚度，为一像素点宽度)；3.3)提取原始图像的每列，计算每列中各个像素点的二阶梯度，以视网膜色素上皮层(rpe)和脉络膜(choroid)的分界线为基准向上搜索每一列第一个二阶梯度为零的像素点；然后，将搜索到的所有像素点连接为一条曲线，搜索曲线上所有间断点并剔除，使用样条曲线对剔除间断点后的曲线进行拟合；使用设置窗口宽度，使用三次多项式savitzky-golay多项式对拟合后的曲线进行平滑，平滑后的曲线作为感光细胞层内节(is)和感光细胞层内节(os)之间的分界线；所述步骤3.3)中的间断点为出现的可去间断点、跳跃间断点、振荡间断点和无穷间断点。3.4)建立能量函数，根据内界膜(ilm)层处理获得神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线；3.5)构建内外能量函数，并根据感光细胞层内节(is)和感光细胞层内节(os)之间的分界线及神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线处理获得外丛状层(opl)和外核层(ονl)之间的分界线；3.6)构建权重函数，并根据神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线及外丛状层(opl)和外核层(onl)之间的分界线处理获得内核层(inl)与外丛状层(opl)之间的分界线以及内丛状层(ipl)和内核层(inl)之间的分界线。所述步骤3.2)中，具体实施中，调整所述n值使得获得的内界膜(ilm)层在视网膜色素上皮层(rpe)和脉络膜(choroid)的分界线之上。所述步骤3.4)具体为：3.4.1)提取原始图像的每列，计算每列中各个像素点的纵向灰度梯度g(r,y)，其中y为列中的纵坐标，r为列编号；3.4.2)从内界膜(ilm)层开始向下搜索，搜索每列中纵向灰度梯度的第一个极小点坐标，记为l(r)；3.4.3)同时建立以下能量函数：e(r)＝g(r,y)+|l(r-1)-l(r)|/500+|l(r+1)-l(r)|/250针对各列，以mine(r)为目标寻优，得到能量函数最小时的像素点，组成像素点集合；3.4.4)对像素点集合里的所有像素点进行多项式拟合，拟合得到神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线。所述步骤3.5)具体为：3.5.1)以感光细胞层内节(is)和感光细胞层内节(os)之间的分界线及神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线为第一指定合围区域的上下边界，对原始图像的第一指定合围区域的每列像素点进行高斯梯度滤波，降低散弹噪声强度；3.5.2)然后构建以下公式的内外能量函数：eint(i)＝α[d-|ki-ki-1|]2+β|ki+1-ki+ki-1|2eext(i)＝γ|g’(xi,yi)×ι(xi,yi)|2其中，eint(i)表示内部能量强度，eext(i)表示外部能量强度，g’为原始图像的二阶梯度，i(xi,yi)为像素点(xi,yi)的灰度，ki表示列中的第i个像素点的坐标，xi为像素点ki在图像中的横坐标，yi为像素点ki在图像中的纵坐标；α、β、γ分别表示第一、第二、第三自由参数，通过优化求解，d为相邻像素点之间的平均距离；3.5.3)针对每一列，从列像素集合随机选择n个像素点，求解获得使得总能量eint(i)+eext(i)最小的像素点，将该像素点的总能量作为能量最小值，将能量最小值的像素点从列像素集合排除；初始情况下，列像素集合为列中的所有像素点构成。3.5.4)不断重复上述步骤3.5.3)，从中再次随机选择n个像素点，直至所有像素点的总能量eint(i)+eext(i)之和保持不变；3.5.5)将各列搜索到能量最小值的像素点进行样条曲线拟合，获得外丛状层(opl)和外核层(ονl)之间的分界线。所述步骤3.6)具体为：3.6.1)以神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线及外丛状层(opl)和外核层(onl)之间的分界线为第二指定合围区域的上下边界，对原始图像的第二指定合围区域中各个像素点的纵向梯度；3.6.2)构建以下权重函数，计算相邻两个像素点之间的权重：wab＝2.0001-g(a)-g(b)其中，wab表示第二指定合围区域中像素点a和像素点b之间的权重，像素点a和像素点b相邻，g(a)和g(b)分别为第二指定合围区域中的像素点a和像素点b的纵向梯度绝对值；3.6.3)第二指定合围区域中，从最左列的一像素点到最右列的一像素点的所有路径中搜索找到最短路径，最短路径为长度最短的路径，路径的长度为路径上所有相邻像素点的权重之和，搜索时设定最左列和最右列的像素点梯度均为1；3.6.4)使用弗洛伊德算法(floyd算法)输入路径的长度和权重进行第一次搜索获得最短路径，对最短路径进行移动平均平滑处理，作为内核层(inl)与外丛状层(opl)之间的分界线；3.6.4)将内核层(inl)与外丛状层(opl)之间的分界线从第二指定合围区域去除，再重复使用弗洛伊德算法(floyd算法)输入路径的长度和权重进行第二次搜索获得最短路径，对最短路径进行移动平均平滑处理，作为内丛状层(ipl)和内核层(inl)之间的分界线。本发明采集眼底的oct图像后，完成模版滤波，rpe/choroid灰度分层，ilm层识别，is/os梯度搜索，nfc/gcl特征提取，opl/onl能量函数优化，inl/opl、ipl/inl路径搜索，最终完成oct图像中视网膜的分割。本发明具有的有益效果是：本发明的视网膜分层方法，挖掘多个特征，能够利用普通配置的计算机实现多个层结构的自动分割，对视网膜复杂病症层状结构实现有效的自动探测。由于oct图像中存在散点噪声，本发明只在灰度分辨明显的rpe层使用灰度搜索算法。本发明提出的自动分层方法，避开了预处理手段，分多步对视网膜图像进行分层，在每个a扫描上搜索了曲线特征，然后对各层的初步分层结果进行连续性和完整性判断，修正不符合要求的分割点。在内部细节上，利用图论和基于动态规划的最短路算法，自动寻找出准确的视网膜二层结构，验证了算法的准确性和可靠性。各层边界的结果具有较强的可比性，适合多种临床诊断应用。更重要的是，本发明方法能再在视频速率上部署应用，通过现代计算机gpu的并行计算能力，本发明的插入计算量较少，信号处理和图像渲染所需资源也较少。本发明采用进行了系列化的分层，对图像原有处理步骤少，检测效率具有一定的优势。附图说明图1是本发明方法的流程图。图2是实施例采集的眼底图。图3是图2视网膜样本的分层结果图。具体实施方式以下结合附图及实施例，对本发明进行进一步详细说明。应当理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限定本发明。本发明的实施例如下：1)采集眼底视网膜的oct图像作为原始图像，如图2所示。2)使用权系数矩阵模板[1/9，1/9，1/9；1/6，1/6，1/6；1/9，1/9，1/9；]遍历整个原始图像进行模板滤波；3)采用以下方式进行处理获得各条分界线，进行眼底oct图像中的视网膜分层。3.1)针对oct图像的每一列记录灰度最大值所在像素点，将各列的灰度值最大值所在像素点连接一起，作为视网膜色素上皮层(rpe)和脉络膜(choroid)之间的分界线；3.2)首先使用模板[-1,-1,1,1]对原始图像进行滤波；接着以原始图像从上到下的前10行像素作为背景，计算图像背景处的噪声方差σ，以噪声方差σ的7倍作为阈值，滤波后的原始图像进行二值化处理；然后使用sobel边缘检测算子检测二值化处理后的图像中的边界，使用样条曲线拟合边界作为内界膜(ilm)层；在本实施例中，n值为7，满足上述条件：使得获得的内界膜(ilm)层在视网膜色素上皮层(rpe)和脉络膜(choroid)的分界线之上。3.3)提取原始图像的每列，计算每列中各个像素点的二阶梯度，以视网膜色素上皮层(rpe)和脉络膜(choroid)的分界线为基准向上搜索每一列第一个二阶梯度为零的像素点；然后，将搜索到的所有像素点连接为一条曲线，搜索曲线上所有间断点并剔除，使用样条曲线对剔除间断点后的曲线进行拟合；使用三次多项式savitzky-golay多项式对拟合后的曲线进行平滑，平滑后的曲线作为感光细胞层内节(is)和感光细胞层内节(os)之间的分界线；3.4)建立能量函数，根据内界膜(ilm)层处理获得神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线；3.4.1)提取原始图像的每列，计算每列中各个像素点的纵向灰度梯度g(r,y)，其中y为列中的纵坐标，r为列编号；3.4.2)从内界膜(ilm)层开始向下搜索，搜索每列中纵向灰度梯度的第一个极小点坐标，记为l(r)；3.4.3)同时建立以下能量函数：e(r)＝g(r,y)+|l(r-1)-l(r)|/500+|l(r+1)-l(r)|/250针对各列，以mine(r)为目标寻优，得到能量函数最小时的像素点，组成像素点集合；3.4.4)对像素点集合里的所有像素点进行多项式拟合，拟合得到神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线。3.5)构建内外能量函数，并根据感光细胞层内节(is)和感光细胞层内节(os)之间的分界线及神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线处理获得外丛状层(opl)和外核层(ονl)之间的分界线；3.5.1)以感光细胞层内节(is)和感光细胞层内节(os)之间的分界线及神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线为第一指定合围区域的边界，对原始图像的第一指定合围区域的每列像素点进行高斯梯度滤波，进一步降低散弹噪声强度；3.5.2)然后构建以下公式的内外能量函数：eint(i)＝α[d-|ki-ki-1|]2+β|ki+1-ki+ki-1|2eext(i)＝γ|g’(xi,yi)×ι(xi,yi)|23.5.3)针对每一列，从列像素集合随机选择n个像素点，求解获得使得总能量eint(i)+eext(i)最小的像素点，将该像素点的总能量作为能量最小值，将能量最小值的像素点从列像素集合排除；3.5.4)不断重复上述步骤3.5.3)，从中再次随机选择n个像素点，直至所有像素点的总能量eint(i)+eext(i)之和保持不变；3.5.5)将各列搜索到能量最小值的像素点进行样条曲线拟合，获得外丛状层(opl)和外核层(ονl)之间的分界线。3.6)构建权重函数，并根据神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线及外丛状层(opl)和外核层(onl)之间的分界线处理获得内核层(inl)与外丛状层(opl)之间的分界线以及内丛状层(ipl)和内核层(inl)之间的分界线。3.6.1)以神经纤维层(nfl)和神经节细胞层(gcl)之间的分界线及外丛状层(opl)和外核层(onl)之间的分界线为第二指定合围区域的边界，对原始图像的第二指定合围区域中各个像素点的纵向梯度；3.6.2)构建以下权重函数，计算相邻两个像素点之间的权重：wab＝2.0001-g(a)-g(b)其中，wab表示第二指定合围区域中像素点a和像素点b之间的权重，像素点a和像素点b相邻，g(a)和g(b)分别为第二指定合围区域中的像素点a和像素点b的纵向梯度绝对值；3.6.3)第二指定合围区域中，从最左列的一像素点到最右列的一像素点的所有路径中搜索找到最短路径，最短路径为长度最短的路径，路径的长度为路径上所有相邻像素点的权重之和，搜索时设定最左列和最右列的像素点梯度均为1；3.6.4)使用弗洛伊德算法输入路径的长度和权重进行第一次搜索获得最短路径，对最短路径进行移动平均平滑处理，作为内核层(inl)与外丛状层(opl)之间的分界线；3.6.4)将内核层(inl)与外丛状层(opl)之间的分界线从第二指定合围区域去除，再重复使用弗洛伊德算法输入路径的长度和权重进行第二次搜索获得最短路径，对最短路径进行移动平均平滑处理，作为内丛状层(ipl)和内核层(inl)之间的分界线。图2视网膜样本的分层结果如图3所示。如图3所示，获得各条分界线从上到下依次为：内界膜(ilm)层、神经纤维层(nfl)和神经节细胞层(gcl)、内丛状层(ipl)、内核层(inl)、外丛状层(opl)、外核层(ονl)、感光细胞层内节(is)和感光细胞层内节(os)、视网膜色素上皮层(rpe)和脉络膜(choroid)。在本实施例中，按照医学标准获得的眼底视网膜分层和本发明算法实施例获得的结果进行比较，两者差异边界/分界线取平局偏差和标准差见表1，以像素为单位。表1视网膜分层结果层平均误差/像素标准差/像素ilm0.610.18nfl/gcl1.120.04ipl/inl1.360.42inl/opl1.540.21opl/onl1.270.23is/os1.540.20rpe/脉络膜0.320.13实验结果表明，本发明对于眼科oct图像的均有较好分割的效果，能实现相对准确的眼底视网膜分层。对比现有报道的单一图论、灰度等方案，显著提高了分割精度和稳定性，显示了本发明方法的优势。综合来说，本发明能够利用普通配置的计算机实现多个层结构的自动分割，对视网膜复杂病症层状结构实现有效的自动探测，采用进行了系列化的分层，对图像原有处理步骤少，检测效率具有一定的优势。在本发明实施例中，本领域普通技术人员还可以理解，实现上述实施例方法中的全部或部分步骤是可以通过程序来指令相关的硬件来完成，所述的程序可以在存储于一计算机可读取存储介质中，所述的存储介质，包括rom/ram、磁盘、光盘等。以上所述仅为本发明的较佳实施例而已，并不用以限制本发明，凡在本发明的精神和原则之内所作的任何修改、等同替换和改进等，均应包含在本发明的保护范围之内。当前第1页12