基于能谱CT的材料组分彩色表征研究
孔慧华1,2, 连祥媛1, 陈平2, 潘晋孝1,2
1.中北大学理学院, 山西 太原 030051
2.信息探测与处理山西省重点实验室, 山西 太原 030051

作者简介: 孔慧华, 女, 1977年生, 中北大学理学院副教授 e-mail: huihuak@163.com

摘要

基于光子计数探测器的X射线能谱CT, 通过增加能谱分辨率实现了CT图像由灰度向彩色的转变, 提高了材料识别能力。 然而随着能谱通道数量的增加, 通道中的噪声也显著增加, 降低了材料识别的准确性。 为充分利用能谱CT图像的稀疏性以及能谱通道之间图像的相关性, 提出一种多约束窄谱CT迭代重建算法, 可在降低噪声的同时有效保留图像的边缘及细节部分。 进一步利用主成分分析对窄谱CT图像中的能谱信息进行估计, 并建立主成分图像与R, G和B颜色分量之间的映射关系, 最后获取彩色CT图像。 该方法通过材料组分的彩色表征可以有效识别材料, 同时降低图像中的背景噪声。 仿真实验和实际实验结果验证了多约束窄谱CT迭代重建算法的有效性以及利用主成分分析进行材料组分彩色表征的可行性。

关键词: 能谱CT; 稀疏表示; 相关性; 主成分分析; 彩色表征
中图分类号:TP391.4 文献标志码:A
Research on Color Characterization of Material Components Based on Spectral CT
KONG Hui-hua1,2, LIAN Xiang-yuan1, CHEN Ping2, PAN Jin-xiao1,2
1. School of Science, North University of China, Taiyuan 030051, China
2. Shanxi Key Laboratory of Signal Capturing & Processing, Taiyuan 030051, China
Abstract

Photon-counting detector based X-ray spectral computed tomography (CT), realizes the transformation of CT image from gray to color by increasing energy resolution, which increases material identification capability. However, with increasing the number of energy channels, the channel’s noise increases significantly, which decreases the accuracy of material identification. In order to make full use of the sparsity of spectral CT images and the correlation between spectral CT images, a multi constraint narrow-spectral CT iterative reconstruction algorithm is proposed, which can effectively preserve the edges and details of the image while reducing the noise. Furthermore, principal component analysis (PCA) is used to estimate the spectrum information in narrow spectrum CT images, and the mapping relationship between principal component image and color components R, G, B are established. Finally, the color CT image is obtained. This method can effectively identify materials through the color representation of material components and reduce the background noise in the images. The results of simulation and practical experiments show the proposed reconstruction algorithm is effective, and it is feasible to use PCA for the color characterization of material components.

Keyword: Spectral CT; Sparse representation; Correlation; Principal component analysis; Color characterization
引言

基于光子计数探测器(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图像重建可表示为

minfEAfE-PE2(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,  1mM,  1nN,  J=M×N(2)

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

minfE,  D,  {αj}μ2AfE-PE22+λΦ1(fE)+β2Φ2(fE)+Φ3(fE)(3)

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

Φ1(fE)fETV(4)

Φ2(fE)jEjfE-Dαj22+jγjαj0(5)

Φ3(fE)-jcov(EjfE,  Ejf先验)σ(EjfE)σ(Ejf先验)(6)

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

fETV=m=1Mn=1N(xfE)m, n2+(yfE)m, n2(7)

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

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

ρ=cov(EjfE, Ejf先验)σ(EjfE)σ(Ejf先验)=m=1ndn=1nd(umn-u̅)(vmn-v̅)m=1ndn=1nd(umn-u̅)2m=1ndn=1nd(vmn-v̅)2(8)

式(8)中, u̅v̅分别表示图像块EjfEEjf先验中像素的平均灰度值。

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

minfE, D, {αj}μ2AfE-PE22+λΦ1(fE)+β2Φ2(fE)(9)

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

minfEΦ3(fE)(10)

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

第二个子问题是对第一个子问题求解的结果进行负相关系数最小化处理, 采用梯度下降法求解, 设第一个子问题求解的结果为 fE(1), 则对子块Ej f(1)E

umn(l+1)=umn(l)-η-cov(u(l), v)σ(u(l))σ(v)umn(11)

式(11)中, η 为梯度下降算法中的步长, u(0)=Ej fE(1)

由于(子)块之间存在像素重叠的现象, 记录子块提取过程中像素重叠的次数, 将其作为权重, 最后通过加权平均方法取得像素点的平均值。 将子块Ej fE(1)经过式(11)运算的结果记为Ej f^E(1), 则第二个子问题求解的结果见式(12)

fE(2)=jEjf^E(1)jETjEj(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≤ iL)个主成分图像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 先验图像
(a): 小鼠胸腔切片; (b): 临床前小鼠
Fig.1 The prior images
(a): The phantom of mouse thoracic caoity; (b): The preclinical mouse

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

图2 四种算法下小鼠胸腔的重建图像
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC 从上往下依次为第一、 第二和第三能谱通道
Fig.2 The reconstruction images of mouse thoracic cavity by four algorithms
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC From top to bottom, the rows are 1st, 2nd and 3rd energy channel

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

表1 四种算法在三个能谱通道下小鼠胸腔重建图像的数量性评价指标 Table 1 Quantity evaluation of four algorithms for mouse thoracic cavity reconstruction in the three energy channels

为获取材料组分的彩色表征, 首先对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 小鼠胸腔的主成分分析图像
(a): PCA-1; (b): PCA-2; (c): PCA-3; (d): PCA-2的平方; (e): PCA-3的平方
Fig.3 PCA images of mouse thoracic cavity
(a): PCA-1; (b): PCA-2; (c): PCA-3; (d): The square of PCA-2; (e): The square of PCA-3

图4 四种重建算法下小鼠胸腔材料组分的彩色表征
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC
Fig.4 Color characterization of material components for mouse thoracic cavity by four reconstruction algorithms
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC

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

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 四种算法下临床前小鼠的重建图像
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC 自上而下依次为第1、 第6和第13能谱通道
Fig.5 The reconstruction images of preclinical mouse by four algorithms
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC From top to bottom, the rows are 1st, 6th and 13th energy channel

对能谱通道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 临床前小鼠的主成分分析图像
(a): PCA-1; (b): PCA-2; (c): PCA-3; (d): PCA-2的平方; (e): PCA-3的四次方
Fig.6 PCA images of preclinical mouse
(a): PCA-1; (b): PCA-2; (c): PCA-3; (d): The square of PCA-2; (e): The fourth pwoer of PCA-3

图7 四种重建算法下临床前小鼠材料组分的彩色表征
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC
Fig.7 Color characterization of material components for preclinical mouse by four reconstruction algorithms
(a): TV; (b): DL; (c): IPCC; (d): TV-DL-IPCC

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

4 结论

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

参考文献
[1] Wu Weiwen, Yu Haijun, Chen Peijun, et al. Physics in Medicine and Biology, 2020, 65(24): 245006. [本文引用:1]
[2] Nicolas Ducros, Juan Felipe Perez-Juste Abascal, Bruno Sixou, et al. Medical Physics, 2017, 44(9), e174. [本文引用:1]
[3] Zhang Zheng, Chen Buxin, Xia Dan, et al. Medical Image Analysis, 2021, 70: 102030. [本文引用:1]
[4] Kanii Yoshinori, Ichikawa Yasutaka, Nakayama Ryohei, et al. Japanese Journal of Radiology, 2020, 38(3): 215. [本文引用:1]
[5] Wang Shaoyu, Wu weiwen, Feng Jian, et al. Physics in Medicine and . Biology. 2020, 65(24): 245005. [本文引用:1]
[6] Dong Yiqiu, Hansen Per Christian, Kjer Hans Martin. IEEE Transactions on Computational Imaging, 2018, 4(4): 528. [本文引用:1]
[7] Niu Shanzhou, Yu Gaohang, Ma Jianhua et al. Inverse Problems, 2018, 34(2): 024003. [本文引用:1]
[8] Rigie David S, La Riviere Patrick J. Physics in Medicine and Biology, 2015, 60(5): 1741. [本文引用:1]
[9] Xie Huiqiao, Ren Yan, Long Wenting, et al. IEEE Transactions on Biomedical Engineering, 2021, 68(3): 1074. [本文引用:1]
[10] Kong Huihua, Lei Xiaoxue, Lei Lei, et al. IEEE Access, 2020, 8: 133367. [本文引用:1]
[11] Kong Huihua, Liu Rui, Yu Hengyong. Journal of X-ray Science and Technology, 2016, 24(2): 221. [本文引用:1]
[12] Wang Miaoshi, Zhang Yanbo, Liu Rui, et al. Physics in Medicine & Biology, 2016, 61(24): 8699. [本文引用:1]