改進FDK算法的錐形束CT圖像重建

吳澤靜，徐田猛，賈思怡，代穎涵，陸佳華，王遠軍

(上海理工大學 醫療器械與食品學院，上海 200093)

1 前 言

X射線計算機斷層成像(X-ray Computer Tomography，X-CT)技術是一種廣泛應用于工業檢測、醫學檢查、安保檢測等領域的技術，但它只能獲得二維斷層圖像，展現某一截面的生理信息[1].錐束CT(Cone Beam Computer Tomography)因具有精度高、輻射小、成像速度快、重建圖像分辨率高等特點，所以成為目前CT成像技術領域研究熱點方向之一[2].1984年，Feldkamp、Davis、Kress提出了基于錐形光束、圓形軌跡掃描和完全數據條件下的快速重建算法即FDK算法[3]，它是目前認為最有效的三維濾波反投影重建算法，廣泛應用于CT系統中.1993 年，Ge Wang等人將FDK的圓形軌道拓展為任意軌道，提出了 G-FDK算法[4].1999 年，Turbell 將錐束 CT 投影數據重新排列為三維平行束的投影數據，提出了 P-FDK算法[5].2000 年，Grass 等人基于 P-FDK 對截斷后的投影數據進行數據重建，相繼提出了 T-FDK[6]和 HT-FDK[7]算法.三維錐束CT中的FDK 算法因為具有計算速度快，圖像重建質量好，數值運算穩定，便于硬件加速實施等優點成為當今人們最受推崇的圖像重建算法[8].

FDK算法是二維扇束濾波反投影的三維推廣，因此對于采集到投影數據去噪，濾波器顯得至關重要.1971年，后來印度學者Lakshminarayanan和Ramachandran提出了離散形式的R-L濾波函數，其簡單實用，重建圖像輪廓清晰[9].1974年，Shepp和Logan提出了S-L濾波函數[10]，該函數對高頻段重建質量較好，而低頻段質量不高.近些年來，研究者們又提出和改進了一些新的濾波函數，如Hanning濾波函數、Hamming濾波函數、R-L和S-L混合濾波函數等[11].

2005年，Wei Yuchuan 等人基于理想斜變濾波函數的設計理論，引入了δ函數，得出一種形式簡單的濾波函數，記為 NEW 濾波函數[12].經過研究發現此函數旁瓣衰減比傳統的S-L濾波函數要更快，會使重建出來的圖像擁有更高的密度分辨率.基于前人獲得的NEW& R-L混合濾波函數，本文通過研究加權系數對NEW& R-L混合濾波函數的影響，從而找出最佳重建效果的加權系數.另外，本文提出了一種針對富含噪聲的圖像的濾波方法NEW& S-L，通過加權優化，使該濾波方法能更加有效地濾除圖像中的噪聲.

2 FDK算法及其改進

2.1 FDK算法實現步驟

FDK算法是一種針對錐束幾何圓形掃描軌跡提出的近似的重建算法，錐角越小重建的誤差也越小.其方法和扇形掃描的濾波反投影算法類似，具體步驟如下所示：

1)使用加權函數對投影數據進行加權：

pw(β，a，b)=p(β，a，b)×w(a，b)

(1)

2)對加權后的投影數據進行水平方向上的一維濾波：

pf(β，a，b)=pw(β，a，b)*g(a)

(2)

其中g(a)為Ramp濾波函數

3)對不同投影角度下濾波后的投影數據沿X射線反方向進行反投影：

(3)

其中，U(x，y，β)=R+xcosβ+ysinβ.

2.2 常用濾波函數

由于在直接反投影(Back Project，BP)過程中，容易產生星狀偽影，導致圖像的失真[13]，所以在反投影重建之前必須要對投影數據進行一維濾波，然后再對修正后的數據進行無偽影重建，這個過程就是進行濾波反投影的過程.接下來介紹幾種常見的傳統濾波函數：

R-L濾波函數是由印度學者Rammachandran和Lakshm Inarayan提出的，它是利用窗函數截斷理想濾波器，得到近似理想濾波器：

(4)

R-L濾波函數采用序列為：

(5)

R-L濾波器形式簡單，重建出的圖像邊緣輪廓清晰，缺點就是由于截斷了濾波函數就會導致震蕩響應，出現 Gibbs 現象.尤其是當被測物體含有噪聲或者衰減系數變化較大時Gibbs效應會更加明顯，使得重建圖像質量下降[14].

(6)

S-L濾波函數采樣序列：

(7)

S-L函數由于旁瓣衰減的較快，有效地減少了震蕩，所以可以減少Gibbs效應.但是由于S-L函數的主瓣較低，所以降低了重建圖像的空間分辨率，在低頻段部分，S-L的重建效果不如R-L重建的效果.

為了使濾波函數兼具R-L濾波函數主瓣高以及S-L濾波函數旁瓣衰減速度快的特點，一些學者提出將R-L和S-L濾波函數結合在一起.通常，主瓣高且窄，圖像的空間分辨率就高；旁瓣衰減迅速，可提高其密度分辨率.混合濾波器就對此能進行很好的折中，對圖像的噪聲有著很好的抑制，同時也能很大地改善圖像的分辨率.根據實際所需要的頻率響應函數，通過調整R-L濾波器和S-L濾波器的加權系數，達到一定程度上抑制噪聲并提高空間分辨率的目的.其采樣序列表達式：

(8)

其中k1和k2滿足以下關系：k1≥0，k2≥0，k1+k2=1.由上式可知當k1為0時是S-L濾波器，k2為0時是R-L濾波器.經過對1序列峰值及旁瓣衰減快慢進行比較，當k1=k2=0.5時，R-L& S-L混合濾波函數的峰值保持在較高的水平，同時，旁瓣衰減速度較R-L濾波函數也在一定程度變快不少.圖1給出了R-L，S-L，R-L& S-L濾波函數的采樣序列：

圖1 R-L，S-L，以及R-L& S-L濾波函數采樣序列

2.3 本文改進的濾波函數

NEW濾波函數是Wei Yuchuan等人基于δ函數提出的一種重建效果優于R-L濾波器的函數.NEW濾波函數離散序列采樣形式為：

(9)

這種函數比R-L濾波器計算要更加簡便，且此函數不僅能保證空間分辨率無顯著降低，并且有效地減小了Gibbs現象，因為該濾波函數的旁瓣比 S-L 函數收斂的更快，對投影數據中的高頻噪聲的抑制作用會更好.結合傳統的將R-L濾波函數與S-L濾波函數結合在一起的思路，考慮將該NEW濾波函數與R-L濾波函數結合在一起，這樣既能保證所形成的混合濾波函數有較高的主瓣，同時旁瓣衰減速度比傳統的R-L& S-L濾波函數要快，因此此種混合濾波器的重建效果要比傳統R-L& S-L的混合濾波器的重建效果要好.其采樣序列為：

(10)

其中k1和k2為加權系數，滿足k1≥0，k2≥0，且k1+k2=1.顯然，當k1為0時是NEW濾波器，k2為0時是R-L濾波器，下文中所采用的k1=0.7，k2=0.3.不同k1和k2值組成的混合濾波器空間域的主瓣及副瓣的幅值和寬度不同，使用主瓣高且窄的濾波器所重建出圖像的空間分辨率更高，而使用副瓣衰減迅速的濾波器重建出的圖像的密度分辨率更高.因此，可根據混合濾波器的主瓣及副瓣的幅值和寬度來選擇不同的k1和k2值，以滿足不同需求.

R-L& S-L濾波函數和NEW& R-L濾波函數，雖然折合了R-L濾波函數主瓣較高的特性，但旁瓣衰減速度仍然很慢，并且由于旁瓣幅度的原因，不能很好地抑制投影圖像中的噪聲.因此，針對含有較多噪聲的圖像，本文提出了一種混合濾波數，將S-L濾波函數與NEW濾波函數結合在一起，構成NEW& S-L混合濾波器，并通過改變比較加權系數，來抑制高頻振蕩的吉布斯效應，從而對含有噪聲的圖像達到很好的重建效果.以下便是NEW& S-L濾波函數的采樣序列：

(11)

其中λ1≥0，λ2≥0，λ1+λ2=1，當λ1=0，λ2=1時，該函數為S-L濾波函數，當λ1=1，λ2=0 時，該函數為NEW濾波函數.取不同的λ1，λ2加權系數組合，可得到不同的主瓣高度以及不同的旁瓣衰減速度，為使重建圖像空間分辨率和密度分辨率都高，選取的加權系數組合應使旁瓣衰減地較快的同時，主瓣高度也應盡可能的高.圖2第1列是R-L，NEW，R-L& NEW的采樣序列，第2列是是加權系數λ1=0.3時的NEW& S-L濾波函數的采樣序列.

圖2 NEW& S-L濾波函數采樣序列

2.4 FDK算法中的插值方法

FDK算法中最關鍵的一步就是插值算法的選擇，這是因為投影數據是離散采樣得到的，在反投影過程中，會出現像素的投影地址“對不準”的問題，因此需要采用插值法來估算投影地址的投影值.由于使用最近鄰插值將會出現塊狀和鋸齒，因此本文將使用插值性能更好的雙線性插值和雙三次插值對圖像進行重建[15].

1)雙線性插值

雙線性插值是指分別在x，y方向上進行一次線性插值.已知函數f在點Q11=(i，j)，Q12=(i，j+1)，Q21=(i+1，j)和Q22=(i+1，j+1)4個的值.分別在x，y方向上進行線性插值，可得到P(i+u，j+v)處的灰度值為：

f(i+u，j+v)≈(1-u)(1-v)f(Q11)+u(1-v)f(Q21)+

(1-u)vf(Q12)+uvf(Q22)

(12)

2)雙三次插值

雙三次插值是一種比雙線性插值更為復雜的一種插值算法，其不僅考慮反投影插值點周圍4個相鄰投影點灰度值的影響，還考慮了點p(i+u，j+v)周圍16個投影點灰度值的影響.雙三次插值需要選取基函數來計算16個像素點的權重，最終p(i+u，j+v)處的灰度值即為周圍16個點的灰度值的加權疊加[16].基函數S(w)的表達形式如下：

(13)

雙三次插值方法公式為：

f(x，y)=ABC

(14)

其中A，B，C為3個矩陣，其形式分別為：

A=[S(1+u)，S(u)，S(1-u)，S(2-u)]

(15)

(16)

C=[S(1+v)，S(v)，S(1-v)，S(2-v)]T

(17)

3 實驗結果與分析

3.1 傳統濾波函數和NEW濾波函數重建效果

本文基于三維頭模型Shepp-Logan來進行錐束CT的濾波反投影重建，仿真實驗中探測器采樣點數 N 設為 256，測試所使用的電腦配置Intel(R)Core(TM)i5-3210M CPU，主頻為：2.5GHz，運行內存為8GB，測試平臺為MATLAB 2016a.設置重建圖像的大小為 256×256.其中設置探測器通道個數 N=256，采樣點數為 256，源到旋轉中心的距離為 44，源到探測器中心的距離為 88，虛擬探測器長度為 2，采用雙線性插值法，得到三維 Shepp-Logan 頭部模型各中心斷面圖像及0°，45°，90°下的投影數據如圖3所示.

圖3 頭模型中心斷面圖像和投影

為了探究NEW濾波器的重建效果，本文將NEW濾波器重建的結果與R-L濾波器，S-L濾波器重建結果進行了比較，結果如圖4所示.

圖4 3種濾波函數重建結果(第1行為R-L濾波函數重建圖像在x=0，y=0，z=0處的斷面，第2行為S-L濾波重建的結果，第3行為NEW濾波函數的重建結果)

由圖4可知用R-L濾波器的震蕩效應非常明顯，而S-L濾波器的震蕩效應沒有R-L那么明顯，重建出來的圖像邊緣更加平滑些，并且濾除噪聲的效果比R-L要好得多.NEW濾波器的重建效果是3種濾波器中最好的，不僅有效地減弱了Gibbs效應，得到的圖像空間分辨率和密度分辨率相對于前面兩種濾波器都有所提高.

3.2 加權系數k1對NEW&R-L混合濾波函數的影響

為了定量評價重建精度，本文采用歸一化均方距離判據d，和圖像清晰度c來評價仿真重建圖像質量，二者的計算公式分別為：

(18)

(19)

本節主要是研究加權系數k1對NEW&R-L函數的影響，改變加權系數k1，得到加權系數k1與歸一化均方距離d以及圖像清晰度c之間的關系，如圖5所示.

通過對圖5的分析，可以看出，在圖像未添加任何噪聲的情況下，加權系數為0.7時，歸一化均方距離d值最小，說明重建出來的圖像與原圖的差距最小，最接近原圖；另外，由加權系數與圖像清晰度之間的關系可以看出，當k1逐漸增大時，圖像清晰度c也逐漸增大，當k1=0.7左右時，圖像清晰度隨k1的變化越來越緩慢，說明圖像清晰度已經增大到一個較高的水平，再增加k1清晰度也不會有一個較大的變化，因此，為了得到一個較好的重建結果，本文選擇當k1=0.7，k2=0.3時，利用NEW& R-L混合濾波函數對圖像進行重建.為了突出其重建效果，將NEW濾波函數重建結果，R-L& S-L濾波重建結果與其進行比較，實驗結果如表1所示.

圖5 加權系數k1與歸一化均方距離以及圖像清晰度的關系

表1 3種重建結果比較

由表1可知選用加權系數k1=0.7的NEW& R-L混合濾波函數進行重建時，歸一化均方距離d值不僅較小，且重建圖像清晰度也較其他兩種濾波函數高.說明NEW& R-L混合濾波函數能在一定程度上抑制重建時的震蕩效應，并且抑制噪聲的作用更加明顯，空間分辨率和密度分辨率都維持在很高的水平.

3.3 NEW&S-L濾波函數去除噪聲驗證

為得到使濾波效果最佳的加權系數，改變加權系數λ1，得到加權系數λ1與歸一化均方距離d，以及圖像清晰度c之間的關系如圖6所示.

圖6 加權系數λ1與歸一化均方距離和圖像清晰度的關系

由圖6可知，隨著加權系數λ1的增大，歸一化均方距離d先緩慢增加，后增加速度加快，因此根據歸一化均方距離可以選擇加權系數λ1小于0.5的數值，再結合加權系數λ1與圖像清晰度的關系，可知，加權系數λ1小于0.5的幾個數值重建出來的圖像質量差不多，因此，本文選擇以加權系數λ1=0.3為例，對加入均值為0，方差為1，噪聲強度為5%的高斯噪聲的三維頭模型Shepp-Logan進行濾波重建，為突出該新型混合濾波器對噪聲的濾除性能，用R-L& S-L濾波器以NEW& R-L濾波器進行對照，重建結果如表2所示.

表2 R-L& S-L，NEW& R-L和NEW& S-L濾波重建結果

由表2所知，NEW& S-L混合濾波函數濾除噪聲的效果要比R-L& S-L濾波函數以及NEW& R-L濾波函數要好得多.此外，重建出的圖像邊緣平滑清晰，震蕩效應較低，小的細節也能清晰的重建出來，而R-L& S-L濾波函數在小細節的方面重建的不好，甚至不能重建出來.

3.4 不同插值算法對重建圖像質量的影響

為了探究不同插值法對重建圖像的影響，本文將在使用NEW&R-L混合濾波函數的基礎上，分別將雙線性插值方法和雙三次插值方法應用于FDK算法中，從而得到重建圖像.結果如表3所示.

表3 雙線性插值和雙三次插值重建后比較結果

由表3可知，使用雙三次插值可以觀察到更多的細節，并且由于雙三次插值是考慮到周圍16個點的灰度值的影響，所以他能夠減輕雙線性插值導致的邊緣失真現象，重建出來的圖像更加接近真實的圖像.

4 結 語

本文通過對FDK算法中的濾波器的研究，提出了一種針對含有大量噪聲圖像的濾波方法NEW& S-L，并對其加權系數λ1進行了研究，得到濾除噪聲最佳的加權系數.并且在原有NEW& R-L濾波函數的基礎上，對其進行加權優化.實驗結果表明，NEW& S-L濾波函數的濾除噪聲的本領要比其他的濾波函數濾波效果更好.另外，本文比較了雙線性插值和雙三次插值的重建效果，得出了雙三次插值能夠顯示更多的細節并且它能夠減輕雙線性插值所導致的邊緣失真的現象，使得重建圖像更加的接近原圖像.但是使用雙三次插值的計算量較大，重建時間較長.因此接下來的工作主要是如何改進插值算法以及硬件設備，從而使重建時間更短，更能符合實際應用的需求.