基于MAP-EM算法的双能CT直接迭代基材料分解方法 [PDF全文]
(1南京航空航天大学机械结构力学及控制国家重点实验室, 南京210016)

为了提高双能CT基材料分解的精度,降低基材料图像的噪声,提出了基于MAP-EM算法的直接迭代基材料分解方法.结合MAP-EM算法,推导出基材料分解直接迭代求解公式,基于双能投影数据集直接重建基材料分解图像,并对该方法的性能进行了评价和分析.仿真结果表明,所提方法可显著降低分解误差和基材料图像噪声,提高对比噪声比.与基于FBP算法的图像域基材料分解方法相比,该方法可使基材料图像中各材料区域的噪声水平下降57.42%~63.64%,分解误差水平降低31.72%~62.14%,对比噪声比提高1.37%~223.17%.

Direct iterative basis material decomposition method for dual-energy CT based on MAP-EM algorithm
Zhou Zhengdong1,Zhang Xuling1,Xin Runchao1,2,Jia Junshan1,Wei Shisong 1,Mao Ling1
(1State Key Laboratory of Mechanics and Control of Mechanical Structures, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)(2Department of Nuclear Science and Engineering, Nanjing University of Aeronautics and Astronautics, Nanjing 210016, China)

To improve the accuracy of basis material decomposition and reduce the noise in the basis material images for dual-energy CT(computed tomography), a direct iterative basis material decomposition method based on the MAP-EM(maximum a posteriori expectation-maximization)algorithm was proposed. Combined with the MAP-EM algorithm, the direct iterative formulas for basis material decomposition were derived. The basis material decomposition images were directly reconstructed from the dual-energy projection data set. The performance of the method was evaluated and analyzed. The simulation results show that the proposed method can significantly reduce the error of decomposition and the noise of basis material images, and improve the contrast-to-noise ratio. Compared with the image domain basis material decomposition method based on the FBP(filtered back projection)algorithm, the proposed method can reduce the noise of each material region in the basis material images by 57.42% to 63.64%. The decomposition error levels of each material region in the basis material images can be reduced by 31.72% to 62.14%. The contrast-to-noise ratios of each material region in the basis material images can be improved by 1.37% to 223.17%.

引言

1976年Alvarez和Macovski[1]首次提出双能CT(dual-energy CT, DECT)的概念,利用2种不同管电压下的X射线对目标进行扫描成像.利用双能CT可获得伪单能图像,降低射束硬化伪影的影响[2].双能CT通过识别2个能谱下的能量信息,可以提供更多的物质属性,在物质识别方面远远优于传统CT.

目前,基材料分解方法可分为2类:①两步法,包括基于图像域的基材料分解方法[3]和基于投影域的基材料分解方法[4]; ②直接迭代基材料分解方法[5],通过建立投影数据与基材料分解系数图像之间的直接关系,根据投影数据直接求解基材料图像.对于两步法,由于第1步难以实现一一映射(bijection),从而导致信息丢失,而第2步无法对第1步的缺失信息进行弥补[6],因此分解结果不可避免地存在误差.直接迭代基材料分解方法可克服两步法的固有缺陷,提高分解精度[7].此外,相关直接迭代投影函数模型以及其求解方法也被相继提出.前期研究表明,最大后验概率期望最大(maximum a posteriori expectation-maximization, MAP-EM)算法可有效提高能谱CT重建图像的质量[8].

本文针对双能CT提出了基于MAP-EM算法的直接迭代基材料分解方法(direct iterative basis material decomposition based on MAP-EM algorithm, MAP-EM-DIBMD),推导出迭代求解公式,基于双能投影数据集直接重建基材料图像,并与基于FBP算法的图像域基材料分解方法(basis material decomposition in image domain based on FBP algorithm, FBP-IDBMD)[9]进行性能比较与分析.

1 材料与方法1.1 MAP-EM统计重建算法

在实际应用中,由于受到非相干散射、脉冲堆积、电子噪声等因素的影响,投影数据中不可避免会存在噪声,利用Lambert-Beer定律难以对其进行准确描述,而利用探测器接受光子的统计特性,能够比较准确地描述X射线的衰减规律.一定数目的入射光子与目标物体相互作用后被探测器接受,对于每个入射光子可以用二项分布来描述.在入射光子数目足够多的情况下,可以用如下的泊松分布模型来描述[10]:

yk~possion(Zk+ek)(1)

Zk=AkI=∑Ni=1akiIi i=1, 2, …, N; k=1,2,...,K(2)

式中,yk为第k个采集单元探测到的投影数据; poisson()为泊松统计模型; A={aki}为投影矩阵,其中aki为第k个采集单元探测到穿过第i个像素的光子概率,Ak表示投影矩阵A的第k个行向量; I为需要求解的目标图像; ek为第k个采集单元的系统误差; N为像素总数; K为采集单元总数; Zk为第k个采集单元探测到的投影数据期望值.

最大后验概率(maximum a posteriori, MAP)算法通过一组投影测量值y,找到图像I,使得条件概率P(I〖JB<1|〗y)最大,计算方法为[11]

P(Ι〖JB<1|〗y)=(P(y〖JB<1|〗I)P(I))/(P(y))(3)

式中,P(I)为图像的先验知识,一般采用Gibbs先验分布函数,其表达式为

P(I)=(exp(-βU(I)))/c(4)

式中,β为正则化参数; c为常数; U(I)为能量函数.

将Gibbs分布应用到MAP算法中,并且通过期望最大化(expectation maximization, EM)算法进行求解,可得MAP-EM算法的迭代求解公式为[12]

Is+1i=(Isi)/(∑kaki+β(∂U(Is))/(∂Ii))∑k(ykaki)/(∑i'aki'Isi')(5)

式中,Isi为目标图像上第i个像素迭代s次后的像素值; β(∂U(Is))/(∂Ii)为惩罚函数项,用于抑制重建图像中的噪声.

1.2 直接迭代基材料分解方法

根据基材料分解模型,可将目标物体的线性衰减系数表示为多种材料的线性衰减系数加权和[13],即

μ(E,l)≈∑nfm-1xm(l)μm(E)(6)

式中,nf为基材料的数目; xm(l)为第m种基材料图像; l为X射线的路径; μm(E)为第m种基材料在能量E下的线性衰减系数.

在能量E下沿路径l的投影函数离散形式可表示为

P(E,l)=∑nfm=1μm(E)Bm(l)(7)

式中,Bm(l)=∫lxm(l)dl 为第m种基材料图像沿路径l的投影.

将管电压u下的连续X射线能谱离散为Q个能量区间,则目标物体在管电压u下的投影可以表示为各能量区间下投影的加权和,即

P(u,l)=∑Qq=1w(Eq)P(Eq,l)=

Qq=1nfm=1w(Eqm(Eq)Bm(l)=∑nfm=1μm(u)Bm(l)(8)

式中,w(Eq)=np(Eq)〖JB<2/〗∑Qq=1np(Eq)为第q个能量区间内光子数的归一化系数,其中np(Eq)为第q个能区的光子数; μm(u)为给定管电压u下某种材料的平均线性衰减系数.

第m种基材料图像沿路径l的投影可以表示为

Bm(l)=Axm(l)m=1,2,…,nf(9)

将式(9)代入式(8)中,可以得到直接迭代投影函数的离散形式为[7]

P=JRx(10)

J=J1CK,J1=[μ1(u1)… μnf(u1)

 

μ1(uNe)… μnf(uNe)](11)

R=CnfA(12)

式中,x=[x1 x2 … xnf]T为各基材料图像的联合矩阵,维度为nfN; J为联合基材料矩阵; Ne为能量区间的数量; CK为K×K的单位矩阵; R为联合投影算子; Cnf为nf×nf的单位矩阵.

对于双能CT双基材料分解问题,即Ne=2, nf =2,其基材料分解函数模型可表示为

P=JRx=[μ1(uL)A μ2(uL)A

μ1(uH)A μ2(uH)A][x1

x2](13)

记W=JR为变换投影算子矩阵,则有P=Wx.由此可见,基材料分解问题可以转化为图像重建问题,根据变换投影算子W以及高、低管电压下的联合投影矩阵P,利用图像重建算法可以得到基材料图像联合矩阵x.利用MAP-EM统计重建算法对上述分解模型进行求解时,记W={ckn},则由式(5)可以得到MAP-EM-DIBMD算法的迭代求解公式为

xs+11i=(xs1i)/(∑kcki+β(∂U(xs1))/(∂x1i))∑k(y1kcki)/(∑i'cki'xs1i')

xs+12i=(xs2i)/(∑kcki+β(∂U(xs2))/(∂x2i))∑k(y2kcki)/(∑i'cki'xs2i')(14)

式中,x1i和 x2i分别为基材料图像 x1 和x2上第i个像素的像素值; U(xs1)和U(xs2)分别为第s次迭代得到的基材料图像xs1和xs2的能量函数; y1k 和 y2k 分别为2种管电压下第k个采集单元探测到的投影数据.设探测器的探测单元数为t,角采样间隔设置为1°,则k的取值范围为1~360t,i'的取值范围为1~t2.

1.3 性能评估指标

利用噪声水平和对比噪声比(contrast-to-noise ratio, CNR)对基材料图像的质量进行评估,其计算公式为

σ=(1/(N-1)∑Ni=1(zi-z-)2)1/2(15)

R=(|z-c-z-b|)/(σb)(16)

式中,σ为噪声水平; R为对比噪声比; zi为第i个像素的像素值; z-为该ROI内所有像素的平均值; z^-c为对比材料区域内像素的均值; z^-b为背景材料区域内像素的均值; σb为背景材料区域的噪声水平.σb越小,R越大,表明图像的噪声水平越低,图像质量越好.

利用均方误差(MSE)对重建图像的精确度进行定量评估,其计算公式为

M=(∑Ni=1(zreci-zrefi)2)/N(17)

式中,zreci和zrefi分别为重建图像和参考图像上第i个像素的像素值.均方误差越小,表示重建图像与参考图像的误差越小,重建图像越接近于参考图像.

利用误差水平[14]对基材料分解结果的准确性进行定量评估,其计算公式为

δ=1/(E2-E1)∫E2E1(|μ(E)-[b1μ1(E)+b2μ2(E)]|)/(μ(E))dE(18)

式中,E1为探测器所能接受的光子最低能量,在考虑过滤材料的情况下,一般将其设置为20 keV; E2为探测器所能接受的光子最高能量,应不小于最大扫描管电压; μ1(E)、μ2(E)、μ(E)分别为2种基材料以及模体在能量E下的线性衰减系数.

1.4 模拟实验装置

采用Geant4(Geometry and Tracking)对快速切换双能CT系统进行模拟.射线源到模体中心的距离设为670 mm,源到探测器的距离设为1 010 mm,X射线源的扇形角设为8.6°.利用Spektr软件包[15]模拟X射线能谱,能量间隔设为1 keV,将管电压分别设为80和140 kV,电压切换周期设置为0.5 ms[16].采用厚度为0.8 mm铍和8.4 mm铝作为过滤材料,以滤除低能X射线,减少低能辐射损伤.

模体设置为圆柱体,其断层剖面结构如图1所示.模体的几何参数以及材料组成见表1.其中,聚乙烯(PE)和羟基磷灰石(HA)分别用于模拟软组织和骨组织,质量分数为10%的生理盐水用于模拟衰减系数介于2种组织衰减系数之间的物质,利用高质量密度的单质材料铝(Al)来测试基材料分解方法的稳定性[7].为了使模拟实验更接近于真实实验,仿真过程中模体材料的线性衰减系数用其管电压下的平均线性衰减系数来描述(见表2).仿真过程中,每次发射的光子数设为2×106,角采样间隔为1°.

图1 圆柱模体剖面图

图1 圆柱模体剖面图

表1 模体的几何参数和组成成分

表1 模体的几何参数和组成成分

表2 2种管电压下各材料的平均线性衰减系数 cm-1

表2 2种管电压下各材料的平均线性衰减系数 cm-1

1.5 基材料分解理论值计算

选取软组织和骨组织作为2种基材料,其中骨组织的等效材料为质量密度1 g/cm3的羟基磷灰石,软组织的等效材料为质量密度0.97 g/cm3的聚乙烯材料.利用各材料在高低管电压下的平均线性衰减系数值,计算其理论分解系数值为

[b1

b2]=[μ1L μ2L

μ1H μ2H]-1L

μH](19)

式中,b1和b2分别为软组织和骨组织的分解系数值; μ1L、μ1H分别为80 kV低管电压和140 kV高管电压下软组织的平均线性衰减系数; μ2L、μ2H分别为80 kV低管电压和140 kV高管电压下骨组织的平均线性衰减系数; μL、μH分别为80 kV低管电压和140 kV高管电压下模体的平均线性衰减系数.在20~140 keV能量范围内,模体中各材料的线性衰减系数如图2所示.理论分解结果见表3.由于空气的线性衰减系数近似为0,故将其忽略不计.

图2 模体中各材料的线性衰减系数曲线

图2 模体中各材料的线性衰减系数曲线

表3 模体中各材料的理论分解系数值

表3 模体中各材料的理论分解系数值

2 实验结果与分析2.1 正则化参数评价

MAP-EM算法在迭代过程中加入了惩罚项,相当于在重建过程中施加一定的先验约束条件来达到抑制噪声的目的.迭代过程中,正则化参数β的选择对重建图像质量具有较大影响.为了评价正则化参数对分解结果的影响,引入如下2个指标:①重建图像与参考图像的近似程度,通过均方误差M来确定; ②重建图像的噪声水平,通过各模体材料的对比噪声比R来确定.M越小,重建图像越接近于参考图像; R越大则图像的噪声越小.计算均方误差M时,将理论分解图像作为参考图像; 计算R时,将PE设置为背景材料.

在MAP-EM-DIBMD方法中,对于每一个正则化参数,随着迭代次数的增加,均方误差呈现先减小后增大的趋势,因此,对于每个正则化参数都存在一个最小的MSE.图3给出了最小MSE以及不同模体材料的CNR随正则化参数变化的曲线.由图3(a)可知,随着正则化参数β的增大,最小MSE呈现先减小后增大的趋势,当β在0.1附近时,其达到最小值.图3(b)为软组织基材料图像上各模体材料的CNR随正则化参数变化的曲线,当β在0.2附近时,材料Al、HA和生理盐水的CNR达到最大值.骨组织基材料图像上各模体材料的CNR随正则化参数变化的曲线如图3(c)所示,同样当β在0.2附近时,材料Al、HA和生理盐水的CNR达到最大值.考虑到β=0.2时,其对应的M值与β=0.1时的M值相差不大,因此,在后续实验中将β设置为0.2.

2.2 基材料分解图像重建结果

在MAP-EM-DIBMD方法中,将迭代停止条件设置为2次迭代结果之间的MSE差值小于1×10-6.当β=0.2时,MAP-EM-DIBMD方法得到的基材料分解图像见图4(a).图4(b)为FBP-IDBMD方法得到的基材料分解图像,图4(c)为理论方法计算得到的基材料分解图像.

图3 最小均方误差及对比噪声比随正则化参数变化的曲线

图3 最小均方误差及对比噪声比随正则化参数变化的曲线

图4 不同基材料分解方法得到的基材料图像

图4 不同基材料分解方法得到的基材料图像

2.3 结果分析2.3.1 噪声水平分析

表4给出了基于FBP-IDBMD 和MAP-EM-DIBMD两种方法得到的基材料图像的噪声水平.由表可知,与FBP-IDBMD方法相比,MAP-EM-DIBMD方法可使基材料图像的噪声水平下降57.42%~63.64%.

表4 基材料图像的噪声水平

表4 基材料图像的噪声水平

2.3.2 对比噪声比分析

表5给出了基于FBP-IDBMD 和MAP-EM-DIBMD两种方法得到的模体中各材料区域(用作计算CNR的背景材料PE除外)的对比噪声比.表中,R1和R2分别表示FBP-IDBMD和MAP-EM-DIBMD方法得到的对比噪声比.由表可知,与FBP-IDBMD方法相比,MAP-EM-DIBMD方法可使对比噪声比提高1.37%~223.17%.

表5 模体中各材料区域的对比噪声比

表5 模体中各材料区域的对比噪声比

2.3.3 误差水平分析

计算MAP-EM-DIBMD和FBP-IDBMD两种方法得到的不同材料区域2种基材料分解系数b1和b2的均值,并对其误差水平进行分析,各材料区域分解系数均值见表6,分解误差水平见图5.由表可知,与FBP-IDBMD方法相比,MAP-EM-DIBMD方法可使Al、HA、生理盐水、PE的误差水平分别下降62.14%、31.72%、45.75%、45.73%.

表6 模体中各材料的分解系数均值

表6 模体中各材料的分解系数均值

图5 2种方法的δ值比较

图5 2种方法的δ值比较

3 结论

1)针对双能CT,提出了基于MAP-EM算法的直接迭代基材料分解方法, 将基材料分解问题转化为基材料图像重建问题,推导出MAP-EM-DIBMD方法的迭代求解公式,有效提高基材料分解的精度和鲁棒性.

2)仿真实验结果表明,MAP-EM-DIBMD算法显著降低了基材料分解图像的噪声水平和误差水平, 提高了基材料分解图像的对比噪声比,具有重要的实际应用价值.

3)MAP-EM-DIBMD算法采用迭代的方法重建分解图像,耗时较长,可采用有序子集的方法对算法进行加速,在实际应用中还可采用GPU提高算法的计算效率.下一步可将所提方法拓展到光子计数能谱CT并应用于临床研究,对其性能进行进一步评估和改进,为临床精确诊断提供服务.

参考文献