不完全投影下的CT图像重建仿真研究

李红艳 万钟林

摘 要：低剂量辐射是CT成像中的研究热点和难点，而不完全角度投影又是低剂量辐射中的一种有效手段。因此对不完全角度投影的重建问题也是一个研究热点。建立在CS理论基础上的TV-EM算法虽然能够克服传统算法的重建精度不高，去噪能力不明显等特点，但也对边缘进行了过度平滑。本文在惩罚正则项中加入中值先验分布，在抑制噪声的同时更能保留边缘信息。加入中值先验，形成新的TV-EM算法，并将其应用于不完全投影下的重建。通过仿真实验，选择合适的惩罚因子，验证了其超强的抗噪能力和稳健的收敛性。

关键词：不完全投影重建;稀疏惩罚;贝叶斯方法

中图分类号：TP391.9  文献标识码：A  文章编号：1673-260X（2021）02-0006-04

1 引言

随着科技的发展，CT在临床诊断中的作用日显重要。而人们对健康的重视程度也越来越高。因此CT在医学诊断中的扫描所产生的射线辐射也受到人们的广泛关注。毕竟，射线辐射对人体健康是有损害的。据报道，接受超过28次CT扫描的患者致癌几率比平均水平高12%，而儿童在接受CT检查时所受的辐射影响会更大。因此，人们都希望接受足够小的辐射，而能得到精确的CT图像辅助诊断。因此，低剂量辐射也成为医学影像中亟待解决的热点和难点。而有限角度和稀疏角度重建就是降低辐射剂量的一种有效手段。

在不完全角度投影数据条件下，无论是有限角度还是稀疏角度，经典FBP算法的重建结果总是不够理想，虽然重建速度快，但是存在伪影，甚至有一部分信息严重缺失，导致重建图像质量不尽如人意[1]。这就推动了迭代重建算法的研究。2006年，Candes等[2]提出了压缩感知（CS）理论。随后做了大量重建算法的研究，如Sidky等人提出全变分最小化算法[3，4]，Chen等[5]于2008年提出PICCS算法，而Wang等[6]于2010年提出RRD算法，Zhang等[7]于2013年提出ADTVM算法，都是基于CS理论，选取一种或几种稀疏变换作为目标函数，再以成像模型AX=P作为约束条件，再求解出重建图像[8]。在实际成像过程中，由于数据的随机噪声和成像模型的非随机误差，造成成像模型AX=P总是不相容的。所以用成像模型AX=P作为约束条件并不总是最合适的。本文研究的TV-EM（全变差-期望最大）算法[9]，是用贝叶斯方法建立目标函数，以ML-EM算法作为基础，以图像X的稀疏变换提取图像的先验知识，加入图像的中值先验到每一步迭代中，从而改善重建图像的质量[10，11]。

2 获取投影

重建物体的离散模型以网格像素（二维图像）的形式进行刻画，获取投影数据模型如图1所示[12]。探测器检测结果为射线方向的累积值，并假设射线没有宽度。探测器检测的投影值与射线穿过体素的长度相关。实际上投影值是用像素内的线段长度加权的像素值的和。

这里，像素值xj为第j个像素的灰度值。探测器的第i个探测元发出的射线在每个像素内的线段长度记为aij。如果第i条射线与第j个像素不相交，则aij=0。而第i条射线的投影值pi是探测到的光子数。而投影方程可表示为：

pi=aijxj

3 贝叶斯图像重建原理

对于一个成像问题来说，假定每一个图像像素发出的光子数是一个独立的泊松随机变量。而实际上第i个探测元理论上探测到像素j发出的光子数的均值是aijxj。而所有的探测元探测到的pi组成一个样本。那么可建立似然函数[17]：

prob=∏i（∑jaijxj）  （2-1）

取对数似然

ln（prob）=∑i[-∑jaijxj+piln（∑jaijxj）-lnpi！] （2-2）

在图像重建中，ML-EM方法就是寻找一个解，使得探测到的数据最有可能发生，即最大化概率Prob，其中pi为探测到的投影数据，而xj是待求图像。可用如下优化目标来表示：

=argmaxx≥0ln（prob）  （2-3）

为求解目标（2-3），Shepp和Vardi提出ML-EM算法，迭代公式[16]为：

j=∑i  （2-4）

其中pi中是投影数据，xj是图像当前估计值。

4 TV-EM算法

TV约束是基于Donho，Candes和Tao等提出的压缩感知理论，假定图像具有分段光滑性质。选择图像具有稀疏性的TV范数作为正则项，加入（2-3）作为正则约束。

根据Candes和Donoho提出的理论，图像重建目标变为：

=argmaxx≥0ln（prob）+？琢TV（x）  （3-1）

？琢——正则化参数

TV（x）——圖像的TV范数

TV（全变差）范数是图像的梯度的l1范数，但由于使用起来不方便，习惯用l2范数近似代替l1范数，这样就可以对TV范数求导了。但是l2范数在抑制噪声的同时，把边缘信息也过度平滑了。引入TV范数为稀疏惩罚项的TV-EM算法迭代公式[9]为：

j=∑i  （3-2）

其中，TV范数求偏导是离散化后的公式[12]：

=

+

-  （3-3）

参数？着是一个很小的数，为了避免分母为0。

5 中值先验

中值先验在降低噪声的同时，也保持了图像的边缘，它的分布如下[18]：

R（x）=aexp-∑j

a——常数;

？茁——正则化参数;

Mj——领域中像素的中间值。

将（4-1）取对数后加入目标函数（3-1）中，可得

=argmaxx≥0ln（prob）+？琢TV（x）+lnR（x） （4-2）

那么加入中值先验的新的迭代公式如下：

j=∑i （4-3）

不完全投影数据的背景下，每一步迭代中都需计算当前的梯度偏导数（3-3）及中值先验值，综合了稀疏先验和中值先验的方法，希望能在稀疏约束下抑制去躁的同时能够保持边缘的信息。

基于中值先验的TV-EM算法的具体步骤如下：

Step 1.给未知量xj赋初值，xj（1）=xj（0），j=1，2，…，M，t=1。

Step 2.计算所有投影的估计值i=aijxj（i），i=1，2，…，I。

Step 3.计算误差，？驻i=，i=1，2，…，I。

Step 4.计算当前估计值的TV范数偏导数和中值先验。

Step 5.按照（4-3）获取新的像素值，完成一轮迭代。

Step 6.赋值：t=t+1，然后重复（2）到（5）進行新一轮迭代。

6 仿真实验

6.1 重建图像评判准则[19]

本文采用128×128的Shepp-Logan头模型，应用MATLAB7.0模拟获取CT的投影数据。分别根据FBP、ML-EM、TV-EM和本文方法进行重建。并把原始的图像和重建得到的图像进行比较，就可以判断出重建结果的优劣。本文使用PSNR值来进行图像质量的比较，其公式如下：

PSNR=10×log 10

=20×log 10

其中，MSE是原图像与处理图像之间均方误差。

MAXI表示图像颜色的最大数值，8位采样点表示为255。

PSNR值越大，说明重建图像与原图越相似，重建结果越好。

6.2 仿真实验

本文采用128×128的Shepp-Logan头模型，通过MATLAB R2016模拟获取投影数据。分别取稀疏均匀角度20个和30个。并且为了观察算法的去噪效果，我们在投影数据中加入0.05的泊松噪声。参数？着=1e-7。

其中正则化参数取为α=0.11，β=0.01。如图2所示。图2中（a）原图;（b）20个角度下FBP无噪声重建;图（c-e）20个角度投影算法100次迭代的结果;（f-h）30个角度投影各算法100次迭代的结果。依直观感受，在稀疏均匀角度下，经典FBP算法重建的图像存在严重伪影，而且去噪也不理想，而ML-EM算法的去噪能力在稀疏均匀角度下也表现一般，而TV-EM算法的重建效果稳定，去噪能力较强，本文方法重建的结果不仅去躁明显，边缘也更加清晰。

各方法重建结果的PSNR值如表1所示。数据显示，本文方法重建结果的PSNR值较大，更接近原图，结果更好。

7 结语

在传统的ML-EM算法中加入稀疏惩罚项，得到的TV-EM算法虽然能够克服传统算法的重建精度不高，去噪能力不明显等特点，但也对边缘进行了过度平滑。本文在惩罚正则项中加入中值先验分布，在抑制噪声的同时更能保留边缘信息。但是方法的结果依赖于参数的选择。对于参数现在只能依靠经验，所以这是下一步研究的工作。另外，重建速度还有待加快。

——————————

参考文献：

〔1〕蔺鲁萍，王永革.不完全角度CT图像重建的模型与算法[J].北京航空航天大学学报，2017，43（04）：823-830.

〔2〕Donoho D L. Compressed sensing. Information Theory[J]，IEEE Transactions on，2006，52（04）：1289-1306.

〔3〕Sidky E Y， Kao C M， Pan X. Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT [J]. Journal of X-ray Science and Technology， 2006， 14（02）： 119-139.

〔4〕Sidky E Y， Pan X. Image reconstruction in circular cone-beam computed tomography by constrained， total-variation minimization [J].Physics in Medicine and Biology， 2008， 53（17）： 4777-807.

〔5〕Chen G H，Tang J，and Leng S.Prior image consrained compressed sensing（PICCS）：a method to accurately reconstruct dynamic CT images from highly undersampled projection data sets[J]，Med.Phys.， 2008， 35（02）：660-663.

〔6〕Wang L，Li L，Yan B，et al.An algorithm for CT Image Reconstruction from Limited-View Projections[J].Chinese Physics B， 2010， 19（08）：088-106.

〔7〕Zhang H，Wang L，Yan B，et al.Image reconstruction based on alternating direction Total Variation minimization in linear scan Computed Tomography[J].Chinese Physics B，2013， 22（07）：078-107.

〔8〕Bian J，Siewerdsen JH，Han X，et al. Evaluation of Sparse-view Reconstruction from Flat-panel-detector Cone-beam CT[J].Physics in Medicine&Biology. 2011，55（22）： 6575-6599.

〔9〕曾更生.医学图像重建[M].第1版.北京：高等教育出版社，2010.148-156.

〔10〕Buades A，Coll B，Morel JM.A Review of Image Denoising Algorithms，with a New one[J].Siam Journal on Multiscale Modeling&Simulation.2005，4（02）：490-530.

〔11〕Buades A，Coll B，Morel JM.A Non-Local Algorithm for Image Demoising[C].IEEE Computer Society Conference on Computer Vision&Pattern Re ognition.IEEE Computer Society， 2005.60-65.

〔12〕陳建林，闫镔，李磊，等.CT重建中投影矩阵模型研究综述[J].CT理论与应用研究，2014，23（02）：317-328.

〔13〕Wang T，Kievit FM，Veiseh O，et al.Limited-angle cone-beam computed tomography image reconstruction by total variation minimization and piecewise-constant modification[J].Journal of Inverse and Ill-posed Problems.2016， 21（06）：735-754.

〔14〕Natterer F，Wubbeling F.Mathematical Methods in Image Reconstruction[M]. Society for Industrial and Applied Mathematics， 2001.107-108.

〔15〕张蕾，宋文爱，杨顺明.超声散射CT技术发展综述[J].CT理论与应用研究，2011，20（03）：415-424.

〔16〕闫镔，李磊.CT图像重建算法[M].第1版.北京：科学出版社，2017.97-109.

〔17〕廖文熙.PET图像重建算法的研究与优化[D].浙江大学，2010.

〔18〕何玲君.基于混合先验的CT图像重建研究[J].传感器世界，2011，17（03）：16-19.

〔19〕蒋刚毅，黄大江，王旭.图像质量评价方法研究进展[J].电子与信息学报，2010，32（01）：219-226.