东北大学 中荷生物医学与信息工程学院, 辽宁 沈阳 110169
收稿日期:2016-06-08
基金项目:国家自然科学基金资助项目(61372014)。
作者简介:侯晓文(1989-), 男, 山东菏泽人, 东北大学博士研究生;
康雁(1964-), 男, 辽宁沈阳人, 东北大学教授,博士生导师。
摘要:针对重建算法对不同的检测对象需要建立不同的查找表问题, 提出了一种基于最大似然-可分离抛物面型替代函数的重建算法.根据双能CT的物理模型和统计模型建立了对数似然函数, 并以之为目标函数.根据目标函数的凸性, 构造了可分离抛物面型替代函数.实验结果表明, 该算法重建所得各能级图像与原始图像的相关系数大于0.983, 信噪比大于12 dB, 均大于查表法重建结果的相应值, 重建图像质量优于查表法.
关键词:双能CT最大似然凸优化替代函数重建算法
Maximum Likelihood-Separable Paraboloidal Surrogate Function Algorithm for Dual-Energy CT Reconstruction
HOU Xiao-wen, TENG Yue-yang, LIU Yu-jia, KANG Yan
School of Sino-Dutch Biomedical and Information Engineering, Northeastern University, Shenyang 110169, China
Corresponding author: KANG Yan, professor, E-mail: kangyan@bmie.neu.edu.cn
Abstract: In the reconstruction algorithm, different lookup table problems need to be set up for different objects. An algorithm based on maximum likelihood-separable paraboloidal surrogate function was proposed. A log-likelihood function, as the objective function, was constructed based on the physical model and statistical model of dual-energy CT. Separable paraboloidal surrogate function was constructed according to the convex characteristic of the objective function. The experiment result shows that, for all of the energy levels, the correlation coefficient between the image reconstructed by this method and the original image is greater than 0.983, and the signal to noise ratio is greater than 12 dB, both are greater than the corresponding values of the result using the look up table algorithm. The quality of image reconstructed by the proposed method was better than the look up table method.
Key Words: dual-energy CTmaximum likelihoodconvex optimizationsurrogate functionreconstruction algorithm
计算机断层成像(computed tomography, CT)技术应用十分广泛[1].但尚存在一些不足, 比如金属伪影[2].双能CT(dual-energy CT, DECT)能够显著降低金属伪影, 并且能够分辨被检测目标的组成成分, 大大改善CT的检测质量[3-6].
DECT概念是1976年由Alvarez和Macovski提出的.物质对X射线的吸收系数是随X射线的能量变化的, 对于大多数物质而言, 这种变化的衰减系数可以由光电子效应和康普顿散射的线性组合精确近似.DECT通过两束能谱不同的X射线扫描目标物体, 得到两套投影数据.通过相应的重建算法, 得到光电子效应和康普顿散射的系数图.自DECT的概念提出后, 国内外研究人员对其重建算法进行了大量研究, 随能级变化吸收系数的线性组合也由原来的光电效应和康普顿散射转换为基物质分解[7].目前国内的双能重建算法主要是查表法(look up table, LUT)[3].查表法原理简单, 易于实现.但是, 为建立查找表, 需要事先知道目标物体的组成成分.此外, 对于不同的检测对象, 需要建立不同的查找表.因而, 实用性不太强.为此, 本文提出一种基于最大似然-可分离抛物面替代函数(maximum likelihood-separable paraboloidal surrogate, ML-SPS)的迭代重建算法.
本文首先介绍ML-SPS算法, 然后利用该算法对GATE得到的模拟数据投影进行重建实验, 并将结果和查表法的重建结果对比, 最后根据实验结果得出结论.
1 模型和算法1.1 双能CT的物理模型和统计模型被扫描物体对X射线的衰减系数随X射线能级大小的变化而改变.据此, Alvarez等提出了精确X射线衰减物理模型的数学表达式[2]:
(1) |
DECT重建中, 需要两个具有不同能谱的X射线扫描目标物体, 得到两套投影数据.一束X射线的能谱所包含的能级数远大于2.为了能够用两套投影数据重建DECT, 需要对衰减系数进行分解.常用的分解方法是双效应分解和双物质分解.本文采用双物质分解, 并取空气和碘作为基物质.如此, 任意一种物质的衰减系数都可以表示为空气和碘的衰减系数的线性组合, 即
(2) |
(3) |
(4) |
(5) |
(6) |
第i条射线穿过目标物体后测量得到投影值yi的概率为
(7) |
(8) |
(9) |
(10) |
(11) |
(12) |
1.3 凸优化方法简化目标函数式(10)所示的目标函数难以直接求解.为简化目标函数的优化过程, 需要对目标函数进行优化转换.目标函数中存在负指数和负对数的形式, 因此考虑采用凸优化的方法对目标函数进行简化[8].
首先需要构建凸组合系数.凸组合系数需要满足两个条件[8]: ①所有系数之和为1;②每个系数都是非负的.
目标函数(10)中既含有指数项又含有对数项, 因此, 对两个部分分别构建两个不同的凸组合系数.对于指数项部分, 构建凸组合系数:
(13) |
(14) |
根据凸函数的性质可得
(15) |
构建凸组合系数:
(16) |
(17) |
利用凸函数的性质可以得到指数项的替代函数:
(18) |
(19) |
1.4 泰勒级数展开进一步简化目标函数替代函数相对于原来的函数需要满足三条性质, 才能保证利用替代函数求得的极值点能够保证原函数收敛.对于函数h(x), 假设在xk处构建其替代函数h′(x, xk), 这三条性质可表示为
(20) |
(21) |
(22) |
(23) |
把替代函数Q2和Q3代入到Q1可以得到可分离的抛物面型替代函数(SPS).对替代函数各个自变量求偏导数, 并令偏导数等于零得到迭代公式:
(24) |
(25) |
(26) |
2) repeat
For j=1, n do:
式(24);
c′j=t*cj;
式(25);
d′j=t*dj;
End for
3) C:=C’; D:= D’;
4) until{stopping criteria}
2 模拟实验结果与分析为了验证ML-SPS算法的重建效果, 利用GATE模拟光子计数型探测器CT扫描的过程, 得到两组投影数据, 并利用这两组数据进行了重建.
2.1 模拟实验参数GATE模拟中, 两束入射X射线的电压分别为80 kVp和140 kVp, 其能谱范围分别为40~80 keV和40~140 keV.被扫描物体是一个有6个等大的圆形区域的方形模体, 如图 1所示.
图 1(Fig. 1)
图 1 模体Fig.1 Phantom |
模体的6个圆形区域自左至右, 从上到下, 依次填充的是水、碳、钙、铝、硫和氯化钠.背景区域填充空气.模体边长为128 mm.
GATE模拟CT扫描过程的几何参数:射线源到模体中心的距离为550 mm, 探测器为直线型光子计数探测器, 通道数为128, 探测器到模体中心距离为90 mm, 共扫描360个角度, 角度间隔为1°.
2.2 重建结果分析针对GATE模拟得到的投影数据分别进行了ML-SPS算法重建和LUT算法重建.其中, LUT算法建立的查找表的精度为10-5, ML-SPS算法中投影矩阵的计算采用Siddon的算法[9], 迭代次数设置为200次.图 2展示了两种方法重建得到的10 keV能级下的图像.
图 2(Fig. 2)
图 2 10 keV图像Fig.2 Image of 10 keV |
对比图 1和图 2可以看出,ML-SPS和LUT算法得到的重建结果与原始图像差别不大, 且两种算法的重建结果相似.
为了更为客观地对比分析ML-SPS算法的重建结果, 分别计算了ML-SPS和LUT算法相对于原始图像的相关系数(correlation coefficient, CC)和信噪比(signal noise ratio, SNR)[10].由于DECT重建结果中包含多个能级的图像(比如文中能谱包含40~140 keV共101个能级的图像), 因此, 需要对每个能级的图像分别求取CC和SNR.为了便于分析所有能级的CC和SNR值, 以能级为横坐标, CC值和SNR值为纵坐标, 画出了CC值和SNR值曲线,如图 3和图 4所示.
图 3(Fig. 3)
图 3 相关系数曲线Fig.3 Curves of correlation coefficient |
图 4(Fig. 4)
图 4 信噪比曲线Fig.4 Curves of SNR |
由图 3和图 4可见, 对于所有能级, ML-SPS算法重建图像的信噪比和相关系数均略大于LUT算法.ML-SPS算法重建质量略优于LUT算法.
3 结论1) 根据DECT物理模型和统计模型的特征, 得出ML目标函数, 并利用其凸函数性质, 通过凸优化和拟泰勒级数展开的方法构建了SPS, 将原目标函数的优化问题转换为线性问题, 并证明了其收敛性.
2) 所提出的ML-SPS算法模拟实验结果略优于传统的LUT算法.
3) 所提出的ML-SPS算法弥补了LUT算法需要事先根据被扫描物体的成分建立查找表的不足, 更具有实用性.
参考文献
[1] | Ying Z, Naidu R, Crawford C R. Dual energy computed tomography for explosive detection[J].Journal of X-Ray Science and Technology, 2006, 14(4): 235–256. |
[2] | Alvarez R E, Macovski A. Energy-selective reconstructions in X-ray computerized tomography[J].Medical Physics, 1976, 21(5): 733–744. |
[3] | 高洋. 双能CT图像重建算法研究[D]. 重庆: 重庆大学, 2012. ( Gao Yang.Study on image reconstruction algorithm of dual-energy CT[D]. Chongqing:Chongqing University, 2012.http://cdmd.cnki.com.cn/Article/CDMD-10611-1012050080.htm) |
[4] | Zhang R Q, Thibault J B, Bouman C A, et al. Model-based iterative reconstruction for dual-energy X-ray CT using a joint quadratic likelihood model[J].IEEE Transaction on Medical Imaging, 2014, 33(1): 117–134.DOI:10.1109/TMI.2013.2282370 |
[5] | Fessler J A, Elbakri I, Sukovic P, et al.Maximum likelihood dual-energy tomographic image reconstruction[C]// Medical Imaging 2002:Image Processing.San Diego:SPIE, 2002:38-49. |
[6] | Long Y, Fessler J A. Multi-material decomposition using statistical image reconstruction for spectral CT[J].IEEE Transaction on Medical Imaging, 2014, 33(8): 1614–1626.DOI:10.1109/TMI.2014.2320284 |
[7] | Wang A S, Hsieh S S, Pelc N J. A review of dual energy CT:principles, applications, and future outlook[J].CT Theory and Applications, 2012, 21(3): 367–386. |
[8] | De Pierro A R. A modified expectation maximization algorithm for penalized likelihood estimation in emission tomography[J].IEEE Transaction on Medical Imaging, 1995, 14(1): 132–137.DOI:10.1109/42.370409 |
[9] | Siddon R. Fast calculation of the exact radiological path for a three-dimensional CT array[J].Medical Physics, 1985, 12(2): 252–255.DOI:10.1118/1.595715 |
[10] | 何玲君. 基于最大似然和罚似然估计的CT统计重建算法研究[D]. 太原: 中北大学, 2011. ( He Ling-jun.The study of CT statistical reconstruction algorithm based on maximum likelihood and penalized likelihood estimates[D]. Taiyuan:North University of China, 2011.http://cdmd.cnki.com.cn/Article/CDMD-10110-1011156236.htm) |