基于導引圖像濾波和截斷全變分的CT重建算法

時鴻雁， 吳愛弟， 冀東江

(天津職業(yè)技術師范大學 理學院， 天津 300222)

0 引言

根據(jù)采集到的投影數(shù)據(jù)，利用重建算法重建出待檢測樣品的內部結構，這項技術稱為計算機斷層成像技術(Computed Tomography，CT)[1]。解析重建算法和迭代重建算法是兩種常用的圖像重建算法。其中濾波反投影重建算法(Filtered Back Projection, FBP)[1]屬于解析重建算法，代數(shù)重建算法(Algebraic Reconstruction Technique, ART)[3]和聯(lián)合代數(shù)重建算法(Simultaneous Algebraic Reconstruction Technique, SART)[4]屬于兩種常見的迭代重建算法。在投影數(shù)據(jù)不完整或噪聲較大的情況下，與解析重建算法相比，迭代重建方法可以重建出更好的結果。迭代重建算法的優(yōu)勢是可以將先驗信息融合到重建圖像過程當中，在實際應用過程中，常常采用正則化技術來引入先驗信息以提高重建圖像的質量。潘曉川[6]等將TV作為正則項引入到ART算法中，重建后得到的圖像質量更好，但是對于具有更多圖像紋理細節(jié)和復雜邊緣結構的圖像，將導致重建的圖像過度平滑。后來Dou[7]等提出了截斷全變分用于圖像平滑，截斷TV去除不重要的細節(jié)的同時還可以保留顯著的邊緣，而且它僅對圖像的部分梯度進行懲罰，由此可以解決TV對較大的梯度幅度進行懲罰而導致的過度平滑的問題。因此，我們可以將截斷TV作為正則項用于圖像重建。

何凱明等[8]提出了導引圖像濾波，它能更好的利用導引圖像的結構，把導引圖像含有的重要特征轉移到輸入的圖像中，并讓濾波后的圖像保留輸入圖像的部分信息。冀東江等[9]提出了基于導引圖像濾波的聯(lián)合代數(shù)重建算法，在圖像重建過程中更新先驗指導圖像，迭代地合并信息，提高了重建圖像的質量。余維等[10]提出了一種稀疏性誘導的動態(tài)引導圖像濾波算法。該算法利用引導圖像濾波，把SART算法重建的圖像作為輸入圖像，TDM-STF算法重建的圖像作為導引圖像，以解決有限稀疏角CT圖像重建問題。

受上述研究的啟發(fā)，我們提出了一種基于導引圖像濾波和截斷全變分的CT重建算法。我們把SART算法重建的圖像當作輸入圖像，基于截斷TV最小化后的圖像當作導引圖像，采用導引圖像濾波來處理稀疏角CT重建問題。這樣既利用了基于截斷TV的圖像重建算法具有保護強邊緣的作用，又發(fā)揮了具有保護邊緣和平滑圖像作用的導引圖像濾波能傳遞圖像特征的優(yōu)勢。本文使用SART算法、SART-TV算法、SART-TTV算法、SART-TV-GIF算法和SART-TTV-GIF算法分別在30個投影、45個投影、和90個投影下重建高分辨率的頭部模型。SART-TTV-GIF算法可以有效地抑制噪聲，減少偽影并更好地還原圖像的邊緣細節(jié)。

1 算法

SART算法

CT圖像重建是一個典型的逆問題，實質上是求解線性方程組:

Af=p

(1)

式中,A=(aij)M×N為系統(tǒng)投影矩陣；aij為第j個像素對第i條射線對應投影值貢獻的權系數(shù)；f=(f1,f2,…,fN)Τ為待重建圖像；p=(p1,p2,…，pM)Τ,pi為第i條射線所對應的投影值;N為圖像像素的個數(shù)。

SART重建算法的基本思想如下:在迭代更新的過程中，先計算每個掃描旋轉角度下經過某一像素點的全部射線的誤差值，再對這些誤差值進行加權平均。該算法的迭代公式可以表示如下[4]:

(2)

1.2 截斷TV

Dou等提出了一種新的正則化方法，稱為截斷TV[7]。它可以有效的減少不重要的細節(jié)，保留顯著的邊緣，且不會造成圖像過于平滑。

截斷TV表示為[7]:

(3)

(4)

式中，u為平滑的結果圖；(ux,uy)為u的梯度。截斷TV只會懲罰幅度小于閾值ε的梯度，而對于大于閾值ε的梯度不會進行懲罰。

1.3 導引圖像濾波

導引圖像濾波已被廣泛用于噪聲減少、偽影去除、圖像增強和圖像重建。它是一種自適應權重過濾器，可以平滑圖像并同時保持邊緣。GIF中的重要假設是濾波器輸出的圖像fout能在一個正方形窗口i內被導引圖像Iguide線性表示。即導引圖像濾波的輸出圖像fout如下所示[8]：

(5)

式中，ak和bk在以像素點k為中心且半徑為r的正方形窗口wk中是常數(shù)。另外，ak和bk系數(shù)是通過最小化下面的能量函數(shù)來求解的：

(6)

1.4 SART-TTV-GIF算法

SART-TTV-GIF算法的求解步驟如下。

(1)SART步驟

SART-TTV-GIF算法的第一個步驟是使用SART算法進行投影與反投影運算來迭代更新圖像：

fn=SART(fn-1)

(7)

(2)TTV步驟

SART-TTV-GIF算法第二個步驟是求解截斷TV的最小值問題。目標函數(shù)如下所示:

(8)

式(8)的左邊第一項是數(shù)據(jù)保真項，第二項是截斷TV。其中截斷TV是非凸面的，并且由于其“截斷的”形狀而難以優(yōu)化。但是可以使用L0正則化來重新表達。以ε為參數(shù)，在(4)中定義的T(u)等價于

φ(u,l)=min{ε|l|0+|u-l|}

(9)

式中，|·|0為零次冪運算符，即如果l≠0則|l|0=1，否則|l|0=0。

式(9)又等價于

(10)

(11)

用拉格朗日乘數(shù)并引入布雷格曼距離，式(11)轉化為：

(12)

上述聯(lián)合最小化問題可以通過將其分解為幾個子問題來交替解決，如下所述：

①固定l1,l2,b1,b2,t1,t2來計算un

(13)

得出歐拉-拉格朗日方程:

fn-un+

(14)

式中，Δ為拉普拉斯算子，式(12)可用高斯-賽德爾迭代算法或者FFT算子快速求解。

由圖8可見，隨著聚乙烯靶高度的增加，中子探測效率增加，中子能量分辨率變化不大。因此，在聚乙烯靶高度方向(圖1中y方向)，應盡量選擇尺寸大的聚乙烯靶，這樣可以在保持能量分辨率的條件下，提高中子探測效率。

②固定un,l1,l2,t1,t2來計算b1,b2

應用收縮算子獲得該子問題的唯一極小值：

(15)

(16)

③固定un,l1,l2,b1,b2來更新t1,t2

(17)

(18)

④固定un,t1,t2,b1,b2來更新l1,l2

(19)

這個子問題的解決方案是

(20)

(21)

(3)GIF步驟

(22)

(23)

(24)

表1 SART-TTV-GIF算法的主要步驟

2 圖像結果的評價

本文是對重建的圖像進行客觀評價，就是通過選取峰值信噪比(PSNR)[11]和均方誤差(MSE)[11]來對重建圖像的降噪能力和精度等方面進行分析。其中，PSNR的值越大，MSE的值越小，重建圖像的質量越高。其計算公式如下[11]:

(25)

(26)

式中，ti,j,ri,j分別為原圖和重建圖像在像素處的像素值；MAX1為原圖的最大像素值。

3 實驗與分析

本文對高分辨率的頭部模型[12-13]進行了仿真模擬,實驗需要在如下環(huán)境下進行：Windows10操作系統(tǒng)，處理器為AMD RyZen 5 2500U with Radeon Vega Mobile Gfx 2.00 GHz，8.00 GB內存，Matlab R2018b，對大小為243×243的圖像進行重建。

3.1 無噪聲實驗

本文在180°范圍內等間隔采樣了90個投影、45個投影和30個投影，分別對SART算法、SART-TV算法、SART-TTV算法、SART-TV-GIF算法和SART-TTV-GIF算法來進行實驗。為了讓不同算法之間的比較更為公平，這五種算法的最優(yōu)參數(shù)和迭代步數(shù)都是通過評價指標和視覺效果的最優(yōu)來選取的。表2所示為這五種算法分別在30個投影、45個投影和90個投影下的參數(shù)值的選擇。

表2 參數(shù)值的選擇

續(xù)表

如圖1～圖4所示，橘色箭頭指出的重建圖像的區(qū)域有明顯的偽影。圖1為頭部腫瘤模型，原始模型中自帶有高斯噪聲[12]。在90個投影采樣下，五種算法重建得到的圖像都相對較好，但SART算法重建得到的圖像邊界模糊，SART-TV算法和SART-TV-GIF算法重建得到的圖像邊界較為清晰，SART-TTV算法和SART-TTV-GIF算法重建得到的圖像邊界最為清晰。在45個投影采樣下，SART算法重建得到的圖像有明顯的噪聲，重建得到的圖像質量不好；SART-TV算法和SART-TV-GIF算法重建得到的圖像有多處邊緣偽影；SART-TTV算法重建得到的圖像存在相對較少的偽影；SART-TTV-GIF算法重建得到的圖像幾乎沒有偽影，質量更好。在30個投影采樣下，SART算法重建得到的圖像有嚴重的噪聲；SART-TV和SART-TV-GIF算法重建得到的圖像有很多偽影；SART-TTV算法重建得到的圖像有較多的偽影；SART-TTV-GIF算法重建得到的圖像幾乎沒有偽影，重建結果更好。在多個稀疏角度下，通過SART-TTV-GIF算法重建的圖像質量更高，邊緣偽影更少，能更好地保護邊界結構。

圖1 參考圖像

圖2 90個投影(SART、SART-TV、SART-TTV、SART-TV-GIF、SART-TTV-GIF)

圖3 45個投影(SART、SART-TV、SART-TTV、SART-TV-GIF、SART-TTV-GIF)

圖4 30個投影(SART、SART-TV、SART-TTV、SART-TV-GIF、SART-TTV-GIF)

表3所示為對高分辨率的頭部模型重建結果的評價。

表3 圖像重建結果的評價

從表3可以清楚地看到，SART-TTV-GIF算法的MSE值最小，而PSNR值最大。所以，在這五種算法中，SART-TTV-GIF算法的重建效果最好。

3.2 加噪聲實驗

在實際的應用中，探測器采集到的投影數(shù)據(jù)通常是含有噪聲的。因此，我們向模擬的投影數(shù)據(jù)中添加了泊松噪聲，入射的X射線光子個數(shù)為1.0×103，添加泊松噪聲的均值為0。添加噪聲后各個算法的參數(shù)值的選擇如表2所示。

如圖5～圖7所示，隨著向模擬投影數(shù)據(jù)中添加泊松噪聲，不同算法重建圖像的質量都有所降低。在90個投影采樣下，五種算法重建得到的圖像都相對較好，但SART算法重建得到的圖像邊界模糊，SART-TV算法和SART-TV-GIF算法重建得到的圖像邊界較為清晰，SART-TTV算法和SART-TTV-GIF算法重建得到的圖像邊界最為清晰。在45個投影采樣下，SART算法重建得到的圖像有明顯的噪聲和偽影，重建得到的圖像質量不好；SART-TV算法和SART-TV-GIF算法重建得到的圖像有很多偽影，SART-TTV算法重建得到的圖像有較多的邊緣偽影；SART-TTV-GIF算法重建得到的圖像幾乎沒有偽影，質量更好。在30個投影采樣下，SART算法重建得到的圖像有嚴重的噪聲和偽影；SART-TV算法和SART-TV-GIF算法重建得到的圖像有較為嚴重的偽影，SART-TTV算法重建得到的圖像有相對較少的偽影；SART-TTV-GIF算法重建得到的圖像則存在少量偽影。在多個稀疏角度下，通過SART-TTV-GIF算法重建的圖像質量最高，偽影更少,能更好的保留邊界結構。

圖5 90個投影(SART、SART-TV、SART-TTV、SART-TV-GIF、SART-TTV-GIF)(含噪聲)

圖6 45個投影(SART、SART-TV、SART-TTV、SART-TV-GIF、SART-TTV-GIF)(含噪聲)

圖7 30個投影(SART、SART-TV、SART-TTV、SART-TV-GIF、SART-TTV-GIF)(含噪聲)

表4所示為對添加泊松噪聲的高分辨率的頭部模型重建結果的評價。表4與表3相比，MSE的值全都變大，PSNR的值全都變小，可以得出對圖像添加噪聲后，各個算法的重建圖像質量普遍偏低。

表4 含噪圖像重建結果的評價

從表4可以清楚地看到，SART-TTV-GIF算法的MSE值最小，而PSNR值最大。所以，在這五種算法中，SART-TTV-GIF算法的重建效果最好。

4 結論

本文利用五種圖像重建算法對高分辨率的頭部模型進行了仿真實驗，對稀疏角投影數(shù)據(jù)進行了圖像重建，并對各種算法的重建結果進行了比較。實驗結果表明，SART-TTV-GIF算法得到的重建圖像質量更高，可以有效抑制噪聲，減少偽影，并能更好地保留邊界結構，并且運行速度很快。但是，SART-TTV-GIF算法包含許多參數(shù)，調試起來很耗時。本文的參數(shù)都是通過反復的實驗來選取的，而且本文僅對高分辨率的頭部模型進行了測試。在未來的工作中，將進一步討論參數(shù)的最優(yōu)選取方式并測試各種模型以測試算法的有效性，同時也會研究將該算法應用到有限角CT重建領域當中。