一种基于核磁共振图像的脑血管硬化检测方法

1.本发明涉及医学影像处理技术领域，具体涉及一种基于核磁共振图像的脑血管硬化检测方法。


背景技术：

2.脑卒中已成为世界范围内第二大致残或致命的疾病，其主要成因为脑动脉瘤、动静脉畸形等严重的脑血管疾病。根据世界卫生组织相关统计，全世界每年有1500万人中风，其中500万人死亡，另有500万人终身残疾。由此可见，针对脑血管疾病的早期诊断、安全治疗以及术后康复检查等方面的工作已刻不容缓。
3.近些年，随着医学影像技术与计算机辅助诊断技术的飞速发展，磁共振血管造影(magnetic resonance angiography，mra)等技术在脑血管疾病的临床诊断和治疗中被普遍采用。在临床诊断中，医生从图像中获取的脑血管中心线信息，可以在脑血管疾病的治疗上起到良好的协助作用。因此面向mra图像的脑血管中心线提取技术成为医学图像处理领域的核心问题之一。由于人们常常是对心脏冠脉血管和视网膜血管进行研究，然而对于脑血管的研究少之又少。


技术实现要素：

4.本发明的目的在于提供一种基于核磁共振图像的脑血管硬化检测方法，通过分析处理三维时间飞跃磁共振血管图像，为医务人员诊断脑血管硬化提供参考。
5.为了实现上述目的，本发明采用的技术方案如下：
6.一种基于核磁共振图像的脑血管硬化检测方法，包括以下步骤：
7.s1：获取检测对象的脑血管图像，并对脑血管图像进行预处理；
8.s2：基于fmm-mrf脑血管分割模型，提取脑血管特征分割图；
9.s3：获取脑血管中心线；
10.s4：根据血管中心线计算得到血管直径分布图；
11.s5：计算血管狭窄程度，并根据脑血管狭窄测量标准评估脑血管硬化程度。
12.进一步的，所述s1中脑血管图像通过三维时间飞跃磁共振血管成像(3d-tof-mra)技术获取。
13.进一步的，所述s1中预处理具体包括以下步骤：
14.s11：对三维时间飞跃磁共振血管图像进行偏置场修正；
15.s12：去除脑血管图像中的颅骨结构；
16.s13：将去除颅骨的血管图像输入基于hessian矩阵特征值的多尺度滤波器，得到增强后的血管特征图。
17.进一步的，所述s2具体包括以下子步骤：
18.s21：采用k-means聚类方法获取有限混合模型的初始参数；
19.s22：采用最大期望算法对有限混合模型的初始参数进行更新优化，进而建立有限
混合模型；
20.s23：将脑血管图像输入有限混合模型得到血管分割结果；
21.s24：血管分割结果作为马尔可夫随机场统计过程中的初始标记场，结合血管特征图进一步优化分割结果。
22.进一步的，所述s3具体包括以下子步骤：
23.s31：通过细化算法获取脑血管骨架线；
24.s32：提取脑血管骨架线血管分支点集；
25.s33：根据局部边缘距离场进行多分支点加权平均，将多个血管分支点合并为一个分支点，进而优化脑血管骨架线；
26.s34：通过优化后的脑血管骨架线转化为脑血管中心线。
27.进一步的，所述s4具体包括以下子步骤：
28.s41：提取脑血管中心线的每一点；
29.s42：采用基于mom评价模型的跟踪算法测量每一点的直径值。
30.进一步的，所述s42具体包括以下子步骤：
31.s421：计算感兴趣区域的每一中心点的血管方向；
32.s422：获取垂直于每一中心点的血管方向的左右边缘点；
33.s423：计算左右边缘点距中心点的距离，即中心点所在血管半径；
34.s424：根据血管半径构建脑血管直径分布图。
35.进一步的，所述s5具体包括：
36.s51：根据脑血管直径分布图获取颅内动脉管径和残存血管最小管径；
37.s52：按照下列公式计算狭窄程度：
38.脑血管狭窄程度＝(颅内动脉管径-残存血管最小管径)/颅内动脉管径
×
100％；
39.s53：参照北美症状性动脉内膜切除术研究协助组制定的脑血管狭窄测量标准，评估脑血管硬化程度。
40.本发明的有益效果：
41.(1)本发明通过对三维时间飞跃磁共振血管图像进行预处理，并且采用多尺度滤波器增强血管目标，能够能加精确的提取目标血管的中心线，同时抑制非血管结构的干扰，有效地去除部分噪声干扰；
42.(2)本发明通过颅骨剔除技术和脑血管特征分割来精简感兴趣区域；由于血管的颅内容积小，传统em算法对fmm参数估计过程是有偏的，由此提出基于血管知识的em算法；将多尺度血管增强结果与fmm中血管容积估计权值相结合，进一步得到具有概率意义的血管特征图，并将其嵌入到mrf过程，构建了双能量约束的mrf血管优化过程，在分割精度和鲁棒性上相对传统统计模型明显提高。
43.(3)本发明针对中心线提取与介入路径规划：在细化算法得到骨架线的基础上，结合边缘距离场解决了分支点不唯一问题，得到精确的脑血管中心线；
44.(4)本发明通过图像处理得到的脑血管中心线，进一步计算血管直径，为医务人员诊断脑血管硬化提供可靠数据支持。
附图说明
45.图1为本发明的整体方案流程图；
46.图2为本发明的s1步骤流程图；
47.图3为本发明的s2步骤流程图；
48.图4为本发明的s3步骤流程图；
49.图5为本发明的s4步骤流程图；
50.图6为本发明的s5步骤流程图。
具体实施方式
51.以下将参照附图和优选实施例来说明本发明的实施方式，本领域技术人员可由本说明书中所揭露的内容轻易地了解本发明的其他优点与功效。本发明还可以通过另外不同的具体实施方式加以实施或应用，本说明书中的各项细节也可以基于不同观点与应用，在没有背离本发明的精神下进行各种修饰或改变。应当理解，优选实施例仅为了说明本发明，而不是为了限制本发明的保护范围。
52.需要说明的是，以下实施例中所提供的图示仅以示意方式说明本发明的基本构想，遂图式中仅显示与本发明中有关的组件而非按照实际实施时的组件数目、形状及尺寸绘制，其实际实施时各组件的型态、数量及比例可为一种随意的改变，且其组件布局型态也可能更为复杂。
53.如图1所示，一种基于核磁共振图像的脑血管硬化检测方法，包括以下步骤：
54.s1：获取检测对象的脑血管图像，并对脑血管图像进行预处理；
55.s2：基于fmm-mrf脑血管分割模型，提取脑血管特征分割图；
56.s3：获取脑血管中心线；
57.s4：根据血管中心线计算得到血管直径分布图；
58.s5：计算血管狭窄程度，并根据脑血管狭窄测量标准评估脑血管硬化程度。
59.如图2所示，所述s1中脑血管图像通过三维时间飞跃磁共振血管成像技术获取。
60.s11：对三维时间飞跃磁共振血管图像进行偏置场修正；本实施例利用simpleitk中的n4biasfieldcorrectionimagefilter方法实现了tof-mra图像中的偏置场的修正，采用otsu阈值分割算法获取一个掩膜；用于控制偏置场修正范围，在算法执行中，算法迭代16次，可以得到稳定的偏置场修正结果。
61.s12：去除脑血管图像中的颅骨结构；本实施例中，对于3d-tof-mra图像，不仅包含脑组织以及分布在其表层间隙的脑血管，也含有头骨皮层、眼晴以及背景空气，所以为了排除无关信息的干扰，可采用fsl中的bet工具对图像进行去除颅骨。
62.s13：将去除颅骨的血管图像输入基于hessian矩阵特征值的多尺度滤波器，得到增强后的血管特征图。在三维图像中，对于脑血管结构，本实施例将脑血管结构响应函数定义如下:
63.64.其中，ra是一个特征值几何比，用来区分盘形结构和管状结构组织。表示如下:
[0065][0066]
rb也是一个几何比，用来区分球形结构和管状结构组织。表示如下:
[0067][0068]
为了区分背景像素，引进s，
[0069][0070]
其中a、b和c分别为ra、rb和s的参数，用于调节血管结构相似函数对图像敏感度。在通常情况下，a和b取固定值0.5就可以获得很好的效果，而c的值由图像的灰度范围决定，本实施例取图像的hessian矩阵最大范数值的一半即:
[0071][0072]
当三维图像中的像素点不是血管结构而是噪声时，该点对应hessian矩阵的特征值λ1、λ2和λ3都为一个接近零的值；此时，血管响应函数v0(s)的值会很小。而当图像中的像素点属于盘形结构或者球形结构时，ra和rb分别约等于零；血管响应函数v0(s)的值也会接近零。只有当图像中的像素点位于管状结构中也就是血管中时，血管响应函数v0(s)才会取得一个较大值。由此可知，通过此函数可以有效地增强图像中的血管结构。
[0073]
如图3所示，所述s2具体包括以下子步骤：
[0074]
s21：采用k-means聚类方法获取有限混合模型的初始参数，设置聚类个数为3，初始聚类中心为峰值的四分之一点、峰值点、以及峰值点的两倍点，并根据聚类后各类的均值和方差作为fmm的初始参数；
[0075]
s22：采用最大期望算法对有限混合模型的初始参数进行更新优化，进而建立有限混合模型；
[0076]
s23：将脑血管图像输入有限混合模型得到血管分割结果；
[0077]
s24：血管分割结果作为马尔可夫随机场统计过程中的初始标记场，结合血管特征图进一步优化分割结果；本实施例中，fmm的分割结果作为mrf过程中的初始标记场，并结合血管特征图构建了一个全新的六邻域势能函数，其定义如下:
[0078][0079]
其中γ是体素的索引值，γ*∈η
γ
,是六领域内的体素的索引值，x
γ
，x
γ
*是标记场中体素对应的标记值。标记值来自集合l＝{l
ν
，lb}，l
ν
为血管标记值，lb为背景标记值。vf
γ
是血管特征图中索引值为γ的点的概率值。其中比例系数满足α1+α2＝1，且参数分别设为0.5。此外函数和定义为:
[0080][0081][0082]
如图4所示，所述s3具体包括以下子步骤：
[0083]
s31：通过细化算法获取脑血管骨架线；
[0084]
s32：提取脑血管骨架线血管分支点集；本实施例中将脑血管骨架线上的点集再次划分为三个子集:
[0085]
血管端点集s1：在其26领域内包含的骨架线点仅为2个，即，该点与周围点的连接度为1；
[0086]
血管路径点集s2：在其26领域内包含的骨架线点仅为3个，即，该点与周围点的连接度为2，骨架线上绝大多数的点属于该点集；
[0087]
血管分支点集s3：在其26领域内包含的骨架线点大于或等于4个，即，该点与周围点的连接度大于2，该点在血管系统中起到了血管分流或合流作用，也是介入路径中重点考察点。
[0088]
s33：根据局部边缘距离场进行多分支点加权平均，将多个血管分支点合并为一个分支点，进而优化脑血管骨架线；对于血管目标，采用欧儿里得距离场。对于血管分支处的n个可能分支点进行加权合并为一个分支点:
[0089][0090]
其中d
pi
为欧式距离场中点pi的值，其中z是正则项:
[0091][0092]
s34：通过优化后的脑血管骨架线转化为脑血管中心线。本实施例中，得到优化后的骨架线后，还需要对骨架线进一步整理。将完整的血管网络定义为vessels＝{(p
start
，p
end
，l
path
)i；i＝1,2，...n}，其中(p
start
，p
end
，l
path
)i表示一段血管，且该段血管路径上无血管分支，其中p
start
，p
end
∈s1∪s3来自血管端点集或血管分支点集，队列l
path
∈s2所有点来自血管路径点集，具体转化过程如下：
[0093]
(1)将血管端点集s1或血管分支点集s3合并定义k＝(s1∪s3，degree)，其中degree记录了各个点的连接度。
[0094]
(2)初始化起始点队列q；初始化队列l；骨架线体素标记空间skel，其中背景标记为1，骨架线点为0；i为血管段识别号，令i＝0。
[0095]
(3)若队列q和集合k均为空，则结束转化过程。否则，从集合k中任取一点p，若点p的degree＝m，则向队列q中加入m个p。集合k中删除点p。
[0096]
(4)将队列q中队头元素点q
top
出队列，令i＝i+1，(p
start
，p
end
，l
path
)i中p
start
＝q
top
，以q
top
开始，从其26邻域内选择一未访问点l
p
对其访问并记录skel(l
p
)＝1，点l
p
进入队列l
path
，依次逐步向下遍历。
[0097]
(5)当l
p
∈k，令p
end
＝l
p
，停止遍历，将集合k中点l
p
的degree＝degree-1，若点l
p
的degree不为0，则将相应数量的点l
p
加入队列q。
[0098]
(6)若队列q非空，则重复步骤(4)。
[0099]
(7)若队列q空，则重复步骤(3)。
[0100]
经过以上步骤，从而得到了每根血管的序列排布，以及每根血管之间的连接关系。
[0101]
对于一些细长的血管，由于局部体素太少(可能才2-4个体素)，无法真正的得到血管中心，为此利用原始数据的血管增强图对每根血管(p
start
，p
end
，l
path
)i中l
path
序列点位置进行微调。计算每根血管的骨架线点l
path
(k)的切线方向，并获得垂直于切线方向的分割掩
模与血管增强图的局部横截面s。设横断面中所有血管体素的血管增强集合为ω，那么微调后的中心点l
path
(k)
*
为:
[0102][0103]
其中s(r)为局部横截面内体素的坐标，z是一个正则项:
[0104][0105]
对调整后的各段脑血管进行插值平滑处理，得到最终的脑血管中心线。
[0106]
如图5所示，所述s4具体包括以下子步骤：
[0107]
s41：提取脑血管中心线的每一点；
[0108]
s42：采用基于mom评价模型的跟踪算法测量每一点的直径值。
[0109]
所述s42具体包括以下子步骤：
[0110]
s421：计算感兴趣区域的每一中心点的血管方向；
[0111]
s422：获取垂直于每一中心点的血管方向的左右边缘点；
[0112]
s423：计算左右边缘点距中心点的距离，即中心点所在血管半径；
[0113]
s424：根据血管半径构建脑血管直径分布图。
[0114]
本实施例中，基于mom评价模型的跟踪过程如下:
[0115]
1)取中心线上的起始点作为当前中心点p，计算点p的血管方向θ
p
。
[0116]
2)在垂直于血管的方向θ上，以中心点p为中心、以距离d同时向左右两边进行搜索，对应的左右边缘点分别是p1和p2。以中心点左边的搜索为例,此时对应的边缘点为p1，当d＝d0。时,在距离d
0-0.5和距离d0+0.5的像素点上分别沿垂直于θ方向上向两边同时取m个像素点(在实际计算中m＝2)，连同搜索线上的像素点共有2m+1个像素点。假设分别为f1(t-m)，
…
，f1(t-1)，f1(t)，f1(t+1)，
…
，f1(t+m)和g1(t-m)，
…
，g1(t-1)，g1(t)，g1(t+1)，
…
，g1(t+m)。
[0117]
对于右边的搜索同样可以得到：f2(t-m)，
…
，f2(t-1)，f2(t)，f2(t+1)，
…
，f2(t+m)和g2(t-m)，
…
，g2(t-1)，g2(t)，g2(t+1)，
…
，g2(t+m)。
[0118]
3)设
[0119][0120]
其中宽度变量d是mom的参数变量。随着宽度变量d的变化mom值的大小也会变化。由于在血管边缘处灰度值会发生跳跃性的变化(梯度值局部最大)，因此对于计算出各个mom值，最大的mom值所对应的d值就是当前中心点处血管的半径值，2d即为直径值，此时称为达到最优。在实际计算中，为了更加精确，将最优化时左右两边缘点p1和p2之间的距离作为当前中心点处的直径值。
[0121]
4)选取中心线上点p的下一个点作为当前中心点，再转到步骤1)，继续进行跟踪，直到到达中心线的终点为止。
[0122]
如图6所示，所述s5具体包括：
[0123]
s51：根据脑血管直径分布图获取颅内动脉管径和残存血管最小管径；
[0124]
s52：按照下列公式计算狭窄程度：
[0125]
脑血管狭窄程度＝(颅内动脉管径-残存血管最小管径)/颅内动脉管径
×
100％；
[0126]
s53：参照北美症状性动脉内膜切除术研究协助组制定的脑血管狭窄测量标准，评估脑血管硬化程度。
[0127]
以上实施例仅是为充分说明本发明而所举的较佳的实施例，本发明的保护范围不限于此。本技术领域的技术人员在本发明基础上所作的等同替代或变换，均在本发明的保护范围之内。

技术特征：
1.一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，包括以下步骤：s1：获取检测对象的脑血管图像，并对脑血管图像进行预处理；s2：基于fmm-mrf脑血管分割模型，提取脑血管特征分割图；s3：获取脑血管中心线；s4：根据血管中心线计算得到血管直径分布图；s5：计算血管狭窄程度，并根据脑血管狭窄测量标准评估脑血管硬化程度。2.根据权利要求1所述的一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，所述s1中脑血管图像通过三维时间飞跃磁共振血管成像技术获取。3.根据权利要求2所述的一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，所述s1中预处理具体包括以下步骤：s11：对三维时间飞跃磁共振血管图像进行偏置场修正；s12：去除脑血管图像中的颅骨结构；s13：将去除颅骨的血管图像输入基于hessian矩阵特征值的多尺度滤波器，得到增强后的血管特征图。4.根据权利要求3所述的一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，所述s2具体包括以下子步骤：s21：采用k-means聚类方法获取有限混合模型的初始参数；s22：采用最大期望算法对有限混合模型的初始参数进行更新优化，进而建立有限混合模型；s23：将脑血管图像输入有限混合模型得到血管分割结果；s24：血管分割结果作为马尔可夫随机场统计过程中的初始标记场，结合血管特征图进一步优化分割结果。5.根据权利要求1所述的一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，所述s3具体包括以下子步骤：s31：通过细化算法获取脑血管骨架线；s32：提取脑血管骨架线血管分支点集；s33：根据局部边缘距离场进行多分支点加权平均，将多个血管分支点合并为一个分支点，进而优化脑血管骨架线；s34：通过优化后的脑血管骨架线转化为脑血管中心线。6.根据权利要求5所述的一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，所述s4具体包括以下子步骤：s41：提取脑血管中心线的每一点；s42：采用基于mom评价模型的跟踪算法测量每一点的直径值。7.根据权利要求6所述的一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，所述s42具体包括以下子步骤：s421：计算感兴趣区域的每一中心点的血管方向；s422：获取垂直于每一中心点的血管方向的左右边缘点；s423：计算左右边缘点距中心点的距离，即中心点所在血管半径；s424：根据血管半径构建脑血管直径分布图。
8.根据权利要求1所述的一种基于核磁共振图像的脑血管硬化检测方法，其特征在于，所述s5具体包括：s51：根据脑血管直径分布图获取颅内动脉管径和残存血管最小管径；s52：按照下列公式计算狭窄程度：脑血管狭窄程度=(颅内动脉管径-残存血管最小管径)/颅内动脉管径
×
100％；s53：参照北美症状性动脉内膜切除术研究协助组制定的脑血管狭窄测量标准，评估脑血管硬化程度。

技术总结
本发明提供一种基于核磁共振图像的脑血管硬化检测方法，涉及医学影像处理技术领域，具体包括以下步骤：S1：获取检测对象的脑血管图像，并对脑血管图像进行预处理；S2：基于FMM-MRF脑血管分割模型，提取脑血管特征分割图；S3：获取脑血管中心线；S4：根据血管中心线计算得到血管直径分布图；S5：计算血管狭窄程度，并根据脑血管狭窄测量标准评估脑血管硬化程度。本发明通过图像处理得到的脑血管中心线，进一步计算血管直径，为医务人员诊断脑血管硬化提供可靠数据支持。供可靠数据支持。供可靠数据支持。


技术研发人员：张赣 刘旺 王清琛 李梦雅
受保护的技术使用者：四川省中医药科学院中医研究所
技术研发日：2023.07.21
技术公布日：2023/10/11