脑灌注影像的处理方法、系统、设备与存储介质

1.本发明涉及医学影像处理技术领域，尤其涉及一种脑灌注影像的处理方法、系统、设备与存储介质。


背景技术：

2.利用ct和mri扫描的灌注影像是一段时间内的动态序列，需要对其进行一系列处理才能够得到相关参数进行分析。一般来说，通过对选定感兴趣层面进行连续动态扫描，获得所选层面的每一像素的时间-浓度曲线，并通过数学模型处理得到：脑血供容积(cerebral blood volume，cbv)、脑血供速率(cerebral blood flow，cbf)、平均通过时间(mean transit time，mtt)、残留函数达峰时间(time to max enhancement，tmax)等血流动力学参数和灌注图像表现。
3.目前处理方法主要分为人工处理和计算机辅助处理。人工处理方法要求阅片人具有较高的理论分析水平，并且流程复杂、操作耗时，存在主观差异，可操作性差。因此，计算机辅助分析具有较高的现实意义。然而，现有计算机辅助分析方法大多无法实现全自动处理，需要人工选定感兴趣区域以选择动脉输入函数；除此之外，在计算灌注参数时无法识别如脑室等无效区域，从而导致定量计算存在较大误差。


技术实现要素：

4.本发明的目的是提供一种脑灌注影像的处理方法、系统、设备与存储介质，能够针对脑灌注影像，自动进行预处理、函数选择和参数计算，给出定量分析结果。
5.本发明的目的是通过以下技术方案实现的：
6.一种脑灌注影像的处理方法，包括：
7.获取原始脑灌注影像数据，并识别数据类型；
8.对所述原始脑灌注影像数据进行预处理，并结合数据类型进行时域插值；
9.通过自动检测算法从时域插值处理后的脑灌注影像数据中自动选择动脉输入函数；
10.基于时域插值处理后的脑灌注影像数据和动脉输入函数，计算灌注参数；
11.从预处理后的原始脑灌注影像数据提取出颅内脑组织区域掩膜，并结合灌注参数对所述预处理后的原始脑灌注影像数据进行后处理和定量统计。
12.一种脑灌注影像的处理系统，用于实现前述的方法，该系统包括：
13.数据获取与类型识别模块，用于获取原始脑灌注影像数据，并识别数据类型；
14.数据预处理与时域插值处理模块，用于对所述原始脑灌注影像数据进行预处理，并结合数据类型进行时域插值；
15.动脉输入函数自动选择模块，用于通过自动检测算法从时域插值处理后的脑灌注影像数据中自动选择动脉输入函数；
16.灌注参数计算模块，用于基于时域插值处理后的脑灌注影像数据和动脉输入函
数，计算灌注参数；
17.影像后处理和定量统计模块，用于从预处理后的原始脑灌注影像数据提取出颅内脑组织区域掩膜，并结合灌注参数对所述预处理后的原始脑灌注影像数据进行后处理和定量统计。
18.一种处理设备，包括：一个或多个处理器；存储器，用于存储一个或多个程序；
19.其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现前述的方法。
20.一种可读存储介质，存储有计算机程序，其特征在于，当计算机程序被处理器执行时实现前述的方法。
21.由上述本发明提供的技术方案可以看出，适用于ct及mri灌注影像数据，能够自动进行预处理，选定动脉输入函数，并自动去除无效区域，从而提供了更准确、便捷、快速的灌注参数定性及定量分析，对提升灌注影像阅片质量具有积极影响。
附图说明
22.为了更清楚地说明本发明实施例的技术方案，下面将对实施例描述中所需要使用的附图作简单地介绍，显而易见地，下面描述中的附图仅仅是本发明的一些实施例，对于本领域的普通技术人员来讲，在不付出创造性劳动的前提下，还可以根据这些附图获得其他附图。
23.图1为本发明实施例提供的一种脑灌注影像的处理方法的流程图；
24.图2为本发明实施例提供的运动校准结果的示意图；
25.图3为本发明实施例提供的静脉输出函数的定位结果示意图；
26.图4为本发明实施例提供的动脉输入函数定位结果示意图；
27.图5为本发明实施例提供的动脉输入函数、静脉输出函数和全局平均的时间-浓度曲线示意图；
28.图6为本发明实施例提供的一种脑灌注影像的处理系统的示意图；
29.图7为本发明实施例提供的一种处理设备的示意图。
具体实施方式
30.下面结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明的保护范围。
31.首先对本文中可能使用的术语进行如下说明：
32.术语“包括”、“包含”、“含有”、“具有”或其它类似语义的描述，应被解释为非排它性的包括。例如：包括某技术特征要素(如原料、组分、成分、载体、剂型、材料、尺寸、零件、部件、机构、装置、步骤、工序、方法、反应条件、加工条件、参数、算法、信号、数据、产品或制品等)，应被解释为不仅包括明确列出的某技术特征要素，还可以包括未明确列出的本领域公知的其它技术特征要素。
33.下面对本发明所提供的一种脑灌注影像的处理方法进行详细描述。本发明实施例
中未作详细描述的内容属于本领域专业技术人员公知的现有技术。本发明实施例中未注明具体条件者，按照本领域常规条件或制造商建议的条件进行。本发明实施例中所用试剂未注明生产厂商者，均为可以通过市售购买获得的常规产品。
34.如图1所示，一种脑灌注影像的处理方法，主要包括如下步骤：
35.步骤1、获取原始脑灌注影像数据，并识别数据类型。
36.本步骤中，获取原始脑灌注影像数据后，可以根据相关的序列描述，判别数据类型为轴向扫描数据或是容积螺旋扫描数据。
37.此外，还可以根据标签数据读取相应用户的信息、序列长度、拍摄间隔、像素空间体积、影像尺寸等信息。其中，此处的用户信息用于软件的数据库存储查询以及前端显示；序列长度、影像尺寸设定了后续算法中的变量大小；拍摄间隔用于进行时域插值；像素空间体积被直接用于计算对应的体积。
38.步骤2、对所述原始脑灌注影像数据进行预处理，并结合数据类型进行时域插值。
39.1)数据预处理。
40.本发明实施例中，数据预处理主要是针对所述原始脑灌注影像数据进行噪声去除和运动校准，减小拍摄过程中因拍摄设备和病人运动导致的误差，随后，对影像进行颅骨去除处理；示例性的，可通过下述方式实现：
41.a)对所述原始脑灌注影像数据中的最大值和最小值进行限制，并进行双边滤波。
42.b)对于双边滤波结果，基于像素灰度互信息，使用梯度下降法进行运动校准，变换类型设置为三维欧拉变换；如图2所示，展示了运动校准结果的示例。
43.c)对于运动校准结果，基于像素灰度，使用形态学操作提取并去除颅骨区域。
44.2)时域插值处理。
45.根据拍摄扫描方法，可以将脑灌注影像的类型分为轴向扫描数据和容积螺旋扫描数据。其中轴向扫描数据每层影像的动态扫描间隔相同，而螺旋扫描数据的扫描间隔不同。因此，为了统一对两种数据进行参数计算，需要进行时域插值以统一采样间隔，而插值可以通过多种方式实现。可选的，本发明使用了二次多项式插值，并进行高斯平滑。此外，由于轴向扫描数据是等间隔扫描，还可以使用fft通过频域进行插值。
46.本领域技术人员可以理解，虽然两种类型都可进行时域插值处理，但是由于轴向扫描数据与螺旋扫描数据的拍摄顺序不一样，因此，需要结合数据类型进行时域插值处理。例如，以层的位置为标号，轴向扫描数据按123123重复，而螺旋扫描数据按123321重复。
47.步骤3、通过自动检测算法从时域插值处理后的脑灌注影像数据中自动选择动脉输入函数。
48.本步骤中，通过自动检测算法自动选择静脉输出函数(venousoutputfunction，vof)，再结合选择的静脉输出函数选择动脉输入函数(arterialinputfunction，aif)。
49.动脉输入函数和静脉输出函数的选取对于灌注参数的计算具有重要影响，相比于现有的半自动方法需要操作者手动选取感兴趣区域，本方法设计的自动检测算法能够实现全自动aif及vof选取。
50.下面针对通过自动检测算法自动选择静脉输出函数与动脉输入函数的优选方式进行介绍。
51.1)通过自动检测算法自动选择静脉输出函数。
52.a)对每个像素的时间-浓度曲线求取积分，获得每个像素的总流量。
53.b)对于每层图像，选取总流量最大的前n个像素，并对选取的n个像素的时间-浓度曲线求取平均，获得局部静脉输出函数。
54.c)对每层图像的局部静脉输出函数进行伽马函数拟合，计算相关系数：
55.r(x,y)＝cov(x,y)/var[x]var[y]
[0056]
其中x为局部静脉输出函数，y为拟合得到的伽马函数，r(x,y)为相关系数，cov(x,y)为协方差，var[
·
]为方差。
[0057]
d)遍历所有局部静脉输出函数，选出相关系数大于n1的候选局部静脉输出函数；计算每一候选局部静脉输出函数的总流量、达峰时间和半峰宽度，寻找达峰时间晚于全局平均达峰时间，且总流量与半峰宽度比值最大的候选局部静脉输出函数，作为全局静脉输出函数；其中，n1为设定值。
[0058]
e)输出全局静脉输出函数，并在对应的位置绘制定位结果。
[0059]
优选的，在寻找全局静脉输出函数时，为了减少计算量并避免偶然误差，还引入了位置先验知识。具体地，观察到正常灌注影像的拍摄具有相同的头部朝向，而函数所在静脉通常位于后脑，因此仅取影像的下半部分进行计算，不考虑影像的上半部分。而步骤d)中的全局平均达峰时间则是基于完整灌注影像(即时域插值处理后的脑灌注影像数据)的颅内区域取平均时间-浓度曲线后计算得到的，用于反映全局的流量情况。
[0060]
2)通过自动检测算法自动选择动脉输入函数。
[0061]
a)计算所述全局静脉输出函数的总流量、峰值流量和半峰宽度。
[0062]
b)计算每个像素的总流量和达峰时间。
[0063]
c)对于每层图像，筛选总流量大于全局静脉输出函数的总流量n2％，且达峰时间至少早于全局平均达峰时间n3秒的像素；对于筛选出的像素，根据计算出的峰值流量、达峰时间和半峰宽度计算相应得分；其中，n2与n3均为设定值。
[0064]
d)对于每层图像，筛选得分大于最大得分n4％的像素，计算平均时间-浓度曲线，进行伽马函数拟合得到拟合结果和相关系数；其中，n4为设定值。
[0065]
e)遍历所有拟合结果，筛选相关系数大于n1的拟合结果，取拟合结果中峰值流量与达峰时间比值最大的位置作为全局动脉输入函数。
[0066]
f)输出全局动脉输入函数，在对应位置绘制定位结果。
[0067]
示例性的，可以设置：n＝5，n1＝0.9，n2＝20，n3＝2，n4＝90。需要说明的是，对于各参数具体数值大小本发明不做限定，实际应用中，可以由用户根据实际情况或者经验自行设定。
[0068]
目前，动脉输入函数比较公认的特性为：具有较大的总流量；达峰时间早、峰值流量高；流量上升、下降迅速。因此，步骤c)中的得分需要能够反应这些特性。可选的，本发明实施中，计算的得分为峰值流量的平方与达峰时间和半峰宽度乘积的比值，表示为：
[0069][0070]
由于动脉输入函数和静脉输出函数的形状类似伽马函数，同时为了减少噪声和回流影响，对时间-浓度曲线进行伽马函数拟合是判断像素是否可以作为目标点的重要依据。而直接进行伽马函数拟合涉及到非线性优化问题，耗时长且计算困难，不利于影像的快速
处理。本发明使用了简化的伽马函数拟合算法，对伽马函数和数据进行变换，从而将其转换为线性优化问题，大幅提升计算效率。
[0071]
如图3所示，展示了静脉输出函数的定位结果示例图，图中用数字“1”标记了相应位置；如图4所示，展示了动脉输入函数定位结果示例图，图中用数字“2”标记了相应位置；图3与图4中l表示location字段，用于指示影像所处层的位置。
[0072]
如图5所示，展示了动脉输入函数、静脉输出函数和全局平均的时间-浓度曲线示例图。
[0073]
步骤4、基于时域插值处理后的脑灌注影像数据和动脉输入函数，计算灌注参数。
[0074]
本步骤的优选实施方式如下：
[0075]
1)依据造影剂稀释理论(indicator-dilution theory)，计算每个像素的残留函数k(t)。
[0076]
每个像素的时间-浓度曲线是由动脉输入函数和该像素的残留函数卷积得到的：
[0077]cvoi
(t)＝(c
art
*k)(t)
[0078]
其中，c
art
(t)为动脉输入函数，c
voi
(t)为像素的时间-浓度曲线，k(t)像素的残留函数。
[0079]
目前，对于残留函数的计算通常使用基于反卷积的方法，根据循环矩阵的构造方式可以进一步细分为多种方法。可选地，本发明例中使用了奇异值分解(svd)方法。
[0080]
2)基于残留函数，计算每个像素的灌注参数，表示为：
[0081]
cbf＝maxk(t),
[0082][0083][0084]
tmax＝argmaxk(t)
[0085]
其中，cbf为脑血供速率，mtt为平均通过时间，cbv为脑血供容积，tmax为残留函数达峰时间，t为时刻。
[0086]
步骤5、从预处理后的原始脑灌注影像数据提取出颅内脑组织区域掩膜，并结合灌注参数对所述预处理后的原始脑灌注影像数据进行后处理和定量统计。
[0087]
本步骤优选实施方式如下：
[0088]
1)对所述预处理后的原始脑灌注影像数据进行灰度分析，获得颅内脑组织区域掩膜。
[0089]
预处理后的原始脑灌注影像数据即为前述步骤2中数据预处理后得到的影像数据，也就是对原始脑灌注影像数据进行噪声去除和运动校准、以及颅骨去除处理后的影像数据。
[0090]
此处的目的是去除颅内区域中的血管等无效区域，这些区域不应在分析的考虑范围内。这些区域的时间-浓度曲线通常表现出与脑组织曲线较大的差异，因此可以通过灰度分析进行筛选。可选地，本发明例中计算了每个像素的峰值流量，并对峰值流量进行了双阈值分割，将满足条件的像素视为脑组织，从而得到颅内脑组织掩膜。
[0091]
2)基于灌注参数绘制灌注参数伪彩图。
[0092]
本发明实施例中，根据之前计算的灌注参数，进行伪彩图绘制，便于定性分析。可选地，本发明例中对灌注参数计算结果进行了高斯滤波，以减小拍摄和计算噪声的影响。
[0093]
3)基于预设阈值对灌注参数伪彩图进行划分，统计满足阈值条件的区域，并基于颅内脑组织区域掩膜去除统计区域中的无效区域。
[0094]
本发明实施例中，基于预设阈值对灌注参数伪彩图进行划分时，阈值大小可根据灌注参数的大小以及实际需求或者经验进行设定，本发明不对具体阈值大小进行限定。
[0095]
示例性的，可以选择了tmax》4s、6s、8s、10s，以及cbv《全局中位数30％作为阈值。
[0096]
如前所述，在灌注参数计算时没有对血管等无效区域进行区分，得到的结果虽然对定性分析的伪彩图影响不大，但对定量分析结果则带来较大误差。基于颅内脑组织区域掩膜，对无效区域进行了去除。此外，还依据连通域分析，去除了满足阈值但面积过小(即面积未超时设定要求)的分割区域，从而减小了噪声干扰。最终的容积定量结果则由满足阈值分割的像素点数与步骤1得到的单个像素容积(即像素空间体积)相乘得到。
[0097]
本发明实施例上述方案适用于ct和mri影像，根据二者的成像原理，对执行前述步骤2，获得时域插值处理后的脑灌注影像数据后，需要首先将图像强度转换为造影剂浓度，得到用于步骤3中自动检测算法的时间-浓度曲线。可选地，本发明例中选用ct灌注影像作为示例，转化公式为：
[0098][0099]
其中，μ(j)为原始图像强度，c(j)为灌注造影剂浓度，μ0为基准强度，选取在造影剂到达前，相应位置前b个采样点的均值，角标i,j均指时间索引即时域采样点。
[0100]
相比于其他现有方法，本发明实施例上述方案，通过设计算法实现了全自动选取动脉输入函数及静脉输出函数，并自动判别血管、脑室等无效区域，从而提供了更准确、便捷、快速的灌注参数定性及定量分析，对提升灌注影像阅片质量具有积极影响。
[0101]
本发明实施例上述方案可以适用于多种场景，例如，应用于科研机构或者医疗机构等，通过上述方案自动对灌注影像进行处理，以辅助完成后续的科研工作、或者医护人员上岗培训或考核等。
[0102]
本发明另一实施例还提供一种脑灌注影像的处理系统，其主要用于实现前述实施例提供的方法，如图6所示，该系统主要包括：
[0103]
数据获取与类型识别模块，用于获取原始脑灌注影像数据，并识别数据类型；
[0104]
数据预处理与时域插值处理模块，用于对所述原始脑灌注影像数据进行预处理，并结合数据类型进行时域插值；
[0105]
动脉输入函数自动选择模块，用于通过自动检测算法从时域插值处理后的脑灌注影像数据中自动选择动脉输入函数；
[0106]
灌注参数计算模块，用于基于时域插值处理后的脑灌注影像数据和动脉输入函数，计算灌注参数；
[0107]
影像后处理和定量统计模块，用于从预处理后的原始脑灌注影像数据提取出颅内脑组织区域掩膜，并结合灌注参数对所述预处理后的原始脑灌注影像数据进行后处理和定量统计。
[0108]
所属领域的技术人员可以清楚地了解到，为描述的方便和简洁，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即将系统的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。
[0109]
本发明另一实施例还提供一种处理设备，如图7所示，其主要包括：一个或多个处理器；存储器，用于存储一个或多个程序；其中，当所述一个或多个程序被所述一个或多个处理器执行时，使得所述一个或多个处理器实现前述实施例提供的方法。
[0110]
进一步的，所述处理设备还包括至少一个输入设备与至少一个输出设备；在所述处理设备中，处理器、存储器、输入设备、输出设备之间通过总线连接。
[0111]
本发明实施例中，所述存储器、输入设备与输出设备的具体类型不做限定；例如：
[0112]
输入设备可以为触摸屏、图像采集设备、物理按键或者鼠标等；
[0113]
输出设备可以为显示终端；
[0114]
存储器可以为随机存取存储器(random access memory，ram)，也可为非不稳定的存储器(non-volatile memory)，例如磁盘存储器。
[0115]
本发明另一实施例还提供一种可读存储介质，存储有计算机程序，当计算机程序被处理器执行时实现前述实施例提供的方法。
[0116]
本发明实施例中可读存储介质作为计算机可读存储介质，可以设置于前述处理设备中，例如，作为处理设备中的存储器。此外，所述可读存储介质也可以是u盘、移动硬盘、只读存储器(read-only memory，rom)、磁碟或者光盘等各种可以存储程序代码的介质。
[0117]
以上所述，仅为本发明较佳的具体实施方式，但本发明的保护范围并不局限于此，任何熟悉本技术领域的技术人员在本发明披露的技术范围内，可轻易想到的变化或替换，都应涵盖在本发明的保护范围之内。因此，本发明的保护范围应该以权利要求书的保护范围为准。