基于能谱CT的材料组分彩色表征研究

孔慧华, 连祥媛, 陈 平, 潘晋孝

1.中北大学理学院，山西 太原 030051 2.信息探测与处理山西省重点实验室，山西 太原 030051

引 言

基于光子计数探测器(photon-counting detector, PCD)的X射线能谱CT，通过设置能谱阈值，可以同时获取多个窄谱通道的投影数据，在材料组分识别方面具有优异表现[1-2]。然而，由于窄谱通道中光子数骤减，噪声水平增加，显著降低了分解材料图像的信噪比，影响了能谱CT的实际应用。因此，降低能谱CT的重建噪声对提高材料分解精度具有重要的现实意义。

在能谱CT中，每个窄谱通道中的图像都可以通过稀疏变换，如离散余弦变换和字典等被稀疏表示。基于全变差(total variation, TV)最小化和字典学习(dictionary learning, DL)的图像重建算法已被成功应用于CT[3-4]和能谱CT中[5-6]。对于能谱CT，不同能谱通道的重建图像具有很强的相关性。许多算法同时考虑了图像的稀疏性和不同能谱通道之间图像的相关性，将低秩先验信息、多通道梯度向量或张量与图像稀疏性结合起来，有效提高了重建图像的质量[7-8]。

在材料分解和识别方面，目前主要有前处理和后处理两种方法。这些方法虽然可以对材料组分进行定量分析，但需要知道探测器响应或者被分析材料的先验信息。主成分分析(principal components analysis, PCA)是一种很好的多元数据分析技术，可以应用于任意维度的图像数据，因此可以用来处理多能谱CT数据[9]。

基于以上分析，首先在能谱CT图像稀疏性的基础上，融入窄谱图像间的相关性，提出一种多约束窄谱CT迭代重建算法。然后利用PCA估计能谱信息，通过对各个主成分图像进行分析，将其函数映射为彩色图像的R，G，B分量，最后获取材料组分的彩色表征。

1 多约束窄谱CT迭代重建算法

能谱CT将X射线宽能谱分布划分成多个不相重叠的窄谱通道，一次扫描可以得到多个窄谱通道的投影数据{PE}，每个通道对应一个窄谱图像fE。窄谱CT图像重建可表示为

(1)

式(1)中，A=(aij)I×J是投影矩阵，I表示投影射线总数，J表示图像像素总数；fE是窄谱重建图像；PE是窄谱投影数据。为了处理方便，重建图像fE=(fj)J×1也可以表示为fE=(fm，n)M×N，其中M，N为图像的行与列，其相互关系可表示为式(2)

fj=fm，n，j=(m-1)×N+n，

1≤m≤M，1≤n≤N，J=M×N

(2)

为了充分利用窄谱CT图像的稀疏性以及窄谱通道间图像的相关性，将图像的TV，DL和图像块之间的相关系数(image patch correlation coefficient，IPCC)整合到重建算法中，提出一种基于先验信息的多约束窄谱CT迭代重建方法，记为TV-DL-IPCC，目标函数的定义如式(3)

(3)

式(3)中，μ,λ,β为平衡保真项与正则化项的参数，Φ1(fE)，Φ2(fE)，Φ3(fE)分别定义为式(4)—式(6)

Φ1(fE)≜‖fE‖TV

(4)

(5)

(6)

Φ1(fE)表示图像fE的TV范数，表达式为

(7)

Φ2(fE)表示图像fE的基于字典学习的正则化项，图像块EjfE是由算子Ej∶RJ→RNd从fE中提取的尺寸为Nd=nd×nd的小图像块，其左上角位置为第j个像素。字典D∈RNd×K是一个矩阵，每一列dk∈RNd称为一个原子，k=1，2，…，K。一般来说，字典是冗余的，即Nd≪K。αj是在字典D下图像块EjfE的稀疏表示，γj为基于字典学习的正则化参数。

对企业i来说，在积分交易价格pφ一定时，βi越大意味着研发效率越低，企业续航技术研发的积极性越低。然而，当积分交易价格提高时，企业续航技术研发的积极性提高，通过续航能力提升获取积分进而出售获利。当企业间进行研发合作时，目标函数是总体利润最大化，因而会考虑本企业对另一企业的技术溢出，溢出率越高，花费相同的成本可以更大程度地提高汽车产品的续航能力；而在进行研发竞争时，企业只考虑自身的利润最大化，不考虑本企业对另一企业的技术溢出。

Φ3(fE)表示能谱图像之间的相关性，式(6)中相关系数前面取负号，是为了与TV和DL正则化项中的非零系数最小化一致。虽然不同能谱通道下的重建图像{fE}的衰减系数并不相同，但其结构信息却是高度相关的。因此，在能谱CT重建中可以用已知的高质量重建图像f先验作为参考，发挥先验信息的作用来约束所有能谱通道的重建图像。选择fE和f先验的相关系数作为衡量通道之间图像相关性的正则化项。为了保证微小结构的相似性，采用块策略，类似于字典学习，首先在窄谱图像和先验图像中提取块EjfE和Ejf先验，分别记为u=(umn)nd×nd和v=(vmn)nd×nd，则图像块之间的相关系数ρ可以表达为式(8)

(8)

目标函数(3)可以利用交替最小化迭代方法求解，将其分为两个子问题，第一个子问题为[见式(9)]

(9)

第二个子问题为[见式(10)]

(10)

第一个子问题的求解通过图像重建和字典学习两个步骤完成。在图像重建阶段，采用Split-Bregman算法求解不可微正则化问题，首先引入辅助变量，将目标函数分为L1范数和L2范数两部分，然后采用Bregman算法求解。在字典学习阶段，为了尽量减小先验信息对重建图像质量的影响，采用局部自适应字典，对重建后的图像进行字典学习。式(9)的具体求解流程见参考文献[10]。

(11)

(12)

注意，若(3)式中只含正则化项Φ1(fE)，则为TV算法；若只含Φ2(fE)，则为DL算法，若只含Φ3(fE)，则为IPCC算法。

2 基于主成分分析的材料组分彩色表征

PCA将一组高度相关的变量转换为少数几个互不相关的综合变量，可以最大限度地表示数据中的方差。本研究感兴趣的是将能谱CT数据中不同类别材料之间的对比方差最大化，因此选取主成分对窄谱CT图像进行分析。

将每一个窄谱通道的重建图像fEl看做一个含有J个分量的向量Cl，l=1，2，…，L，所有窄谱通道下的图像{fE1，fE2，…，fEL}构成矩阵C，其中clj表示第l个能谱通道下第j个像素的衰减系数。主成分分析首先计算这组数据的协方差矩阵S=(sij)L×L；然后计算协方差矩阵的特征值和对应的特征向量，设其特征值为λ1≥λ2≥…≥λL，对应的特征向量为ξi=(ξi1，ξi2，…，ξiL)，i=1，2，…，L。则第i(1≤i≤L)个主成分图像PCA-i见式(13)主成分图像

zi=ξi1fE1+ξi2fE2+…+ξiLfEL

(13)

第一主成分在不同材料对比度上有最大方差，其余主成分的方差依次减少。最后建立主成分图像与R，G和B色彩分量之间的映射关系，获取各组分的彩色表征。

3 实验部分

本研究的主要目的是对材料组分进行彩色表征，将通过仿真实验和实际实验验证提出算法的有效性。为了验证所提出的多约束窄谱CT迭代重建算法的性能，选择TV，DL和IPCC算法作为比较算法。所有的算法都是用Matlab和C++的混合编程模式实现的，接口在Matlab中，大规模计算部分在C++中实现，并通过MEX函数进行编译。

3.1 仿真实验

仿真实验选取由MOBY软件生成的数字小鼠胸腔切片作为测试模型，模型尺寸为20 mm×20 mm，分辨率为512×512。实验中，采用等距扇形束扫描，扫描半径为100 mm，位于物体中心的虚拟探测器长度为20 mm，共有320个探测器单元。设置电压为50 kVp，将该宽能谱分为三个能谱通道: {17 keV—28 keV}，{29 keV—35 keV}，{36 keV—50 keV}。每个通道内光子数为10万，生成具有泊松分布的噪声投影数据，其期望是相应的无噪声情况下接收到的光子数。在[0，21π]内均匀采样，采集的投影视角数为360。利用归一化均方根误差(normalized root mean square error, NRMSE)、结构相似度(structural similarity, SSIM)和峰值信噪比(peak signal to noise ratio, PSNR)对算法的性能进行定量评价。

TV-DL-IPCC算法中的参数是在窄谱通道中通过实验使其性能最优获取的。为了便于比较，TV，DL和IPCC算法中使用与TV-DL-IPCC算法相同的参数，同时每个通道都使用相同的参数，其中μ=100,λ=20,β=10,η=0.02。在IPCC和TV-DL-IPCC算法中，将噪声投影下采用Split-Bregman(SB)算法[11]重建的宽谱图像作为先验图像，如图1(a)所示，该图像与宽谱下无噪声重建图像相比较，其数量性评价指标为NRMSE=0.018 8，SSIM=0.999 8，PSNR=50.470 4。

图1 先验图像

TV，DL，IPCC和TV-DL-IPCC算法在三个能谱通道下经过20次迭代的重建结果如图2所示。从图2可以看出，IPCC算法边缘重建效果很好但没有去噪的功能，TV和DL都可以去除噪声。由于TV假设图像是分片光滑的，导致重建图像存在块状伪影，而因为DL是逐图像块进行处理的，所以平滑效果较好但容易将微弱细节平滑掉。TV-DL-IPCC算法不仅可以有效去除噪声还能加强边缘和细节的重建。

图2 四种算法下小鼠胸腔的重建图像

表1给出了TV，DL，IPCC和TV-DL-IPCC算法在三个能谱通道下的NRMSE，SSIM和PSNR值。从表1可以看出，与其他方法相比TV-DL-IPCC具有较小的NRMSE，较大的SSIM和PSNR。

表1 四种算法在三个能谱通道下小鼠胸腔重建图像的数量性评价指标

为获取材料组分的彩色表征，首先对TV-DL-IPCC算法重建的三个能谱通道下的CT图像进行主成分分析，图3(a)，(b)和(c)给出了三个主成分的图像。第一主成分图像包含了能谱CT图像的99.00%的信息量，表示三个能谱通道CT图像的平均，如图3(a)所示；第二主成分图像在骨骼上取值为负，碘对比剂上取值为正，且包含了少量的背景噪声，如图3(b)所示；第三主成分图像在碘对比剂上取值为负，且包含了大量的背景噪声，如图3(c)所示。通过分析发现背景噪声在零附近波动，为了去掉背景噪声，且突出骨骼和碘对比剂，对第二和第三主成分图像中的像素值进行平方，得到主成分函数的图像，如图3(d)和(e)所示。然后分别将图3(a)，(d)和(e)映射为彩色图像的G，R和B颜色分量，得到小鼠胸腔各组分的彩色表征，如图4所示。

图3 小鼠胸腔的主成分分析图像

从图4可以看出，由IPCC算法获取的彩色CT图像边缘清晰但噪声明显，TV，DL和TV-DL-IPCC算法下的能谱CT图像都很好地实现了材料组分的彩色表征，背景是黑色，软组织是绿色，骨骼根据其硬度在绿色到黄色之间，碘对比剂是紫色。其中TV算法下的彩色表征有一些块状伪影，DL算法下的彩色表征有效去除了块状伪影，但胸腔内的个别软组织也被平滑掉了，而TV-DL-IPCC算法下的彩色表征不仅有效抑制了噪声，且边缘、细节清晰。从图4(a—d)中还可以看到提出的彩色映射函数很好地去除了图像中的背景噪声。

图4 四种重建算法下小鼠胸腔材料组分的彩色表征

3.2 临床前小鼠实验

为了进一步验证提出算法在实际应用中的有效性，将一组真实的小鼠临床前能谱CT投影用于算法的测试。该组数据由美国俞恒永博士所在团队提供，采用新西兰Medipix3探测器[12]。实验用的电压为120 kVp，电流为175 mA。从探测源到系统中心的距离为158 mm，到探测器的距离为255 mm。整个能谱范围内收集了13个能谱通道的投影数据，处理后的投影图像分辨率为360×512，重建图像的分辨率为512×512。在IPCC和TV-DL-IPCC算法中，使用SB算法重建的宽谱图像作为先验图像，如图1(b)所示。与模拟仿真实验类似，算法中使用的参数是通过实验优化获取的，且用于所有能谱通道，其中μ=100，λ=50，β=5，η=0.02。

图5(a—d)分别给出了TV，DL，IPCC和TV-DL-IPCC算法在能谱通道1，6，13下经过20次迭代的重建结果。由图5(a—d)分别可以看到与仿真实验相似的结果，TV重建图像存在块状伪影，DL算法具有较好的去噪效果，IPCC重建图像边缘清晰，但对图像的噪声无能为力。而TV-DL-IPCC算法结合了TV，DL和IPCC的优点，在图像边缘保持以及去噪方面明显优于其他算法。

图5 四种算法下临床前小鼠的重建图像

对能谱通道1，6和13重建的CT图像进行主成分分析，得到的三个主成分图像如图6(a)，(b)和(c)所示。第一主成分图像包含了能谱CT图像的99.49%的信息量，表示三个能谱通道CT图像的平均；第二主成分图像在骨骼上取值为正，且包含了少量的背景噪声；第三主成分图像包含了大量的背景噪声。为了去掉背景噪声，且突出骨骼，对第二和第三主成分图像中的像素值分别进行平方和四次方，得到主成分函数的图像，如图6(d)和(e)所示。然后分别将图6的(a)，(d)和(e)映射为彩色图像的G，R和B颜色分量，得到临床前小鼠各组分的彩色表征，分别如图7(a—d)所示。

图6 临床前小鼠的主成分分析图像

图7 四种重建算法下临床前小鼠材料组分的彩色表征

由图7(a—d)可以看出通过主成分分析，四种算法下的能谱CT图像都很好地实现了材料组分的彩色表征，背景是黑色，软组织是绿色，骨骼根据其硬度在绿色到黄色之间，且很好地去掉了重建图像中的背景噪声。从图7(d)可以看到TV-DL-IPCC算法下的彩色表征不仅边缘清晰，且有效去除了噪声的影响。同时从图7(a)中可以看出TV算法中的块状伪影在彩色映射过程中也被去掉了。

4 结 论

在能谱CT中，每个能谱通道的重建图像都可以被稀疏表示，不同能谱通道下的重建图像具有很强的相关性。利用这些先验信息可以有效地提高重建图像的质量，本研究将TV，DL与IPCC约束项相结合，提出了一种多约束窄谱CT迭代重建算法。实验结果表明，该算法在去噪的同时，很好地保留了图像的边缘和细节特征。为了实现CT图像组分的彩色表征，采用主成分分析的方法对能谱CT图像进行处理，通过建立主成分图像与彩色图像R，G，B分量之间的映射关系，最终获取各材料组分的彩色表征。实验结果表明该方法在获取彩色表征的同时，还可以有效去除噪声的影响。