脑梗塞核心区变化的血管形态学指标的提取与筛选方法与流程

1.本发明涉及梗塞变化预测领域，具体涉及一种脑梗塞核心区变化的血管形态学指标的提取与筛选方法。


背景技术：

2.脑梗死旧称脑梗塞，又称缺血性脑卒中(cerebral ischemic stroke)，是指因脑部血液供应障碍，缺血、缺氧所导致的局限性脑组织的缺血性坏死或软化。脑梗死的临床常见类型有脑血栓形成、腔隙性梗死和脑栓塞等，脑梗死占全部脑卒中的80％，脑梗死不仅给人类健康和生命造成极大威胁，而且给患者、家庭及社会带来极大的痛苦和沉重的负担。
3.目前，对脑梗塞的识别方法，一般都是去医院进行ct检查、mri检查、常规检查或特殊检查。
4.在进行ct检查时，ct显示梗死灶为低密度，可以明确病变的部位、形状及大小，较大的梗死灶可使脑室受压，变形及中线结构移位，但脑梗死起病4～6小时内，只有部分病例可见边界不清的稍低密度灶，而大部分的病例在24小时后才能显示边界较清的低密度灶，且小于5mm的梗死灶。
5.mri检查就是核磁共振成像检查，易得性差，目前无法实现急诊mri检查和复查，临床应用有局限。
6.常规检查就是血、尿、大便常规及肝功能、肾功能、凝血功能、血糖、血脂、心电图等作为常规检查，特殊检查如经颅多普勒超声(tcd)、颈动脉彩色b超、磁共振、血管造影(mra)，数字减影全脑血管造影(dsa)、颈动脉造影，可明确有无颅内外动脉狭窄或闭塞。
7.梗死体积改变之前就开始出现脑血管形态学的变化，脑血管的形态学对血流特征产生影响，如何将临床医生对血管形态的直觉感知转化为特定指标目前尚没有较好的方法。通过分割和量化脑血管形态对预测急性卒中受试者的预后有重要临床意义。有鉴于此，本技术提出一种脑梗塞核心区变化的血管形态学指标的提取与筛选方法。


技术实现要素：

8.本发明的目的在于针对现有技术的不足，提供一种脑梗塞核心区变化的血管形态学指标的提取与筛选方法。
9.为了解决上述技术问题，采用如下技术方案：
10.脑梗塞核心区变化的血管形态学指标的提取与筛选方法，其特征在于包括以下步骤：
11.(1)ctp影像获取：先获取受试者的ctp影像结果。
12.(2)ctp影像血管分割：采用3d unet网络模型对所述ctp影像的脑血管进行分割，所述3d unet网络模型通过跳层连接和下采样层将不同尺度图像信息融合，实现所述脑血管分割。
13.(3)脑血管形态学特征提取：在分割所述脑血管后，进一步量化脑血管形态结构信
息，先提取所述脑血管的3d中心线，提取结果中包含脑血管中心线和各节点信息；然后提取脑血管形态结构特征，得到3d脑血管形态学特征；并根据所述3d脑血管形态学特征进一步提取2d脑血管形态学特征，由此得到受试者的整体脑血管形态学特征及梗塞周围区域脑血管形态学特征。
14.(4)脑血管形态学特征筛选：根据每个所述受试者的整体脑血管形态学特征及梗塞周围区域脑血管形态学特征，建立基于非参数训练的稀疏表示分类器进行分类预测，得到脑血管形态学的影像组学分类预测模型。
15.进一步，所述脑血管形态结构特征包括：血管体积，血管分支数，血管总长度，平均血管长度，平均距离因子df，平均角度soam，边结点个数，网络结点个数，聚类系数c，结构熵h和平均路径长度。
16.进一步，所述平均角度soam用于衡量血管的扭曲度，所述平均角度soam通过四个结点间的3个向量来度量血管3d空间的扭曲度；
17.令t1＝p
k-p
k-1
,t2＝p
k+1-pk,t3＝p
k+2-p
k+1
,则脑血管在pk点的弯曲度dpk和扭度cpk可以表示为：
[0018][0019][0020]
其中，
[0021][0022]
进一步，所述聚类系数c和结构熵h是网络拓扑图的度量参数，所述聚类系数c计算公式如下：
[0023][0024]
其中n表示节点个数，li表示节点i的邻节点个数，ki表示i与邻节点间的边个数；
[0025]
结构熵h计算公式如下：
[0026][0027]
其中表示第i个节点的度。
[0028]
进一步，在所述步骤(3)中，根据所述3d脑血管形态学特征进一步提取2d血管形态学特征的具体过程如下：首先将3d脑血管投影到2d轴状面；然后进行脑血管增强，2d阈值分割和提取血管中心线操作；最后提取2d血管的描述脑血管形态结构特征；通过提取大脑左右两侧血管结构的差异信息来进一步分析脑梗塞，将二者的差值作为对称两侧血管差异特征。
[0029]
进一步，在所述步骤(4)中，在建立得到所述脑血管形态学的影像组学分类预测模型前，通过对受试者的两次mri-dwi图像结果进行分割梗塞核心区，获取梗死核心区的体积，将前后两次对比，得到差值作为不同梗塞变化的分组依据，将所述梗塞变化分为梗塞逆转组和梗塞扩大组。
[0030]
进一步，在所述步骤(4)中，建立基于非参数训练的稀疏表示分类器进行分类预测，得到脑血管形态学的影像组学分类预测模型的过程如下：
[0031]
首先利用稀疏编码和重构计算待测试样与训练集里每类样本间的距离，然后根据k-邻近准则判断测试样本的类别分类过程具体可表示为：
[0032][0033]
其中表示待测样本特征选择后的特征，为训练样本特征集，c为样本总类别数；γ为稀疏表示控制参数；当获得稀疏表示系数后，计算待测试样与训练集里每类样本间的距离，即各类的表示残差：
[0034][0035]
其中δc(
·
)表示选择第c类特征对应的系数。残差最小的类别即为测试样本的类别，因此最终待测样本类别为
[0036]
进一步，在所述步骤(1)中，所述ctp影像的体素分辨率和大小分别为0.486*0.486*1mm3和512*512*319，所述ctp影像提取预测发生梗塞的切片图像进行实验，所述ctp影像的体素灰度由窗口[0,600]裁剪，所述ctp影像的体素灰度用于聚焦脑血管。
[0037]
进一步，在所述步骤(2)中，所述3d unet网络模型在训练阶段输入和输出的ctp的3d图像，进行编码和解码，编码阶段的特征图通过跳层连接传输到解码阶段进行图像特征的卷积融合，得到训练图像数据；然后抽取训练图像数据的血管和背景区域，进行模型训练；最后进行测试，得到整图分割结果。
[0038]
进一步，还包括步骤(5)脑血管形态学特征模型验证：所述脑血管形态学的影像组学分类预测模型建立之后，提取并筛选出受试者的脑血管形态学特征，纳入所述脑血管形态学的影像组学分类预测模型中进行预测，将预测结果与受试者7-10天的复查的mri-dwi图像图像结果进行比较，来评估预测的有效性和准确性。
[0039]
由于采用上述技术方案，具有以下有益效果：
[0040]
本发明为一种脑梗塞核心区变化的血管形态学指标的提取与筛选方法，先获取受试者的ctp影像结果。然后采用3d unet网络模型对cta图像脑血管进行分割。3d unet网络模型通过跳层连接和下采样层将不同尺度图像信息融合，从而实现精确的血管分割。在分割脑血管后，为进一步量化脑血管形态结构信息，提取脑血管的3d中心线，提取结果中包含血管中心线和各节点信息。然后提取脑血管形态结构特征。最终提取受试者的整体脑血管形态学特征及梗塞周围区域脑血管形态学特征。获取梗死核心区的体积，将前后两次对比，得到差值为结局分别进行分类，并纳入两次检查的时间因素。分别根据梗死核心区的体积的增减，拟使用非参数训练的稀疏表示分类器建立决策模型。明确脑血管形态特征因素对结局的影响。同时纳入临床因素进行再次统计，提高预测效能。由此通过脑血管网络化特征指标的建立和提取，可以准确地预测梗塞核心区变化。提取出的特征因子可以作为预测脑梗塞预后的潜在生物标志物。
[0041]
本方法在计算机断层扫描血管造影中分割脑血管，分析和挖掘2d和3d敏感脑血管结构特征，并有效提高细小血管的分割效果。一方面解决了脑血管具有复杂的3d结构分割困难的问题，另一方面解决了一些小血管因造影剂的扩散不均会出现大量连接断裂区域的问题。
附图说明
[0042]
下面结合附图对本发明作进一步说明：
[0043]
图1为本发明中脑梗塞核心区变化的血管形态学指标的提取与筛选方法的流程图。
[0044]
图2为本发明中3d unet网络模型的结构示意图。
[0045]
图3为本发明中脑血管的3d中心线的成像结构示意图。
[0046]
图4为本发明中脑血管的的网络化描述结构示意图。
[0047]
图5为本发明中脑血管形态学特征提取的流程图。
[0048]
图6为本发明中11个脑血管形态结构特征图；
[0049]
图7为本发明中脑血管形态学的影像组学分类预测模型的性能数值表达图。
[0050]
图8为本发明中脑血管形态学的影像组学分类预测模型的roc曲线最佳阈值点图。
[0051]
图9为本发明中脑血管形态学的影像组学分类预测模型的性能表达图。
[0052]
图10为本发明中脑血管形态学的影像组学分类预测模型的对梗塞变化趋势预测的效能曲线表达图。
[0053]
图11为本发明中受试者的梗塞体积变化趋势的示意图。
[0054]
图12为本发明中受试者的两次mri-dwi图像成像中的梗塞体积分割图像图。
具体实施方式
[0055]
为使本发明的目的、技术方案和优点更加清楚明了，下面通过附图及实施例，对本发明进行进一步详细说明。但是应该理解，此处所描述的具体实施例仅仅用以解释本发明，并不用于限制本发明的范围。此外，在以下说明中，省略了对公知结构和技术的描述，以避免不必要地混淆本发明的概念。
[0056]
参看图1-图5所示，脑梗塞核心区变化的血管形态学指标的提取与筛选方法，包括以下步骤：
[0057]
(1)ctp影像获取：先获取受试者的ctp影像结果；
[0058]
所述ctp影像的体素分辨率和大小分别为0.486*0.486*1mm3和512*512*319，为了后续梗塞诊断分析的需要，本实施例提取了160到319个最有可能发生梗塞的切片图像进行实验。所述ctp影像的体素灰度由窗口[0,600](hu)裁剪，所述ctp影像的体素灰度用于聚焦脑血管。
[0059]
具体的，在本实施例中在，将受试者111例，将111个受试者的ctp影像随机分成67:22:22分别用于训练、验证和测试。由放射科医生在itk-snap平台上手动标记，一名负责标记，另一名负责校正。同时在轴向平面、冠状平面和矢状平面上进行标记。
[0060]
(2)ctp影像血管分割：采用3d unet网络模型对所述ctp影像的脑血管进行分割，所述3d unet网络模型通过跳层连接和下采样层将不同尺度图像信息融合，实现所述脑血管分割。
[0061]
具体的，所述3d unet网络模型在训练阶段输入和输出的ctp的3d图像，进行编码和解码，编码阶段的特征图通过跳层连接传输到解码阶段进行图像特征的卷积融合，得到训练图像数据；然后抽取训练图像数据的血管和背景区域，进行模型训练；最后进行测试，得到整图分割结果。
[0062]
作为对本实施例的进一步说明，所述3d unet网络模型在训练阶段输入和输出的3d图像块尺寸均为128
×
128
×
128。编码阶段和解码阶段均包含4层卷积。在编码阶段，每个卷积层依次包含两个3
×3×
3大小的卷积和一个线性整流函数(rectified linear unit,relu)。除最底层外，每个卷积层后面跟一个尺度为2
×2×
2、每个维度步长为2的最大池化层。第一个卷积层输出特征图通道数为32，后面每层特征图通道数随层数增加依次翻倍。在解码阶段，3个反卷积层与解码阶段3个最大池化层对应部署，反卷积的卷积核大小为2
×2×
2，步长为2。每个反卷积层后面依次部署两个3
×3×
3大小卷积核的卷积操作和一个线性整流函数relu操作。除最后一个卷积层外的每个卷积层减半输入特征图通道数。最后一个卷积层使用1
×1×
1的卷积核将网络输出通道数减少为标签图像通道数。在具有相同尺寸特征图的层上，编码阶段的特征图通过跳层连接传输到解码阶段进行图像特征的卷积融合，得到ctp影像的训练图像数据。
[0063]
为了更好地训练网络，本实施例在每个线性整流函数relu前加入批归一化(batch normalization,bn)。
[0064]
作为对本实施例的进一步说明，在所述编码阶段的最底层加入了损失函数，避免模型过拟合。所述损失函数l
dc
的表达如下：
[0065][0066]
其中c∈c表示类别数，u为sofmax层输出结果，v为onehot编码后标签，i表示图像中的像素点。
[0067]
作为对本实施例的进一步说明，参看图2，分别从训练图像数据的血管和背景区域抽取128
×
128
×
128大小的图像块对。每个受试者随机提取100个图像块。所有提取的图像块随机打乱8次输入网络训练，训练过程中批尺寸设置为2，卷积核l2范数约束权重为1e-5，采用自适应矩估计(adaptive moment estimation,adam)优化器，学习率为1e-4。进行模型训练42000次。然后进行测试，在测试阶段，整个图像输入训练模型得到整图分割结果。模型使用基于python编写的张量流(tensorflow,tf)库实现，运算平台为3.00ghz的intel xeon e5-2587w(v4)cpu和nvidia titan xp gpu。
[0068]
(3)脑血管形态学特征提取：参看图3在分割所述脑血管后，进一步量化脑血管形态结构信息，先提取所述脑血管的3d中心线，提取结果中包含脑血管中心线和各节点信息；然后基于所述脑血管的3d中心线，提取脑血管形态结构特征，得到3d脑血管形态学特征。所述脑血管形态结构特征包括：血管体积，血管分支数，血管总长度，平均血管长度，平均距离因子df，平均角度soam，边结点个数，网络结点个数，聚类系数c，结构熵h和平均路径长度。最后根据所述3d脑血管形态学特征进一步提取2d脑血管形态学特征，由此得到受试者的整体脑血管形态学特征及梗塞周围区域脑血管形态学特征。
[0069]
具体的，参看图6，可得到的11个脑血管形态结构特征。
[0070]
作为对本实施例的进一步说明，所述平均距离因子df为血管总长与端点直线距离的比值，所述平均距离因子df用于衡量血管的弯曲度。
[0071]
作为对本实施例的进一步说明，所述平均角度soam用于衡量血管的扭曲度，所述平均角度soam通过四个结点间的3个向量来度量血管3d空间的扭曲度；
[0072]
作为对本实施例的进一步说明，参看图4，
[0073]
令t1＝p
k-p
k-1
,t2＝p
k+1-pk,t3＝p
k+2-p
k+1
,则脑血管在pk点的弯曲度dpk和扭度cpk可以表示为：
[0074][0075][0076]
其中，
[0077][0078]
作为对本实施例的进一步说明，所述聚类系数c和结构熵h是网络拓扑图的度量参数，所述聚类系数c计算公式如下：
[0079][0080]
其中n表示节点个数，li表示节点i的邻节点个数，ki表示i与邻节点间的边个数；
[0081]
结构熵h计算公式如下：
[0082][0083]
其中表示第i个节点的度。
[0084]
作为对本实施例的进一步说明，在所述步骤(3)中，根据所述3d脑血管形态学特征进一步提取2d脑血管形态学特征的具体过程如下：为更详细描述血管形态学特征，除直接提取3d脑血管形态学特征外，本实施例进一步提取血管2d特征。参看图5，其中图5的a表示2d血管投影，图5的b表示2d血管增强，图5的c表示2d血管分割,图5的d表示中心线提取。
[0085]
具体过程如下：首先将3d血管投影到2d轴状面；然后进行血管增强，2d阈值分割和提取血管中心线操作；最后提取2d血管的11个脑血管形态结构特征。通过提取大脑左右两侧血管结构的差异信息来进一步分析脑梗塞，将二者的差值作为对称两侧血管差异特征。
[0086]
具体的，由于脑梗塞的出现改变了大脑血管原始的对称结构，因此可以通过提取大脑左右两侧血管结构差异信息来进一步分析脑梗塞。即首先提取左右两侧血管形态学特征11个，然后将二者的差值作为对称两侧血管差异特征。最终每个受试者提取血管形态学特征44个，其中2d和3d特征各22个。通过提取大脑左右两侧血管结构的差异信息来进一步分析脑梗塞，将二者的差值作为对称两侧血管差异特征。
[0087]
(4)脑血管形态学特征筛选：在建立得到脑血管形态学的影像组学分类预测模型前，通过对受试者的两次mri-dwi图像结果进行分割梗塞核心区，获取梗死核心区的体积，将前后两次对比，得到差值作为不同梗塞变化的分组依据，将所述梗塞变化分为梗塞逆转组和梗塞扩大组。
[0088]
梗塞的变化分为梗塞扩大和梗塞逆转。梗塞扩大是指基线dwi成像未出现高信号的区域，后期出现了梗塞的组织病理学变化。梗塞逆转是指血流再灌注后，dwi高信号区域的神经元损害的组织病理学变化出现的逆转。
[0089]
每个受试者根据初查和复查的结果，以梗塞变化的结局分为梗塞逆转组和梗塞扩大组(图11)，分别对两次mri-dwi图像结果进行分割梗塞核心区。
[0090]
具体的，通过对受试者的两次mri-dwi图像结果进行分割梗塞核心区的方法如下：
如果存在梗塞，梗塞核心由医生使用开源dicom查看软件在b＝1000s/mm2 mri-dwi影像上分割，根据与周围组织和对侧半球相比的高信号，将组织确定为梗塞。根据来自整个分割dwi图像图像的周边信息，包括对称性/双边性、金属、其他致病材料的存在以及颅底附近或脑、骨、空气界面上的典型解剖位置，将类似于梗塞核心的成像伪影与真实梗塞核心区分开来。梗塞延续或梗塞内的实质血肿包括在分割中。如图12所示：一例受试者手动分割完的结果，其中a为冠状面，b为横断面，c为矢状面，分割区域显示为p。分别在冠状面，横断面和矢状面可见p的分割区域，据此计算出梗塞的体积。
[0091]
然后根据每个所述受试者的整体脑血管形态学特征及梗塞周围区域脑血管形态学特征，将梗死核心区体积扩大或逆转预测转换为机器学习二分类问题，并建立基于非参数训练的稀疏表示分类器进行分类预测，得到脑血管形态学的影像组学分类预测模型。
[0092]
作为对本实施例的进一步说明，建立基于非参数训练的稀疏表示分类器进行分类预测，得到脑血管形态学的影像组学分类预测模型的过程如下：
[0093]
首先利用稀疏编码和重构计算待测试样与训练集里每类样本间的距离，然后根据k-邻近准则判断测试样本的类别分类过程具体可表示为：
[0094][0095]
其中表示待测样本特征选择后的特征，为训练样本特征集，c为样本总类别数；γ为稀疏表示控制参数；当获得稀疏表示系数后，计算待测试样与训练集里每类样本间的距离，即各类的表示残差：
[0096][0097]
其中δc(
·
)表示选择第c类特征对应的系数。残差最小的类别即为测试样本的类别，因此最终待测样本类别为
[0098]
(5)脑血管形态学特征模型验证：所述脑血管形态学的影像组学分类预测模型建立之后，提取并筛选出受试者的脑血管形态学特征，纳入所述脑血管形态学的影像组学分类预测模型中进行预测，将预测结果与受试者7-10天的复查的mri-dwi图像结果进行比较，来评估预测的有效性和准确性。
[0099]
按照上述的方法，本发明通过提取闵行中心医院111例受试者的ctp影像。ctp影像提取40个血管特征。前期投入74例受试者训练，后期前瞻性纳入37例受试者测试。
[0100]
预测模型的性能比较如下：具体的，本发明对输入的院111例受试者的ctp影像，采用matlab软件，通过基于非参数训练的稀疏表示分类器进行分类预测，得到mri-dwi图像的影像组学分类预测模型。并使用r语言的precrece包中evalmod函数为指定的模型计算基本评价指标。包括准确率acc(accuracy)、特异性spe(specificity)、敏感性sen(sensitivity)、错误率err(error)、精确率p(precision)、matthews相关系数mcc(matthews correlation coefficient)和f评分(f-score)，用归一化分级(normalized ranks)值之间的散点图作为选择截断值的参考。
[0101]
采用受试者工作特征(receiver operating characteristic curve，roc)曲线来评估四个预测模型的性能，用以比较模型的效能。同时显示6个指标：roc曲线下面积auc(the area under the corresponding curve)、acc、sen、spe、阳性预测值ppv、阴性预测值
npv。它们的计算公式如下：
[0102]
准确率
[0103]
特异性
[0104][0105]
敏感性
[0106]
阳性预测值
[0107]
阴性预测值
[0108]
其中“真正”tp(true positive)表示正确识别梗塞扩大的受试者，“真负”tn(true negative)表示正确识别梗塞逆转的受试者，“假正”fp(false positive)表示将实际上的梗塞逆转被识别为梗塞扩大的受试者，“假负”fn(false negative)表示实际上的梗塞扩大被识别为梗塞逆转的受试者。
[0109]
roc曲线中：纵轴(y轴)是true positive rate(tpr，等价于sensitivity)，即正例样本中被正确判定为正例的样本数比例tpr；横轴(x轴)是false positive rate(fpr，等价于1-specificity)，即负例样本中被错误判定为正例样本的比例fpr，具体计算公式如下：
[0110]
正例样本中被正确判定为正例的样本数比例
[0111][0112]
负例样本中被错误判定为正例样本的比例
[0113][0114]
错误率err定义为分类错误的样本数占样本总数的比例；精确率p(或称查准率)指判断为正例的样本中有多大比例是真的正例；召回率r(recall)指正例样本中有多大比例的正例被发现(判定为正例)，该指标也称为tpr或sen。具体计算公式如下：
[0115][0116]
精确率
[0117]
召回率
[0118]
评分f是精确率p和召回率r的调和平均值；mcc是评估模型的指标，在正例和负例的不同划分对于最终结果没有影响。具体计算公式如下：
[0119]
评分
[0120][0121]
评估模型的指标
[0122]
本发明方法的脑血管形态学的影像组学分类预测模型验证的具体结果参看图7、图8、图9和图10。在标化等级数值等于0.594时，评价指标处于最优状态，acc为0.756达到最高值，err为0.243是最低值，精确率为0.727相对高值，mcc为0.518是最高值，f评分为0.780为最高值。
[0123]
本发明通过上述实施例发现，具体如下。
[0124]
1.根据模型提取血管特征因子，提示三维整体血管的数量，受试者与健康者血管对比数量和血管的分叉数量对比，都是影响梗塞变化的因素。
[0125]
2.在临床上证实了细微的小血管减少确实可以预测梗塞的变化趋势。
[0126]
3.基于血管形态的预后预测虽然取得了一定的效果，整体预测效能较高，后期在临床上对脑血管量化的血管形态学特征指标进一步增加，可提高整体预测率。
[0127]
以上仅为本发明的具体实施例，但本发明的技术特征并不局限于此。任何以本发明为基础，为解决基本相同的技术问题，实现基本相同的技术效果，所作出地简单变化、等同替换或者修饰等，皆涵盖于本发明的保护范围之中。