本发明属于目标分割领域，特别涉及一种基于时间序列sd-oct视网膜图像的ga自动分割方法。
背景技术
sd-oct视网膜图像是一种频域光学相干断层成像图像，它可以有效地呈现视网膜组织层的灰度和结构等变化，临床实验表明sd-oct图像能够用于测量地图状萎缩(ga)的范围。由于在sd-oct图像中ga的范围需要大量的人工手动标注，非常耗时，因此需要一种自动分割ga的方法来快速有效的判断ga的范围。现有的ga分割方法主要是基于二维限制投影图像，通过迭代阈值分割方法和强度分布集来估计ga区域；或者在三维sd-oct图像中通过构建具有深度神经网络的投票系统来自动分割ga。近几年主要出现了如下两种ga自动分割方法：
(1)基于二维限制投影图像分割的方法。该方法通过将三维sd-oct图像限制投影生成二维图像，通过迭代阈值分割方法和强度分布集来估计ga区域，该方法没有排除在三维图像中所限制区域内与ga亮度相似的区域，因此在二维投影图像中对ga的分割造成了影响。
(2)基于深度神经网络的分割的方法。该方法使用具有投票决策策略的vgg16卷积神经网络作为用于ga分割的深度学习模型，但该方法需要大量的训练样本，而且在投票之前需要训练10个深度神经网络模型，非常耗时，不利于医生在少量数据及短时间内判断ga的范围。
综上，现有的ga自动分割方法没有考虑与ga相似区域的影响以及需要在少量数据和短时间内完成ga分割的情况，导致ga的自动分割准确率和效率不高。
技术实现要素：
本发明的目的在于提供一种基于时间序列sd-oct视网膜图像的ga自动分割方法。
实现本发明目的的技术解决方案为：一种基于时间序列sd-oct视网膜图像的ga自动分割方法，包括以下步骤：
步骤1、采集对象i的初始三维sd-oct视网膜图像a，并对三维sd-oct视网膜图像a进行图像处理，建立ga分割模型；
步骤2、采集某一时刻对象i的待分割ga的三维sd-oct视网膜图像b，然后对三维sd-oct视网膜图像b进行步骤1所述的图像处理，并通过步骤1建立的ga分割模型对待分割ga的三维sd-oct视网膜图像b进行ga分割。
本发明与现有技术相比，其显著优点：1)本发明的方法充分利用了不同时间三维sd-oct视网膜图像的ga的相关性，依据目标跟踪的思想提高了ga自动分割的准确率和效率；2)本发明的方法考虑了与ga相似区域的影响，重新获取训练样本的过程使得ga特征在图像中更加显著，避免了三维sd-oct视网膜图像在采集过程中位置偏移对分割结果产生影响，分割结果更加准确；3)本发明的方法只需少量训练数据即可以在短时间内准确分割出ga，分割速度快，效率高。
下面结合附图对本发明作进一步详细描述。
附图说明
图1为本发明基于时间序列sd-oct视网膜图像的ga自动分割方法的流程图。
图2为本发明中图像处理的流程图。
图3为本发明中获取b-scan图像感兴趣区域的流程图。
图4为本发明中获取训练样本的流程图。
图5为本发明中提取训练样本的方向梯度直方图特征的流程图。
图6为本发明实施例中采集的对象i的初始三维sd-oct视网膜图像a'的示意图。
图7为本发明实施例中三维sd-oct视网膜图像a'的内部组织示意图。
图8为本发明实施例中三维sd-oct视网膜图像a'的b-scan图像的ilm层和ilm层向下平移之后形成上下边界的示意图。
图9为本发明实施例中三维sd-oct视网膜图像a'的b-scan图像的感兴趣区域拉平后感兴趣区域图像示意图。
图10为本发明实施例中建立的ga分割模型示意图。
图11为本发明实施例中采集的5月份对象i的三维sd-oct视网膜图像b'的示意图。
图12为本发明实施例中三维sd-oct视网膜图像b'的内部组织示意图。
图13为本发明实施例中三维sd-oct视网膜图像b'的b-scan图像的ilm层和ilm层向下平移之后形成上下边界的示意图。
图14为本发明实施例中三维sd-oct视网膜图像b'的b-scan图像的感兴趣区域拉平后感兴趣区域图像示意图。
图15为本发明实施例中三维sd-oct视网膜图像b'在其b-scan图像上ga自动分割结果与金标准对比示意图。
图16为本发明实施例中三维sd-oct视网膜图像b'在投影图像上ga自动分割结果与金标准对比示意图。
具体实施方式
结合图1，本发明的基于时间序列sd-oct视网膜图像的ga自动分割方法包括以下步骤：
步骤1、采集对象i的初始三维sd-oct视网膜图像a，并对三维sd-oct视网膜图像a进行图像处理，建立ga分割模型。
其中，图像处理具体为：
步骤1-1、对三维sd-oct视网膜图像进行横向扫描，获取n帧b-scan图像。
步骤1-2、对步骤1-1获得的n帧b-scan图像均进行双边滤波去噪，并通过最优三维图搜索方法分割每一帧b-scan图像各自的ilm层和bm层。
步骤1-3、获取步骤1-2中每一帧b-scan图像的感兴趣区域。具体为：
步骤1-3-1、将每一帧b-scan图像的bm层向下平移100个像素，获得新的bm’边界；
步骤1-3-2、将每一帧b-scan图像的ilm层作为上边界，求取上边界与步骤1-3-1获得的bm’边界之间的平均距离l；
步骤1-3-3、将每一帧b-scan图像的ilm层向下平移步骤1-3-2获得的平均距离l，获得新的ilm’边界，将其作为下边界；
步骤1-3-4、将步骤1-3-2的上边界与步骤1-3-3的下边界之间的区域作为每一帧b-scan图像的感兴趣区域。
步骤1-4、将步骤1-3获得的每一帧b-scan图像的感兴趣区域拉平，并以窗口滑动的方法获取三维sd-oct视网膜图像的m个训练样本。具体为：
步骤1-4-1、将步骤1-3-2中的上边界作为第一行，步骤1-3-3中的下边界作为最后一行，由此将每一帧b-scan图像的感兴趣区域拉平，获得每一帧b-scan图像的感兴趣区域图像；
步骤1-4-2、通过w×h的滑动窗口以步长l'遍历步骤1-4-1获得的每一帧b-scan图像的感兴趣区域图像，将每一帧b-scan图像的感兴趣区域图像分为m个训练样本，由此获得三维sd-oct视网膜图像的m个训练样本；其中w为滑动窗口的宽度，h为每一帧b-scan图像感兴趣区域图像的高度。
其中，m的计算公式为：
m＝(w/l')-(w/l'-1)
式中，w为每一帧b-scan的感兴趣区域图像的宽度。
其中，m的计算公式为：
m＝m×n
式中，n为n帧b-scan图像。
步骤1-5、对步骤1-4得到的m个训练样本的方向梯度直方图特征进行提取。具体为：
步骤1-5-1、对步骤1-4得到的m个训练样本进行gamma标准化，获得m个训练样本标准化图像；
步骤1-5-2、求取步骤1-5-1获得的m个训练样本标准化图像中每个像素梯度的幅值和方向；
每个像素梯度的幅值：
每个像素梯度的方向：
式中，f(x,y)为训练样本标准化图像在(x,y)处的像素值；
步骤1-5-3、将步骤1-5-1中的每一个训练样本标准化图像分成n个s×s单元格；其中，所述n与s的关系为：
式中，h'为训练样本标准化图像的高度，w'为训练样本标准化图像的宽度，n为单元格的个数，s为单元格的边长；
步骤1-5-4、为步骤1-5-3的每个单元格构建方向梯度直方图，统计每个单元格的梯度信息；
步骤1-5-5、将p×p个单元格组成一个block块，根据步骤1-5-4统计的每个单元格的梯度信息求取每一个block块的方向梯度直方图向量，并通过l2范数对其进行归一化，由此获得训练样本的方向梯度直方图特征。其中，所述p的计算公式为:
p＝a×s
式中，a为偶数且
其中，建立ga分割模型具体为：通过随机森林算法对步骤1-5提取的方向梯度直方图特征进行训练，获取随机森林模型，将其作为ga分割模型。
步骤2、采集某一时刻对象i的待分割ga的三维sd-oct视网膜图像b，然后对三维sd-oct视网膜图像b进行步骤1所述的图像处理，并通过步骤1建立的ga分割模型对待分割ga的三维sd-oct视网膜图像b进行ga分割。
下面结合实施例对本发明做进一步详细的说明。
实施例
本系统发明以sd-oct视网膜图像作为输入，采用图像处理与机器学习手段自动分割ga区域。结合图1，本发明一种基于时间序列sd-oct视网膜图像的ga自动分割方法，包括以下步骤：
步骤1、采集对象i的初始三维sd-oct视网膜图像a'如图6所示，并对三维sd-oct视网膜图像a'进行图像处理，结合图2具体为：
步骤1-1、对三维sd-oct视网膜图像a'进行横向扫描，获取n帧b-scan图像。本实施例中，n＝128，如图6所示。
步骤1-2、对步骤1-1获得的128帧b-scan图像均进行双边滤波去噪，并通过最优三维图搜索方法分割每一帧b-scan图像各自的ilm层和bm层，如图7所示。
步骤1-3、结合图3，获取步骤1-2中每一帧b-scan图像的感兴趣区域。具体为：
步骤1-3-1、将每一帧b-scan图像的bm层向下平移100个像素，获得新的bm’边界；
步骤1-3-2、将每一帧b-scan图像的ilm层作为上边界，求取上边界与步骤1-3-1获得的bm’边界之间的平均距离l；
步骤1-3-3、将每一帧b-scan图像的ilm层向下平移步骤1-3-2获得的平均距离l，获得新的ilm’边界，将其作为下边界；
步骤1-3-4、将步骤1-3-2的上边界与步骤1-3-3的下边界之间的区域作为每一帧b-scan图像的感兴趣区域。
本实施例中获取三维sd-oct视网膜图像a'的第70帧b-scan图像的感兴趣区域如图8所示。
步骤1-4、结合图4，将步骤1-3获得的每一帧b-scan图像的感兴趣区域拉平，并以窗口滑动的方法获取三维sd-oct视网膜图像的m个训练样本。具体为：
步骤1-4-1、将步骤1-3-2中的上边界作为第一行，步骤1-3-3中的下边界作为最后一行，由此将每一帧b-scan图像的感兴趣区域拉平，获得每一帧b-scan图像的感兴趣区域图像。本实施例中获取三维sd-oct视网膜图像a'的第70帧b-scan图像的感兴趣区域图像如图9所示。
步骤1-4-2、通过w×h的滑动窗口以步长l'遍历步骤1-4-1获得的每一帧b-scan图像的感兴趣区域图像，将每一帧b-scan图像的感兴趣区域图像分为m个训练样本，由此获得三维sd-oct视网膜图像的m个训练样本；其中w为滑动窗口的宽度，h为每一帧b-scan图像感兴趣区域图像的高度。
其中，m的计算公式为：
m＝(w/l')-(w/l'-1)
式中，w为每一帧b-scan的感兴趣区域图像的宽度。本实施例中，w＝512，w＝64，l'＝8，因此可得m＝57。
其中，m的计算公式具体为：
m＝m×n
式中，n为n帧b-scan图像。本实施例中，m＝57×128＝7296。
步骤1-5、结合图5，对步骤1-4得到的m个训练样本的方向梯度直方图特征进行提取。具体为：
步骤1-5-1、对步骤1-4得到的m个训练样本进行gamma标准化，获得m个训练样本标准化图像；
步骤1-5-2、求取步骤1-5-1获得的m个训练样本标准化图像中每个像素梯度的幅值和方向；
每个像素梯度的幅值：
每个像素梯度的方向：
式中，f(x,y)为训练样本标准化图像在(x,y)处的像素值；
步骤1-5-3、将步骤1-5-1中的每一个训练样本标准化图像分成n个s×s单元格；其中，n与s的关系为：
式中，h'为训练样本标准化图像的高度，w'为训练样本标准化图像的宽度，n为单元格的个数，s为单元格的边长。本实施例中取s＝8，为方便特征提取，将样本高度调整为统一高度h'＝128，而w'＝64，因此n＝128。
步骤1-5-4、为步骤1-5-3的每个单元格构建方向梯度直方图，统计每个单元格的梯度信息；
步骤1-5-5、将p×p个单元格组成一个block块，根据步骤1-5-4统计的每个单元格的梯度信息求取每一个block块的方向梯度直方图向量，并通过l2范数对其进行归一化，由此获得训练样本的方向梯度直方图特征。其中，p的计算公式为:
p＝a×s
式中，a为偶数且
建立的ga分割模型如图10所示，包含了分割模型的相关参数。
步骤2、采集某一时刻对象i的待分割ga的三维sd-oct视网膜图像b'如图11所示，本实例中选取的某一时刻为5月份。然后对三维sd-oct视网膜图像b'进行步骤1所述的图像处理，结合图2具体为：
步骤1-1、对三维sd-oct视网膜图像b'进行横向扫描，获取n帧b-scan图像。本实施例中，n＝128，如图11所示。
步骤1-2、对步骤1-1获得的128帧b-scan图像均进行双边滤波去噪，并通过最优三维图搜索方法分割每一帧b-scan图像各自的ilm层和bm层，如图12所示。
步骤1-3、结合图3，获取步骤1-2中每一帧b-scan图像的感兴趣区域。具体为：
步骤1-3-1、将每一帧b-scan图像的bm层向下平移100个像素，获得新的bm’边界；
步骤1-3-2、将每一帧b-scan图像的ilm层作为上边界，求取上边界与步骤1-3-1获得的bm’边界之间的平均距离l；
步骤1-3-3、将每一帧b-scan图像的ilm层向下平移步骤1-3-2获得的平均距离l，获得新的ilm’边界，将其作为下边界；
步骤1-3-4、将步骤1-3-2的上边界与步骤1-3-3的下边界之间的区域作为每一帧b-scan图像的感兴趣区域。
本实施例中获取三维sd-oct视网膜图像b'的第70帧b-scan图像的感兴趣区域如图13所示。
步骤1-4、结合图4，将步骤1-3获得的每一帧b-scan图像的感兴趣区域拉平，并以窗口滑动的方法获取三维sd-oct视网膜图像的m个训练样本。具体为：
步骤1-4-1、将步骤1-3-2中的上边界作为第一行，步骤1-3-3中的下边界作为最后一行，由此将每一帧b-scan图像的感兴趣区域拉平，获得每一帧b-scan图像的感兴趣区域图像。本实施例中获取三维sd-oct视网膜图像b'的第70帧b-scan图像的感兴趣区域图像如图14所示。
步骤1-4-2、通过w×h的滑动窗口以步长l'遍历步骤1-4-1获得的每一帧b-scan图像的感兴趣区域图像，将每一帧b-scan图像的感兴趣区域图像分为m个训练样本，由此获得三维sd-oct视网膜图像的m个训练样本；其中w为滑动窗口的宽度，h为每一帧b-scan图像感兴趣区域图像的高度。
其中，m的计算公式为：
m＝(w/l')-(w/l'-1)
式中，w为每一帧b-scan的感兴趣区域图像的宽度。本实施例中，w＝512，w＝64，l'＝8，因此可得m＝57。
其中，m的计算公式具体为：
m＝m×n
式中，n为n帧b-scan图像。本实施例中，m＝57×128＝7296。
步骤1-5、结合图5，对步骤1-4得到的m个训练样本的方向梯度直方图特征进行提取。具体为：
步骤1-5-1、对步骤1-4得到的m个训练样本进行gamma标准化，获得m个训练样本标准化图像；
步骤1-5-2、求取步骤1-5-1获得的m个训练样本标准化图像中每个像素梯度的幅值和方向；
每个像素梯度的幅值：
每个像素梯度的方向：
式中，f(x,y)为训练样本标准化图像在(x,y)处的像素值；
步骤1-5-3、将步骤1-5-1中的每一个训练样本标准化图像分成n个s×s单元格；其中，n与s的关系为：
式中，h'为训练样本标准化图像的高度，w'为训练样本标准化图像的宽度，n为单元格的个数，s为单元格的边长。本实施例中取s＝8，为方便特征提取，将样本高度调整为统一高度h'＝128，而w'＝64，因此n＝128。
步骤1-5-4、为步骤1-5-3的每个单元格构建方向梯度直方图，统计每个单元格的梯度信息；
步骤1-5-5、将p×p个单元格组成一个block块，根据步骤1-5-4统计的每个单元格的梯度信息求取每一个block块的方向梯度直方图向量，并通过l2范数对其进行归一化，由此获得训练样本的方向梯度直方图特征。其中，p的计算公式为:
p＝a×s
式中，a为偶数且
通过上述建立的ga分割模型如图10对待分割ga的三维sd-oct视网膜图像b'进行ga分割，得到三维sd-oct视网膜图像b'在其b-scan图像上ga自动分割结果与金标准对比如图15所示，得到三维sd-oct视网膜图像b'在投影图像上ga自动分割结果与金标准对比如图16所示。由图15、16可知，本发明ga自动分割结果与金标准具有高度的一致性、稳定性和鲁棒性。
本发明能够在短时间内训练少量样本的情况下自动的分割出ga，ga分割效率高、稳定性和鲁棒性好，对后续ga定量分析具有重要意义。