CN114332087B - Octa图像的三维皮层表面分割方法及系统 - Google Patents

Abstract

本发明涉及OCTA图像的三维皮层表面分割方法及系统，包括：对OCTA图像双边滤波得到滤波后图像；去除滤波后图像中的盖玻片得到目标三维图像；以目标三维图像正中间的二维纵向剖面作为起始平面，分割皮层边界得到起始平面的边界分割线；以起始平面为基准，往前后两侧依次选取预设厚度的三维图像数据块，直至将整个目标三维图像划分完毕；以靠近起始平面的三维图像数据块开始，结合起始平面的边界分割线，依次从内向外计算各三维图像数据块的皮层边界并组成稀疏的皮层表面；对稀疏的皮层表面进行平滑处理得到稀疏的目标皮层表面；对稀疏的目标皮层表面扩充得到密集的三维皮层表面。本发明无需手动分割皮层表面，且三维皮层表面分割精度高。

Description

OCTA图像的三维皮层表面分割方法及系统

技术领域

本发明属于医疗影像处理技术领域，具体涉及OCTA图像的三维皮层表面分割方法及系统。

背景技术

OCT(optical coherence tomography)成像技术主要利用干涉法探测样品背向散射的相干光信号，再通过后处理获取样品深度方向的信息并加以结构成像，具有非侵入、无辐射损伤、高灵敏度、高分辨率等优势。历经二十多年的创新发展，OCT成像技术日趋成熟，针对不同的应用场景也有相应的特异性改进。其中，OCT血管造影成像(OCT-basedangiography，OCTA)是OCT技术的一种功能性扩展，可对一定深度范围内(1-2mm)的微血管网络结构进行高分辨率的可视化成像，同时无需任何外源性的荧光造影剂辅助。OCTA在血管疾病领域有着其独特的优势，可以直接观察病变血管及其旁支随时间的动态变化情况，同时能够获得受损血管深处毛细血管和组织信息。目前，OCTA成像技术已被用于眼睛、皮肤、大脑及肠道等多个部位器官的血管成像临床实践和科学研究。

由于相关OCTA成像技术起步不久，目前针对该类图像的分析处理方法较少，相关人员仍需要通过手动处理。其中，分割去除血管皮层外的非目标区域是处理过程中的关键步骤之一。OCTA技术主要是生成三维图像，而现有条件下除了人为手动剔除外，可使用的边界分割方法还是用于二维剖面图像，并不能用于完成三维皮层边界分割。再者，血管阻塞病变后会先后出现膨胀和塌陷的现象，在图像中会有反映，常规方法在遇到病变血管区域将无法有效对皮层边界进行分割。另外，通过计算血管皮层边界和盖玻片间的体积变化情况，可量化反应血管膨胀或塌陷程度，以供临床或科研人员分析判断疾病进展情况，人为剔除非目标区域将导致相关特征参数的缺失。

发明内容

基于现有技术中存在的上述缺点和不足，本发明的目的之一是至少解决现有技术中存在的上述问题之一或多个，换言之，本发明的目的之一是提供满足前述需求之一或多个的OCTA图像的三维皮层表面分割方法及系统。

为了达到上述发明目的，本发明采用以下技术方案：

OCTA图像的三维皮层表面分割方法，包括以下步骤：

S1、对三维OCTA图像进行双边滤波，得到滤波后图像；

S2、去除滤波后图像中的盖玻片，得到目标三维图像；

S3、以目标三维图像正中间的二维纵向剖面作为起始平面，分割皮层边界，得到起始平面的边界分割线；

S4、以起始平面为基准，往前后两侧依次选取预设厚度的三维图像数据块，直至将整个目标三维图像划分完毕；其中，相邻三维图像数据块有一定重叠；

S5、以靠近起始平面的三维图像数据块开始，结合起始平面的边界分割线，依次从内向外计算各三维图像数据块的皮层边界，并组成稀疏的皮层表面；

S6、对稀疏的皮层表面进行平滑处理，得到稀疏的目标皮层表面；

S7、对稀疏的目标皮层表面扩充得到密集的三维皮层表面。

作为优选方案，所述步骤S2，包括以下步骤：

S21、利用最大类间方差法将滤波后图像转换为0、1表示的掩膜；

S22、依次使用开运算、关运算和腐蚀运算对掩膜进行处理，得到目标掩膜；

S23、间隔选取目标掩膜的数个二维纵向剖面，以横坐标从左往右的顺序，依次从上往下开始判断每列对应像素点的掩膜数据；若为1，则提取相应目标像素点的坐标位置信息，开始下一列判断；若为0，则继续往下判断直至判断的像素点数量超过目标数量阈值，开始下一列判断；

S24、根据所有目标像素点的坐标位置信息拟合二阶曲面，曲面上方的所有杂质数据均剔除，得到剔除上层边缘面的图像；

S25、对剔除上层边缘面的图像执行步骤S21-S24以剔除下层边缘面，得到目标三维图像。

作为优选方案，所述步骤S3中，分割皮层边界的过程包括以下步骤：

S31、将起始平面的二维灰度图像转换为路径成本图；

S32、基于路径成本图和动态规划算法得到成本综合最小的从左往右的连续路径，作为边界分割线。

作为优选方案，所述步骤31包括以下步骤：

S311、将起始平面的二维灰度图像与核为3的y轴方向Sobel算子卷积，得到图像A；

将起始平面的二维灰度图像经过核为5的Canny算子提取边缘，再与长为3、宽为1的矩形结构化元素进行开运算，得到图像B；

S312、将图像A和图像B转换为路径成本图C：

C＝1-[ω*normalize(A)+(1-ω)*normalize(B)]

其中，ω为图像A的权重，1-ω为图像B的权重，normalize(A)为对图像A进行归一化至[0，1]区间，normalize(B)为对图像B进行归一化至[0，1]区间。

作为优选方案，所述步骤S32包括：

设路径成本图C的长为M，高为N，各像素点(i，j)的成本记为C(i，j)，i＝1,2,3,…,M，j＝1,2,3,…,N；

以左端的N个像素点分别作为候选的N个起始点，从左到右依次寻找路径的下一个像素点，直至找到路径右端的像素点；其中，路径中相邻两像素点的成本分别记为C(i'，j')、C(i'+1，j”)，需满足以下条件：min C(i'+1，j”)，j”∈[j'-1，j'+1]；i'依次取值为1,2,3,…,M-1；

分别统计各路径所有像素点的成本总和，筛选得到成本综合最小的目标路径；

采用Savitzky-Golay卷积平滑算子对目标路径进行平滑处理，得到边界分割线。

作为优选方案，所述步骤S5，包括以下步骤：

S51、将预设厚度为H的三维图像数据块上每一层二维图像上对应的像素点的像素值取平均，得到目标二维图像：

其中，I(k，l，h)为三维图像数据块上各层二维图像的像素点的像素值，二维图像的层数为H；

S52、对目标二维图像依次执行所述步骤S31和S32，得到候选边界线L1；

S53、基于与起始平面所处的三维图像数据块的相邻三维图像数据块，根据起始平面的边界分割线起始点以及相邻三维图像数据块对应的目标二维图像转换的路径成本图，得到成本综合最小的从左往右的连续路径作为候选边界线L2；以此类推，得到各个三维图像数据块的候选边界线L2；

S54、将同一三维图像数据块对应的候选边界线L1与候选边界线L2融合得到三维图像数据块的边界分割线；

S55、将所有三维图像数据块的边界分割线汇合，得到稀疏的皮层表面。

作为优选方案，所述步骤S54，包括：

判断候选边界线L1与候选边界线L2同一横坐标的两个点(x，y₁)、(x，y₂)的纵坐标的差值是否处于预设范围内；

若是，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为(y₁+y₂)/2；

若否，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为y：

其中，y₀为相邻三维图像数据块的边界分割线的同一横坐标的点的纵坐标。

作为优选方案，所述步骤S6中，对稀疏的皮层表面依次进行中值滤波。

作为优选方案，所述步骤S7，包括：

将三维图像数据块的边界分割线作为其中间层的边界分割线；

基于相邻两个三维图像数据块的中间层的边界分割线，两中间层之间其他层的边界分割线通过两中间层的边界分割线等差递增或等差递减获得；

将所有三维图像数据块的所有层的边界分割线经过平滑得到密集的三维皮层表面。

本发明还提供OCTA图像的三维皮层表面分割系统，应用如上任一方案所述的三维皮层表面分割方法，所述三维皮层表面分割系统包括：

双边滤波模块，用于对三维OCTA图像进行双边滤波，得到滤波后图像；

盖玻片去除模块，用于去除滤波后图像中的盖玻片，得到目标三维图像；

皮层边界分割模块，用于以目标三维图像正中间的二维纵向剖面作为起始平面，分割皮层边界，得到起始平面的边界分割线；

图像划分模块，用于以起始平面为基准，往前后两侧依次选取预设厚度的三维图像数据块，直至将整个目标三维图像划分完毕；其中，相邻三维图像数据块有一定重叠；

皮层表面构建模块，用于以靠近起始平面的三维图像数据块开始，结合起始平面的边界分割线，依次从内向外计算各三维图像数据块的皮层边界，并组成稀疏的皮层表面；

平滑模块，用于对稀疏的皮层表面进行平滑处理，得到稀疏的目标皮层表面；

扩充模块，用于对稀疏的目标皮层表面扩充得到密集的三维皮层表面。

本发明与现有技术相比，有益效果是：

本发明的OCTA图像的三维皮层表面分割方法及系统，无需手动分割皮层表面，且三维皮层表面分割精度高。

附图说明

图1是本发明实施例1的OCTA图像的三维皮层表面分割方法的流程图；

图2是本发明实施例1的OCTA图像的三维皮层表面分割方法的分割效果图；

图3是本发明实施例1的OCTA图像的三维皮层表面分割系统的模块构架图。

具体实施方式

为了更清楚地说明本发明实施例，下面将对照附图说明本发明的具体实施方式。显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他的附图，并获得其他的实施方式。

实施例1：

如图1所示，本实施例的OCTA图像的三维皮层表面分割方法，包括以下步骤：

S1、对三维OCTA图像进行双边滤波，得到滤波后图像；

具体地，利用双边滤波器对三维OCTA图像进行滤波，亮度差参数选定为5，空间距离参数选定为75，以突出图像中的边界信息。

S2、去除滤波后图像中的盖玻片，得到目标三维图像；

具体地，上述步骤S2包括以下步骤：

S21、利用最大类间方差法将滤波后图像转换为0、1表示的掩膜；

S22、依次使用开运算、关运算和腐蚀运算对掩膜进行处理，得到目标掩膜；

S23、间隔选取目标掩膜的数个二维纵向剖面，以横坐标从左往右的顺序，依次从上往下开始判断每列对应像素点的掩膜数据；若为1，则提取相应目标像素点的坐标位置信息，开始下一列判断；若为0，则继续往下判断直至判断的像素点数量超过目标数量阈值，开始下一列判断；

S24、根据所有目标像素点的坐标位置信息拟合二阶曲面，曲面上方的所有杂质数据均剔除，得到剔除上层边缘面的图像；由于盖玻片为规则物体，因此拟合二阶曲面即可剔除杂质数据；

S25、对剔除上层边缘面的图像执行步骤S21-S24以剔除下层边缘面，得到目标三维图像。由于盖玻片具有一定厚度，所以剔除盖玻片的上层边缘面后，按照上述相同的步骤即可去除盖玻片的下层边缘面。

S3、以目标三维图像正中间的二维纵向剖面作为起始平面，分割皮层边界，得到起始平面的边界分割线；

具体地，上述步骤S3中，分割皮层边界的过程包括以下步骤：

S31、将起始平面的二维灰度图像转换为路径成本图；

S32、基于路径成本图和动态规划算法得到成本综合最小的从左往右的连续路径，作为边界分割线。

其中，上述步骤31包括以下步骤：

S311、将起始平面的二维灰度图像与核为3的y轴方向Sobel算子卷积，得到图像A；

将起始平面的二维灰度图像经过核为5的Canny算子提取边缘，再与长为3、宽为1的矩形结构化元素进行开运算，得到图像B；

S312、将图像A和图像B转换为路径成本图C：

C＝1-[ω*normalize(A)+(1-ω)*normalize(B)]

其中，ω为图像A的权重，1-ω为图像B的权重，normalize(A)为对图像A进行归一化至[0，1]区间，normalize(B)为对图像B进行归一化至[0，1]区间。ω优选取值为0.5。

上述步骤S32具体包括：

设路径成本图C的长为M，高为N，各像素点(i，j)的成本记为C(i，j)，i＝1,2,3,…,M，j＝1,2,3,…,N；

以左端的N个像素点分别作为候选的N个起始点，从左到右依次寻找路径的下一个像素点，直至找到路径右端的像素点；其中，路径中相邻两像素点的成本分别记为C(i'，j')、C(i'+1，j”)，需满足以下条件：min C(i'+1，j”)，j”∈[j'-1，j'+1]；i'依次取值为1,2,3,…,M-1；

由于存在候选的N个起始点，相应的路径也有N条；

分别统计各路径所有像素点的成本总和，筛选得到成本综合最小的目标路径；

采用Savitzky-Golay卷积平滑算子对目标路径进行平滑处理，得到边界分割线。

其中，得到目标路径之后，还可以从右往左回溯目标路径(回溯过程可参考上述相邻像素点的寻找过程)，直至回溯完成，得到目标路径；如此设计，在前统计的N条路径无需存储相应的像素点的位置信息，节约内存。

S4、基于目标三维图像，以起始平面为基准，往前后两侧依次选取预设厚度的三维图像数据块，直至将整个目标三维图像划分完毕；其中，相邻三维图像数据块有一定重叠；

具体地，目标三维图像的宽为S(即有S层剖面)，已知三维图像数据块厚度层数H为10，相邻数据块重叠层数overlap为2，起始平面为最中间层记作Ss(Ss＝[S/2])；

首先，获取第1层至第Ss-1层之间的三维图像数据块，每个数据块包含的层数依次为第Ss-10到第Ss-1层、第Ss-18到第Ss-9层、第Ss-26到第Ss-17层、……、第1层到第10层。

为方便描述分别记为{Ss-10,Ss-1}、{Ss-18,Ss-9}、{Ss-26,Ss-17}、……、{1,10}。

其次，获取第Ss+1至第S层之间的三维图像数据块，每个数据块包含的层数依次为第Ss+1到第Ss+10层、第Ss+9到第Ss+18层、第Ss+17到第Ss+26层、……、第Ss-9层到第S层。

为方便描述分别记为{Ss+1,Ss+10}、{Ss+9,Ss+18}、{Ss+17,Ss+26}、……、{Ss-9,S}。

S5、以靠近起始平面的三维图像数据块开始，结合起始平面的边界分割线，依次从内向外计算各三维图像数据块的皮层边界，并组成稀疏的皮层表面；

具体地，上述步骤S5包括以下步骤：

S51、将预设厚度为H的三维图像数据块上每一层二维图像上对应的像素点的像素值取平均，得到目标二维图像：

其中，I(k，l，h)为三维图像数据块上各层二维图像的像素点的像素值，二维图像的层数为H；

S52、对目标二维图像依次执行上述步骤S31和S32，得到候选边界线L1；

S53、基于与起始平面所处的三维图像数据块的相邻三维图像数据块，根据起始平面的边界分割线起始点以及相邻三维图像数据块对应的目标二维图像转换的路径成本图，得到成本综合最小的从左往右的连续路径作为候选边界线L2；以此类推，得到各个三维图像数据块的候选边界线L2；

S54、将同一三维图像数据块对应的候选边界线L1与候选边界线L2融合得到三维图像数据块的边界分割线；

上述步骤S54，包括：

判断候选边界线L1与候选边界线L2同一横坐标的两个点(x，y₁)、(x，y₂)的纵坐标的差值是否处于预设范围内；

若是，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为(y₁+y₂)/2；

若否，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为y：

其中，y₀为相邻三维图像数据块的边界分割线的同一横坐标的点的纵坐标；

确定所有横坐标点对应的纵坐标，经过平滑处理得到三维图像数据块的边界分割线。

S55、将所有三维图像数据块的边界分割线汇合，得到稀疏的皮层表面。

S6、对稀疏的皮层表面进行平滑处理，得到稀疏的目标皮层表面；

具体地，上述步骤S6包括：对稀疏的皮层表面进行中值滤波；

另外，还可以计算中值滤波后的皮层表面的各像素点的二阶导数，取绝对值得到二维数据L；

从上往下依次判断二维数据L中每行数值，如果不大于阈值(例如10)，则保持不变；如果大于阈值，则取前后两行像素点的坐标值均值作为该像素点的坐标；实现进一步平滑。

S7、对稀疏的目标皮层表面扩充得到密集的三维皮层表面。

上述步骤S7包括：

将三维图像数据块的边界分割线作为其中间层的边界分割线；

基于相邻两个三维图像数据块的中间层的边界分割线，两中间层之间其他层的边界分割线通过两中间层的边界分割线等差递增或等差递减获得；

具体的等差递增或等差递减根据两中间层的边界分割线的各点的坐标值进行等差递增或等差递减确定其他层的边界线即可。例如，两中间层的边界分割线同一横坐标点的纵坐标分别为1和9，两中间层之间还有三层，则三层的同一横坐标点的纵坐标依次为3、5和7。

将所有三维图像数据块的所有层的边界分割线经过平滑得到密集、完整、平滑的三维皮层表面。

如图2所示，经过本实施例的三维皮层表面分割方法所得的皮层表面连续、平整且精确；而传统逐层分割皮层表面，十分耗时，受到皮层上方杂质和噪声影响，所得表面十分粗糙不连续。其中，传统逐层分割皮层表面的过程为对上述步骤S1-S2得到的目标三维图像逐层剖面，每一层剖面进行上述步骤S3得到每一层对应的边界分割线，最终组成三维皮层表面，具体边界分割线的分割过程可以参考上述步骤S3的具体过程。

本发明还提供OCTA图像的三维皮层表面分割系统，应用本实施例上述的三维皮层表面分割方法。

如图3所示，本实施例的三维皮层表面分割系统包括双边滤波模块、盖玻片去除模块、皮层边界分割模块、图像划分模块、皮层表面构建模块、平滑模块和扩充模块。

本实施例的双边滤波模块用于对三维OCTA图像进行双边滤波，得到滤波后图像。

具体地，利用双边滤波器对三维OCTA图像进行滤波，亮度差参数选定为5，空间距离参数选定为75，以突出图像中的边界信息。

本实施例的盖玻片去除模块用于去除滤波后图像中的盖玻片，得到目标三维图像。

具体地，去除滤波后图像中的盖玻片的过程包括：

(1)利用最大类间方差法将滤波后图像转换为0、1表示的掩膜；

(2)依次使用开运算、关运算和腐蚀运算对掩膜进行处理，得到目标掩膜；

(3)间隔选取目标掩膜的数个二维纵向剖面，以横坐标从左往右的顺序，依次从上往下开始判断每列对应像素点的掩膜数据；若为1，则提取相应目标像素点的坐标位置信息，开始下一列判断；若为0，则继续往下判断直至判断的像素点数量超过目标数量阈值，开始下一列判断；

(4)根据所有目标像素点的坐标位置信息拟合二阶曲面，曲面上方的所有杂质数据均剔除，得到剔除上层边缘面的图像；由于盖玻片为规则物体，因此拟合二阶曲面即可剔除杂质数据；

(5)对剔除上层边缘面的图像执行上述剔除上层边缘面的步骤(1)-(4)以剔除下层边缘面，得到目标三维图像。由于盖玻片具有一定厚度，所以剔除盖玻片的上层边缘面后，按照上述相同的步骤即可去除盖玻片的下层边缘面。

本实施例的皮层边界分割模块用于以目标三维图像正中间的二维纵向剖面作为起始平面，分割皮层边界，得到起始平面的边界分割线。

具体地，分割皮层边界的过程包括以下步骤：

(a)将起始平面的二维灰度图像转换为路径成本图；

(b)基于路径成本图和动态规划算法得到成本综合最小的从左往右的连续路径，作为边界分割线。

其中，上述步骤(a)包括以下步骤：

(一)将起始平面的二维灰度图像与核为3的y轴方向Sobel算子卷积，得到图像A；

将起始平面的二维灰度图像经过核为5的Canny算子提取边缘，再与长为3、宽为1的矩形结构化元素进行开运算，得到图像B；

(二)将图像A和图像B转换为路径成本图C：

C＝1-[ω*normalize(A)+(1-ω)*normalize(B)]

其中，ω为图像A的权重，1-ω为图像B的权重，normalize(A)为对图像A进行归一化至[0，1]区间，normalize(B)为对图像B进行归一化至[0，1]区间。ω优选取值为0.5。

上述步骤(b)具体包括：

设路径成本图C的长为M，高为N，各像素点(i，j)的成本记为C(i，j)，i＝1,2,3,…,M，j＝1,2,3,…,N；

以左端的N个像素点分别作为候选的N个起始点，从左到右依次寻找路径的下一个像素点，直至找到路径右端的像素点；其中，路径中相邻两像素点的成本分别记为C(i'，j')、C(i'+1，j”)，需满足以下条件：min C(i'+1，j”)，j”∈[j'-1，j'+1]；i'依次取值为1,2,3,…,M-1；

由于存在候选的N个起始点，相应的路径也有N条；

分别统计各路径所有像素点的成本总和，筛选得到成本综合最小的目标路径；

采用Savitzky-Golay卷积平滑算子对目标路径进行平滑处理，得到边界分割线。

其中，得到目标路径之后，还可以从右往左回溯目标路径(回溯过程可参考上述相邻像素点的寻找过程)，直至回溯完成，得到目标路径；如此设计，在前统计的N条路径无需存储相应的像素点的位置信息，节约内存。

本实施例的图像划分模块用于以起始平面为基准，往前后两侧依次选取预设厚度的三维图像数据块，直至将整个目标三维图像划分完毕；其中，相邻三维图像数据块有一定重叠。

具体地，目标三维图像的宽为S(即有S层剖面)，已知三维图像数据块厚度层数H为10，相邻数据块重叠层数overlap为2，起始平面为最中间层记作Ss(Ss＝[S/2])；

首先，获取第1层至第Ss-1层之间的三维图像数据块，每个数据块包含的层数依次为第Ss-10到第Ss-1层、第Ss-18到第Ss-9层、第Ss-26到第Ss-17层、……、第1层到第10层。

为方便描述分别记为{Ss-10,Ss-1}、{Ss-18,Ss-9}、{Ss-26,Ss-17}、……、{1,10}。

其次，获取第Ss+1至第S层之间的三维图像数据块，每个数据块包含的层数依次为第Ss+1到第Ss+10层、第Ss+9到第Ss+18层、第Ss+17到第Ss+26层、……、第Ss-9层到第S层。

为方便描述分别记为{Ss+1,Ss+10}、{Ss+9,Ss+18}、{Ss+17,Ss+26}、……、{Ss-9,S}。

本实施例的皮层表面构建模块用于以靠近起始平面的三维图像数据块开始，结合起始平面的边界分割线，依次从内向外计算各三维图像数据块的皮层边界，并组成稀疏的皮层表面。

具体地，将预设厚度为H的三维图像数据块上每一层二维图像上对应的像素点的像素值取平均，得到目标二维图像：

其中，I(k，l，h)为三维图像数据块上各层二维图像的像素点的像素值，二维图像的层数为H；

对目标二维图像依次执行上述步骤(a)和(b)，得到候选边界线L1；

基于与起始平面所处的三维图像数据块的相邻三维图像数据块，根据起始平面的边界分割线起始点以及相邻三维图像数据块对应的目标二维图像转换的路径成本图，得到成本综合最小的从左往右的连续路径作为候选边界线L2；以此类推，得到各个三维图像数据块的候选边界线L2；

将同一三维图像数据块对应的候选边界线L1与候选边界线L2融合得到三维图像数据块的边界分割线；

上述融合的具体过程包括：

判断候选边界线L1与候选边界线L2同一横坐标的两个点(x，y₁)、(x，y₂)的纵坐标的差值是否处于预设范围内；

若是，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为(y₁+y₂)/2；

若否，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为y：

其中，y₀为相邻三维图像数据块的边界分割线的同一横坐标的点的纵坐标；

确定所有横坐标点对应的纵坐标，经过平滑处理得到三维图像数据块的边界分割线。

将所有三维图像数据块的边界分割线汇合，得到稀疏的皮层表面。

本实施例的平滑模块用于对稀疏的皮层表面进行平滑处理，得到稀疏的目标皮层表面。

具体地，对稀疏的皮层表面进行中值滤波；

另外，还可以计算中值滤波后的皮层表面的各像素点的二阶导数，取绝对值得到二维数据L；

从上往下依次判断二维数据L中每行数值，如果不大于阈值(例如10)，则保持不变；如果大于阈值，则取前后两行像素点的坐标值均值作为该像素点的坐标；实现进一步平滑。

本实施例的扩充模块用于对稀疏的目标皮层表面扩充得到密集的三维皮层表面。

具体地，将三维图像数据块的边界分割线作为其中间层的边界分割线；

基于相邻两个三维图像数据块的中间层的边界分割线，两中间层之间其他层的边界分割线通过两中间层的边界分割线等差递增或等差递减获得；

具体的等差递增或等差递减根据两中间层的边界分割线的各点的坐标值进行等差递增或等差递减确定其他层的边界线即可。例如，两中间层的边界分割线同一横坐标点的纵坐标分别为1和9，两中间层之间还有三层，则三层的同一横坐标点的纵坐标依次为3、5和7。

将所有三维图像数据块的所有层的边界分割线经过平滑得到密集、完整、平滑的三维皮层表面。

本发明的OCTA图像的三维皮层表面分割系统，实现自动分割皮层表面，且三维皮层表面分割精度高。

实施例2：

本实施例的OCTA图像的三维皮层表面分割方法与实施例1的不同之处在于：

在对稀疏的皮层表面进行平滑处理只采用中值滤波进行平滑即可，无需计算二阶导数，满足不同应用的需求；

其他步骤可以参考实施例1；

相应地，本实施例的OCTA图像的三维皮层表面分割系统的平滑模块也只需采用中值滤波进行平滑即可；

其他构架可以参考实施例1。

以上所述仅是对本发明的优选实施例及原理进行了详细说明，对本领域的普通技术人员而言，依据本发明提供的思想，在具体实施方式上会有改变之处，而这些改变也应视为本发明的保护范围。

Claims (2)

1.OCTA图像的三维皮层表面分割方法，其特征在于，包括以下步骤：

S1、对三维OCTA图像进行双边滤波，得到滤波后图像；

S2、去除滤波后图像中的盖玻片，得到目标三维图像；

S3、以目标三维图像正中间的二维纵向剖面作为起始平面，分割皮层边界，得到起始平面的边界分割线；

S4、以起始平面为基准，往前后两侧依次选取预设厚度的三维图像数据块，直至将整个目标三维图像划分完毕；其中，相邻三维图像数据块有一定重叠；

S5、以靠近起始平面的三维图像数据块开始，结合起始平面的边界分割线，依次从内向外计算各三维图像数据块的皮层边界，并组成稀疏的皮层表面；

S6、对稀疏的皮层表面进行平滑处理，得到稀疏的目标皮层表面；

S7、对稀疏的目标皮层表面扩充得到密集的三维皮层表面；

所述步骤S2，包括以下步骤：

S21、利用最大类间方差法将滤波后图像转换为0、1表示的掩膜；

S22、依次使用开运算、关运算和腐蚀运算对掩膜进行处理，得到目标掩膜；

S23、间隔选取目标掩膜的数个二维纵向剖面，以横坐标从左往右的顺序，依次从上往下开始判断每列对应像素点的掩膜数据；若为1，则提取相应目标像素点的坐标位置信息，开始下一列判断；若为0，则继续往下判断直至判断的像素点数量超过目标数量阈值，开始下一列判断；

S24、根据所有目标像素点的坐标位置信息拟合二阶曲面，曲面上方的所有杂质数据均剔除，得到剔除上层边缘面的图像；

S25、对剔除上层边缘面的图像执行步骤S21-S24以剔除下层边缘面，得到目标三维图像；

所述步骤S3中，分割皮层边界的过程包括以下步骤：

S31、将起始平面的二维灰度图像转换为路径成本图；

S32、基于路径成本图和动态规划算法得到成本综合最小的从左往右的连续路径，作为边界分割线；

所述步骤31包括以下步骤：

S311、将起始平面的二维灰度图像与核为3的y轴方向Sobel算子卷积，得到图像A；

将起始平面的二维灰度图像经过核为5的Canny算子提取边缘，再与长为3、宽为1的矩形结构化元素进行开运算，得到图像B；

S312、将图像A和图像B转换为路径成本图C：

C＝1-[ω*normalize(A)+(1-ω)*normalize(B)]

其中，ω为图像A的权重，1-ω为图像B的权重，normalize(A)为对图像A进行归一化至[0，1]区间，normalize(B)为对图像B进行归一化至[0，1]区间；

所述步骤S32包括：

设路径成本图C的长为M，高为N，各像素点(i，j)的成本记为C(i，j)，i＝1,2,3,…,M，j＝1,2,3,…,N；

以左端的N个像素点分别作为候选的N个起始点，从左到右依次寻找路径的下一个像素点，直至找到路径右端的像素点；其中，路径中相邻两像素点的成本分别记为C(i'，j')、C(i'+1，j”)，需满足以下条件：min C(i'+1，j”)，j”∈[j'-1，j'+1]；i'依次取值为1,2,3,…,M-1；

分别统计各路径所有像素点的成本总和，筛选得到成本综合最小的目标路径；

采用Savitzky-Golay卷积平滑算子对目标路径进行平滑处理，得到边界分割线；

所述步骤S5，包括以下步骤：

S51、将预设厚度为H的三维图像数据块上每一层二维图像上对应的像素点的像素值取平均，得到目标二维图像：

其中，I(k，l，h)为三维图像数据块上各层二维图像的像素点的像素值，二维图像的层数为H；

S52、对目标二维图像依次执行所述步骤S31和S32，得到候选边界线L1；

S53、基于与起始平面所处的三维图像数据块的相邻三维图像数据块，根据起始平面的边界分割线起始点以及相邻三维图像数据块对应的目标二维图像转换的路径成本图，得到成本综合最小的从左往右的连续路径作为候选边界线；以此类推，得到各个三维图像数据块的候选边界线L2；

S54、将同一三维图像数据块对应的候选边界线L1与候选边界线L2融合得到三维图像数据块的边界分割线；

S55、将所有三维图像数据块的边界分割线汇合，得到稀疏的皮层表面；

所述步骤S54，包括：

判断候选边界线L1与候选边界线L2同一横坐标的两个点(x，y₁)、(x，y₂)的纵坐标的差值是否处于预设范围内；

若是，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为(y₁+y₂)/2；

若否，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为y：

其中，y₀为相邻三维图像数据块的边界分割线的同一横坐标的点的纵坐标；

所述步骤S6中，对稀疏的皮层表面进行中值滤波；

所述步骤S7，包括：

将三维图像数据块的边界分割线作为其中间层的边界分割线；

基于相邻两个三维图像数据块的中间层的边界分割线，两中间层之间其他层的边界分割线通过两中间层的边界分割线等差递增或等差递减获得；

将所有三维图像数据块的所有层的边界分割线经过平滑得到密集的三维皮层表面。

2.OCTA图像的三维皮层表面分割系统，应用如权利要求1所述的三维皮层表面分割方法，其特征在于，所述三维皮层表面分割系统包括：

双边滤波模块，用于对三维OCTA图像进行双边滤波，得到滤波后图像；

盖玻片去除模块，用于去除滤波后图像中的盖玻片，得到目标三维图像；具体包括以下步骤：

S21、利用最大类间方差法将滤波后图像转换为0、1表示的掩膜；

S22、依次使用开运算、关运算和腐蚀运算对掩膜进行处理，得到目标掩膜；

S23、间隔选取目标掩膜的数个二维纵向剖面，以横坐标从左往右的顺序，依次从上往下开始判断每列对应像素点的掩膜数据；若为1，则提取相应目标像素点的坐标位置信息，开始下一列判断；若为0，则继续往下判断直至判断的像素点数量超过目标数量阈值，开始下一列判断；

S24、根据所有目标像素点的坐标位置信息拟合二阶曲面，曲面上方的所有杂质数据均剔除，得到剔除上层边缘面的图像；

S25、对剔除上层边缘面的图像执行步骤S21-S24以剔除下层边缘面，得到目标三维图像；

皮层边界分割模块，用于以目标三维图像正中间的二维纵向剖面作为起始平面，分割皮层边界，得到起始平面的边界分割线；分割皮层边界的过程包括以下步骤：

S31、将起始平面的二维灰度图像转换为路径成本图；

S32、基于路径成本图和动态规划算法得到成本综合最小的从左往右的连续路径，作为边界分割线；

所述步骤31包括以下步骤：

S311、将起始平面的二维灰度图像与核为3的y轴方向Sobel算子卷积，得到图像A；

将起始平面的二维灰度图像经过核为5的Canny算子提取边缘，再与长为3、宽为1的矩形结构化元素进行开运算，得到图像B；

S312、将图像A和图像B转换为路径成本图C：

C＝1-[ω*normalize(A)+(1-ω)*normalize(B)]

其中，ω为图像A的权重，1-ω为图像B的权重，normalize(A)为对图像A进行归一化至[0，1]区间，normalize(B)为对图像B进行归一化至[0，1]区间；

所述步骤S32包括：

设路径成本图C的长为M，高为N，各像素点(i，j)的成本记为C(i，j)，i＝1,2,3,…,M，j＝1,2,3,…,N；

以左端的N个像素点分别作为候选的N个起始点，从左到右依次寻找路径的下一个像素点，直至找到路径右端的像素点；其中，路径中相邻两像素点的成本分别记为C(i'，j')、C(i'+1，j”)，需满足以下条件：min C(i'+1，j”)，j”∈[j'-1，j'+1]；i'依次取值为1,2,3,…,M-1；

分别统计各路径所有像素点的成本总和，筛选得到成本综合最小的目标路径；

采用Savitzky-Golay卷积平滑算子对目标路径进行平滑处理，得到边界分割线；

图像划分模块，用于以起始平面为基准，往前后两侧依次选取预设厚度的三维图像数据块，直至将整个目标三维图像划分完毕；其中，相邻三维图像数据块有一定重叠；

皮层表面构建模块，用于以靠近起始平面的三维图像数据块开始，结合起始平面的边界分割线，依次从内向外计算各三维图像数据块的皮层边界，并组成稀疏的皮层表面；具体包括以下步骤：

S51、将预设厚度为H的三维图像数据块上每一层二维图像上对应的像素点的像素值取平均，得到目标二维图像：

其中，I(k，l，h)为三维图像数据块上各层二维图像的像素点的像素值，二维图像的层数为H；

S52、对目标二维图像依次执行所述步骤S31和S32，得到候选边界线L1；

S53、基于与起始平面所处的三维图像数据块的相邻三维图像数据块，根据起始平面的边界分割线起始点以及相邻三维图像数据块对应的目标二维图像转换的路径成本图，得到成本综合最小的从左往右的连续路径作为候选边界线；以此类推，得到各个三维图像数据块的候选边界线L2；

S54、将同一三维图像数据块对应的候选边界线L1与候选边界线L2融合得到三维图像数据块的边界分割线；

S55、将所有三维图像数据块的边界分割线汇合，得到稀疏的皮层表面；

所述步骤S54，包括：

判断候选边界线L1与候选边界线L2同一横坐标的两个点(x，y₁)、(x，y₂)的纵坐标的差值是否处于预设范围内；

若是，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为(y₁+y₂)/2；

若否，则对应的三维图像数据块的边界分割线上同一横坐标的点的纵坐标为y：

其中，y₀为相邻三维图像数据块的边界分割线的同一横坐标的点的纵坐标；

平滑模块，用于对稀疏的皮层表面进行平滑处理，得到稀疏的目标皮层表面；其中，对稀疏的皮层表面进行中值滤波；

扩充模块，用于对稀疏的目标皮层表面扩充得到密集的三维皮层表面，具体包括：

将三维图像数据块的边界分割线作为其中间层的边界分割线；

基于相邻两个三维图像数据块的中间层的边界分割线，两中间层之间其他层的边界分割线通过两中间层的边界分割线等差递增或等差递减获得；

将所有三维图像数据块的所有层的边界分割线经过平滑得到密集的三维皮层表面。

Images