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

李紅艷 萬鐘林

摘 要：低劑量輻射是CT成像中的研究熱點和難點，而不完全角度投影又是低劑量輻射中的一種有效手段。因此對不完全角度投影的重建問題也是一個研究熱點。建立在CS理論基礎(chǔ)上的TV-EM算法雖然能夠克服傳統(tǒng)算法的重建精度不高，去噪能力不明顯等特點，但也對邊緣進(jìn)行了過度平滑。本文在懲罰正則項中加入中值先驗分布，在抑制噪聲的同時更能保留邊緣信息。加入中值先驗，形成新的TV-EM算法，并將其應(yīng)用于不完全投影下的重建。通過仿真實驗，選擇合適的懲罰因子，驗證了其超強(qiáng)的抗噪能力和穩(wěn)健的收斂性。

關(guān)鍵詞：不完全投影重建;稀疏懲罰;貝葉斯方法

中圖分類號：TP391.9? 文獻(xiàn)標(biāo)識碼：A? 文章編號：1673-260X（2021）02-0006-04

1 引言

隨著科技的發(fā)展，CT在臨床診斷中的作用日顯重要。而人們對健康的重視程度也越來越高。因此CT在醫(yī)學(xué)診斷中的掃描所產(chǎn)生的射線輻射也受到人們的廣泛關(guān)注。畢竟，射線輻射對人體健康是有損害的。據(jù)報道，接受超過28次CT掃描的患者致癌幾率比平均水平高12%，而兒童在接受CT檢查時所受的輻射影響會更大。因此，人們都希望接受足夠小的輻射，而能得到精確的CT圖像輔助診斷。因此，低劑量輻射也成為醫(yī)學(xué)影像中亟待解決的熱點和難點。而有限角度和稀疏角度重建就是降低輻射劑量的一種有效手段。

在不完全角度投影數(shù)據(jù)條件下，無論是有限角度還是稀疏角度，經(jīng)典FBP算法的重建結(jié)果總是不夠理想，雖然重建速度快，但是存在偽影，甚至有一部分信息嚴(yán)重缺失，導(dǎo)致重建圖像質(zhì)量不盡如人意[1]。這就推動了迭代重建算法的研究。2006年，Candes等[2]提出了壓縮感知（CS）理論。隨后做了大量重建算法的研究，如Sidky等人提出全變分最小化算法[3，4]，Chen等[5]于2008年提出PICCS算法，而Wang等[6]于2010年提出RRD算法，Zhang等[7]于2013年提出ADTVM算法，都是基于CS理論，選取一種或幾種稀疏變換作為目標(biāo)函數(shù)，再以成像模型AX=P作為約束條件，再求解出重建圖像[8]。在實際成像過程中，由于數(shù)據(jù)的隨機(jī)噪聲和成像模型的非隨機(jī)誤差，造成成像模型AX=P總是不相容的。所以用成像模型AX=P作為約束條件并不總是最合適的。本文研究的TV-EM（全變差-期望最大）算法[9]，是用貝葉斯方法建立目標(biāo)函數(shù)，以ML-EM算法作為基礎(chǔ)，以圖像X的稀疏變換提取圖像的先驗知識，加入圖像的中值先驗到每一步迭代中，從而改善重建圖像的質(zhì)量[10，11]。

2 獲取投影

重建物體的離散模型以網(wǎng)格像素（二維圖像）的形式進(jìn)行刻畫，獲取投影數(shù)據(jù)模型如圖1所示[12]。探測器檢測結(jié)果為射線方向的累積值，并假設(shè)射線沒有寬度。探測器檢測的投影值與射線穿過體素的長度相關(guān)。實際上投影值是用像素內(nèi)的線段長度加權(quán)的像素值的和。

這里，像素值xj為第j個像素的灰度值。探測器的第i個探測元發(fā)出的射線在每個像素內(nèi)的線段長度記為aij。如果第i條射線與第j個像素不相交，則aij=0。而第i條射線的投影值pi是探測到的光子數(shù)。而投影方程可表示為：

pi=aijxj

3 貝葉斯圖像重建原理

對于一個成像問題來說，假定每一個圖像像素發(fā)出的光子數(shù)是一個獨(dú)立的泊松隨機(jī)變量。而實際上第i個探測元理論上探測到像素j發(fā)出的光子數(shù)的均值是aijxj。而所有的探測元探測到的pi組成一個樣本。那么可建立似然函數(shù)[17]：

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

取對數(shù)似然

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

在圖像重建中，ML-EM方法就是尋找一個解，使得探測到的數(shù)據(jù)最有可能發(fā)生，即最大化概率Prob，其中pi為探測到的投影數(shù)據(jù)，而xj是待求圖像?？捎萌缦聝?yōu)化目標(biāo)來表示：

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

為求解目標(biāo)（2-3），Shepp和Vardi提出ML-EM算法，迭代公式[16]為：

j=∑i? （2-4）

其中pi中是投影數(shù)據(jù)，xj是圖像當(dāng)前估計值。

4 TV-EM算法

TV約束是基于Donho，Candes和Tao等提出的壓縮感知理論，假定圖像具有分段光滑性質(zhì)。選擇圖像具有稀疏性的TV范數(shù)作為正則項，加入（2-3）作為正則約束。

根據(jù)Candes和Donoho提出的理論，圖像重建目標(biāo)變?yōu)椋?/p>

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

？琢——正則化參數(shù)

TV（x）——圖像的TV范數(shù)

TV（全變差）范數(shù)是圖像的梯度的l1范數(shù)，但由于使用起來不方便，習(xí)慣用l2范數(shù)近似代替l1范數(shù)，這樣就可以對TV范數(shù)求導(dǎo)了。但是l2范數(shù)在抑制噪聲的同時，把邊緣信息也過度平滑了。引入TV范數(shù)為稀疏懲罰項的TV-EM算法迭代公式[9]為：

j=∑i? （3-2）

其中，TV范數(shù)求偏導(dǎo)是離散化后的公式[12]：

=

+

-? （3-3）

參數(shù)？著是一個很小的數(shù)，為了避免分母為0。

5 中值先驗

中值先驗在降低噪聲的同時，也保持了圖像的邊緣，它的分布如下[18]：

R（x）=aexp-∑j

a——常數(shù);

？茁——正則化參數(shù);

Mj——領(lǐng)域中像素的中間值。

將（4-1）取對數(shù)后加入目標(biāo)函數(shù)（3-1）中，可得

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

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

j=∑i （4-3）

不完全投影數(shù)據(jù)的背景下，每一步迭代中都需計算當(dāng)前的梯度偏導(dǎo)數(shù)（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.計算當(dāng)前估計值的TV范數(shù)偏導(dǎo)數(shù)和中值先驗。

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

Step 6.賦值：t=t+1，然后重復(fù)（2）到（5）進(jìn)行新一輪迭代。

6 仿真實驗

6.1 重建圖像評判準(zhǔn)則[19]

本文采用128×128的Shepp-Logan頭模型，應(yīng)用MATLAB7.0模擬獲取CT的投影數(shù)據(jù)。分別根據(jù)FBP、ML-EM、TV-EM和本文方法進(jìn)行重建。并把原始的圖像和重建得到的圖像進(jìn)行比較，就可以判斷出重建結(jié)果的優(yōu)劣。本文使用PSNR值來進(jìn)行圖像質(zhì)量的比較，其公式如下：

PSNR=10×log 10

=20×log 10

其中，MSE是原圖像與處理圖像之間均方誤差。

MAXI表示圖像顏色的最大數(shù)值，8位采樣點表示為255。

PSNR值越大，說明重建圖像與原圖越相似，重建結(jié)果越好。

6.2 仿真實驗

本文采用128×128的Shepp-Logan頭模型，通過MATLAB R2016模擬獲取投影數(shù)據(jù)。分別取稀疏均勻角度20個和30個。并且為了觀察算法的去噪效果，我們在投影數(shù)據(jù)中加入0.05的泊松噪聲。參數(shù)？著=1e-7。

其中正則化參數(shù)取為α=0.11，β=0.01。如圖2所示。圖2中（a）原圖;（b）20個角度下FBP無噪聲重建;圖（c-e）20個角度投影算法100次迭代的結(jié)果;（f-h）30個角度投影各算法100次迭代的結(jié)果。依直觀感受，在稀疏均勻角度下，經(jīng)典FBP算法重建的圖像存在嚴(yán)重偽影，而且去噪也不理想，而ML-EM算法的去噪能力在稀疏均勻角度下也表現(xiàn)一般，而TV-EM算法的重建效果穩(wěn)定，去噪能力較強(qiáng)，本文方法重建的結(jié)果不僅去躁明顯，邊緣也更加清晰。

各方法重建結(jié)果的PSNR值如表1所示。數(shù)據(jù)顯示，本文方法重建結(jié)果的PSNR值較大，更接近原圖，結(jié)果更好。

7 結(jié)語

在傳統(tǒng)的ML-EM算法中加入稀疏懲罰項，得到的TV-EM算法雖然能夠克服傳統(tǒng)算法的重建精度不高，去噪能力不明顯等特點，但也對邊緣進(jìn)行了過度平滑。本文在懲罰正則項中加入中值先驗分布，在抑制噪聲的同時更能保留邊緣信息。但是方法的結(jié)果依賴于參數(shù)的選擇。對于參數(shù)現(xiàn)在只能依靠經(jīng)驗，所以這是下一步研究的工作。另外，重建速度還有待加快。

——————————

參考文獻(xiàn)：

〔1〕藺魯萍，王永革.不完全角度CT圖像重建的模型與算法[J].北京航空航天大學(xué)學(xué)報，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〕曾更生.醫(yī)學(xué)圖像重建[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理論與應(yīng)用研究，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技術(shù)發(fā)展綜述[J].CT理論與應(yīng)用研究，2011，20（03）：415-424.

〔16〕閆鑌，李磊.CT圖像重建算法[M].第1版.北京：科學(xué)出版社，2017.97-109.

〔17〕廖文熙.PET圖像重建算法的研究與優(yōu)化[D].浙江大學(xué)，2010.

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

〔19〕蔣剛毅，黃大江，王旭.圖像質(zhì)量評價方法研究進(jìn)展[J].電子與信息學(xué)報，2010，32（01）：219-226.