低剂量脑灌注CT图像恢复方法：基于先验图像约束扩散张量

牛善洲，刘 宏，刘沛沄，张梦真，邱 洋，黎 钰，谢国强，刘国良，卢绍辉

1赣南师范大学数学与计算机科学学院，江西 赣州 341000；赣南医学院2医学信息工程学院，3第一附属医院，江西 赣州341000

脑卒中（俗称脑中风）是脑部缺血造成的脑组织损伤，因其发病率高、致残率高、死亡率高和复发率高的特点，已成为危害人类健康最为严重的疾病之一［1,2］。脑灌注CT成像技术是脑卒中检查的杰出代表，可以快速、准确、无创地获取常规CT检查所不能获得的三维脑灌注血流动力学参数，已在急性脑中风诊断方面取得了巨大成功［3-7］。然而，脑灌注CT检查中需要对感兴趣层面进行重复动态扫描，使得患者接受X射线的辐射剂量过高，导致患者出现严重脱发、皮肤灼伤、头痛等问题［8-10］。因此，为了使脑灌注CT成像更具有临床生命力，优化控制X射线辐射剂量已经成为灌注CT成像领域亟待解决的关键性课题。

降低管电流或管电压是实现低剂量脑灌注CT成像最简单有效的方式。但是，降低管电流或管电压会导致传统滤波反投影算法重建的图像中含有大量的噪声和伪影，进而影响到脑血流参数图的准确性［11-13］。目前，低剂量脑灌注CT成像方法主要分为两类：低剂量脑灌注CT图像重建方法和低剂量脑灌注CT图像后处理方法。低剂量脑灌注CT重建方法计算量大，耗费时间长，限制了其在临床实践中的应用［14,15］。图像后处理技术因其简单且易于操作已被广泛应用于低剂量脑灌注CT图像恢复。例如，有学者提出的基于低秩和全变分正则化的低剂量脑灌注CT图像的恢复方法［16］，基于各向异性扩散滤波的低剂量脑灌注CT图像的恢复方法［17］，基于双边滤波的低剂量脑灌注CT图像的恢复方法［18］。

脑灌注CT图像之间丰富的结构相关性可以为低剂量脑灌注CT图像恢复提供先验信息。例如，学者提出了PSRR策略，该方法将先前标准剂量扫描得到的CT图像进行配准后用于低剂量图像滤波，该方法对图像配准的精度要求很高［19］。Ma等［20］提出了先验图像引导的非局部均值滤波方法，该方法降低了目标图像和先验图像在配准精度上的要求。Supanich 等［21］提出了HYPR方法，该方法先构造了一个合成图像，并利用该图像与低剂量图像之间的结构相关性对低剂量噪声图像进行滤波。由于先验图像与低剂量脑灌注CT序列图像的衰减值不一致，两幅图像之间的差异会导致滤波后的图像与目标图像产生较大误差。为了解决上述问题，本文提出了基于先验图像约束扩散张量（PICDT）的低剂量脑灌注CT图像恢复方法，该方法可以将高质量先验图像中的互补结构信息引入到低剂量脑灌注CT图像恢复过程中，并且保持目标图像的结构和强度特征。数值体膜和临床数据实验结果表明本文提出的方法在有效去除低剂量脑灌注CT图像噪声和伪影的同时可以较好地保持图像的结构细节特征。

1 方法

1.1 脑灌注CT成像

脑灌注CT成像过程如下：首先进行常规CT平扫，确定脑灌注扫描层面，静脉快速团注对比剂，当对比剂首次通过脑血管及脑实质时，对选定层面进行连续动态CT扫描并重建出对应的时间序列图像，继而获得每个像素的时间-浓度曲线，并根据时间-浓度曲线可以得到脑血流动力学参数图像。低剂量脑灌注CT图像恢复问题可以写成如下形式：

其中，x是待恢复图像，{x1,x2,x3,...,xK}是获取的低剂量脑灌注CT序列图像，xk是第k帧图像，K是序列图像的帧数，β＞0 是正则化参数，R(x)是正则化函数。

1.2 各向异性扩散滤波

各向异性扩散［22］可以写成如下的优化问题：

其中，f是函数空间Ω 上的光滑函数，φ(·)是势函数，的Euler-Lagrange方程为：

其中，div 是散度算子。因此，方程（2）可以通过如下的梯度下降算法求解：

其中，n是迭代步数的指数，τ＞0 是步长。

1.3 基于扩散张量的各项异性扩散滤波

图像u的结构张量［22］可定义为如下的半正定对称矩阵：

其中，uσ=kσ*u，kσ为高斯卷积核，J11,J12,J22是结构张量Jρ中的元素。结构张量（5）的两个正交的特征向量为v1,v2，v1平行于

基于上述结构张量的特征向量，图像u的扩散张量D可以写为：

特征值γ1,2决定了扩散张量在两个特征向量方向v1,2上的平滑强度。γ1,2可以由如下公式给出：

其中，ε＞0 控制法向量扩散的强度。各向异性扩散过程(3)可以用扩散张量D写成如下形式：

1.4 基于先验图像约束扩散张量的低剂量脑灌注CT图像恢复

扩散张量将先验图像μ中的几何结构特征嵌入到待恢复的低剂量图像xk中，利用高质量先验图像μ中的互补信息来增强目标图像xk的质量，同时保持目标图像xk的重要特征。我们构建先验图像约束的扩散张量：

其中，Dxk是图像xk的扩散张量，Dμ是图像μ的扩散张量，s∈[0,1] 决定了两个扩散张量的相对贡献。在本研究中，s由如下公式给出［23］：

其中，abs(·)为绝对值函数，v1(μ)·v1(xk)是图像μ和图像xk归一化特征向量的内积。

综上所述，基于PICDT的低剂量脑灌注CT图像恢复方法的计算步骤如下：

初始化：迭代初始值设置为待恢复的低剂量CPCT图像xk，迭代步长设置为τ；

步骤1：分别计算第n次迭代后目标图像和先验图像μ的扩散张量根据式（10）计算得到新的扩散张量

1.5 对比方法

为了验证PICDT方法的有效性，我们将PICDT方法与滤波反投影（FBP）算法、未使用先验图像的扩散张量（DT）方法和先验图像约束的全变分（PICTV）方法进行比较。当公式（10）中的s=1时，PICDT方法即退化为DT方法。PICTV方法的目标函数为：

其中，x是待恢复图像，{x1,x2,x3,...,xK}是获取的低剂量脑灌注CT序列图像，xk是第k帧图像，K是序列图像的帧数，β＞0 是正则化参数，μ是先验图像，TV(x-μ)是x-μ的全变分。

1.6 参数选择

PICDT方法的参数有步长τ，高斯卷积核kσ，控制法向量强度ε。对于步长τ的选取，我们采用试误法，实验结果表明步长在1.0×10-2≤τ≤1.0×10-1的区间内恢复后的图像效果较好。我们跟据Kazantsev等［23］的工作确定了高斯卷积核kσ的大小。ε决定了扩散强度，并将ε设置为FBP图像梯度直方图的90%［24］。DT方法的参数选取方式与PICDT方法一致。PICTV方法中的参数β决定了恢复图像的光滑程度，其值与待恢复图像的噪声水平有关。参数β的选取也采用试误法，我们将恢复后的图像与参考图像进行比较后选取最优的β值。在数值体膜实验中，PICDT方法的τ=0.08，DT方法的τ=0.03，PICTV方法的β=9.3×10-5。在临床数据实验中，PICDT 方法的τ=0.055，DT 方法的τ=0.04，PICTV方法的β=15。

2 结果

2.1 数值体膜实验

我们使用Shepp-Logan体膜仿真得到低剂量脑灌注CT 图像，相应的成像几何参数为：（1）投影角度为1160，每个投影角度的探测元个数为672，每个探测器单元的大小为1.4077 mm；（2）放射源到探测器的距离为1040 mm；（3）射线源到旋转中心的距离为570 mm。图像维数为512×512 。X 射线的入射光子数目为1×105。

图1给出了Shepp-Logan体膜实验结果。其中，第1行是体膜图像，第2行是FBP算法的图像，第3行是DT方法恢复的图像，第4行是PICTV方法恢复的图像，第5行是PICDT方法恢复的图像；从左到右第1列为第#5帧图像，第2列为第#15帧图像，第3列为第#25帧图像。图2为图1中第#15帧图像的感兴趣区域（图1中红色方框所示）的放大图。从图2可以看出，PICDT方法恢复的图像与其它方法恢复的图像相比拥有更多的图像结构特征。进一步，图3给出了图1中结果的UQI［25］值，FBP，DT，PICTV，PICDT 方法的UQI 平均值分别为0.6900、0.7909、0.7980、0.8083。PICDT方法恢复的图像的UQI值最高，即最接近体膜图像。

图3 Shepp-Logan体膜结果的UQI值Fig.3 UQI values of the Shepp-Logan phantom results.

脑血流动力学参数能够精准反映脑组织的灌注情况，为临床诊断和治疗提供重要的参考价值。我们根据文献［26］中的灌注模型，使用截断奇异值分解算法［27］求解脑血流量（CBF）参数图像。图4为Shepp-Logan体膜数据的CBF参数图像，其中，图4A为体膜图像的CBF参数图像，图4B为FBP算法的CBF参数图像，图4C为DT方法的CBF参数图像，图4D为PICTV方法的CBF参数图像，图4E为PICDT方法的CBF参数图像。从图4中可以看出，PICDT 方法的CBF 参数图像与体膜的CBF参数图像最接近。为进一步客观评价本文方法，我们用FSIM（feature similarity indexing method）［28］值和SSIM（structured similarity indexing method）［28］值 对CBF参数图像进行定量分析。图5给出了图4中CBF参数图像的FSIM值和SSIM值，FBP算法的CBF参数图像的FSIM值和SSIM值分别为0.7146，0.5868；DT方法的CBF 参数图像的FSIM 值和SSIM 值分别为0.7296，0.8091；PICTV方法的CBF参数图像的FSIM值和SSIM 值分别为0.7355，0.8214；PICDT 方法的CBF参数图像的FSIM 值和SSIM 值分别为0.7383，0.9247。PICDT方法的CBF图像的FSIM和SSIM值最高，即最接近体膜的CBF参数图像。

2.2 临床数据实验

本实验数据采集自临床脑灌注CT扫描，图像大小为512×512。选取正常剂量的脑灌注CT图像作为金标准，并用2.1节中的方法仿真得到低剂量脑灌注CT图像。

图6为临床数据实验结果，第1行是常规剂量图像，第2行是低剂量脑灌注CT图像，第3行是DT方法恢复的图像，第4 行是PICTV 方法恢复的图像，第5 行是PICDT方法恢复的图像；从左到右第1列为第#5帧图像，第2列为第#15帧图像，第3列为第#25帧图像和第4列为第#35帧图像。图7为图6中第#25帧图像的感兴趣区域（图6中红色方框所示）的放大图。从图7可以看出，PICDT方法恢复的图像与DT方法恢复的图像相比拥有更多的图像结构特征。图8给出了图6中结果的UQI值，FBP，DT，PICTV，PICDT方法恢复图像的UQI平均值分别为0.8414、0.8605、0.8750、0.8836。PICDT方法得到的图像的UQI值最高，即最接近常规剂量的脑灌注CT图像。

图8 临床数据实验结果的UQI值Fig.8 UQI values of the clinical data results.

图9为临床数据的CBF参数图像。其中，图9A为常规剂量图像的CBF参数图像，图9B为FBP算法的CBF参数图像，图9C为DT方法的CBF参数图像，图9D为PICTV方法的CBF参数图像，图9E为PICDT方法的CBF参数图像。从图9中可以看出，我们提出的PICDT方法的CBF参数图像与常规剂量图像的CBF参数图像最接近。为进一步评价本文方法的有效性，我们使用FSIM值和SSIM值对脑血流动力学参数图进行定量分析。图10为由图9中结果的FSIM值和SSIM值。图10中FBP所得CBF参数图像的FSIM值和SSIM值分别为0.7769，0.5299；DT所得CBF参数图像的FSIM值和SSIM 值分别为0.8456，0.6863；PICTV 所得CBF参数图像的FSIM值和SSIM值分别为0.8566，0.7129；PICDT所得CBF参数图像的FSIM值和SSIM值分别为0.8689，0.7731。从图10中可以看出PICDT方法恢复的CBF参数图像的FSIM值和SSIM值均高于FBP、DT和PICTV方法，即最接近常规剂量得到的CBF参数图像。

图9 临床实验图像的CBF参数图Fig.9 CBF maps obtained from clinical dataset.A:CBF maps of normal-dose CT images;B:CBF maps of FBP images;C:CBF maps of DT images;D:CBF maps of PICTV images;E:CBF maps of PICDT images.

图10 临床数据的CBF参数图的FSIM值和SSIM值Fig.10 FSIM and SSIM values of the images in Fig.9.

3 讨论

目前，使用先验图像引导的低剂量脑灌注CT图像恢复方法［19-21］是利用序列图像之间的结构相关性，将高质量图像与低剂量脑灌注CT图像配准后进行滤波，却忽略了高质量图像与目标图像之间强度的差异，从而导致恢复图像中出现错误的组织结构［20］。本文中，为了解决先验图像与目标图像强度不一致问题，我们提出了基于先验图像约束扩散张量的低剂量脑灌注CT图像恢复方法。本文方法首先分别计算先验图像和目标图像的扩散张量，并构造一个新的联合扩散张量，然后利用该扩散张量将先验图像中的几何结构特征引入到低剂量脑灌注CT图像恢复过程中，在有效去除低剂量脑灌注CT图像噪声和伪影的同时可以保持图像的结构细节特征。

在Shepp-Logan体膜实验中，从图1中可以看出，PICDT方法在噪声抑制效果上佳，并且图2进一步表明PICDT 方法在保持结构细节特征方面也有良好的表现。图3定量结果表明，DT方法的平均UQI值比FBP算法的平均UQI 值提高了14%，PICTV 方法的平均UQI 值比FBP 算法的平均UQI 值提高了16%，而PICDT方法的平均UQI值比FBP方法的平均UQI值提高了17%。从图4可以看出，FBP方法的CBF参数图像存在大量的噪声，DT方法和PICTV方法的CBF参数图像虽然有所改善，但其与体膜图像的仍然存在一定差距，而PICDT方法的CBF参数图像与体膜图像的CBF参数图像极为接近。我们还分别计算了各方法的FSIM和SSIM进一步表明了PICDT方法的CBF参数图像恢复效果最好（图5）。在临床数据实验中，从图6可以看出，PICDT方法在噪声抑制方面表现更加优秀。从图7可以看出，虽然PICTV方法在抑制噪声方面的表现比DT方法更优，但PICTV方法因先验图像和低剂量图像之间的强度差异产生了部分伪影并抹去了部分特征，从而导致PICTV方法的感兴趣区域图像恢复效果在视觉上比DT 方法更差，而PICDT 方法未产生伪影且与FBP，DT和PICTV方法相比图像恢复效果更加优秀。图8为各方法的UQI，其中，DT方法的平均UQI值比FBP算法的平均UQI值提高了2%，PICTV方法的平均UQI值比FBP算法的平均UQI值提高了4%，而PICDT方法的平均UQI值比FBP方法的平均UQI值提高了5%。从图9可以看出，PICDT方法的CBF参数图像最接近常规剂量图像的CBF参数图像。图10定量结果表明，PICDT方法的CBF参数图像具有最高的FSIM和SSIM，这进一步表明，PICDT方法恢复的低剂量图像的CBF参数图像与高剂量图像的CBF参数图像在结构和细节方面更为接近。

本文方法仍然有一定的缺陷，对不同的目标图像需要人工选择正则化参数，限制了其更为广泛的应用。如何自适应选择正则化参数一直以来都是CT成像中公开的难问题，本文下一步研究工作便是如何自适应选择正则化参数。