基于三維加權和雙域濾波的螺旋CT偽影校正

牛曉偉，孔慧華，邸云霞

(1.中北大學 數學學院，太原 030051; 2.信息探測與處理山西省重點實驗室，太原 030051)

0 引言

螺旋計算機斷層掃描技術(CT， computed tomography)自問世以來，具有了更快的掃描方式和更高的軸向分辨率，而越來越多的高分辨率錐束探測器需要使用專用的三維錐束重建算法，以提供更好的圖像質量，螺旋CT已經廣泛應用在了工業無損檢測、臨床診斷等領域[1-3]。螺旋CT的掃描軌跡是邊旋轉邊平移，即螺旋掃描軌跡，掃描模式中的視角實際上是分布在一條三維曲線上，給定一定數量的視角，視角的取樣密度取決于沿z軸方向的幾何覆蓋，z方向的覆蓋范圍越長，視角的采樣密度就越低。

對于目前的單源多排探測器CT(MDCT, multi-row detector CT)掃描儀，螺距通常被限制在1.5以下，這種限制是由于如果使用目前最先進的重建算法來重建圖像，螺距大于1.5時就會出現圖像偽影。因此，在高螺距螺旋CT重建問題中存在兩個內在的科學挑戰：一是角度樣本區間大，如圓形掃描中的視角欠采樣重建問題；二是視角范圍被截斷，如有限視角重建問題。也就產生了兩類圖像偽影：有限視角偽影和沿z軸方向的數據不一致造成的偽影，這些偽影基本上都是由掃描的三維性質和近似算法的使用引起的[4-5]。

錐束偽影是螺旋CT中常見的一種偽影，產生的一個根本原因是錐角和螺距過大而引起的共軛射線不一致，并且偽影會隨著重建平面與中心平面之間的距離增大而變嚴重[6-7]。螺旋FDK算法需進行投影數據重排為平行數據，這會降低空間分辨率，造成損失，僅在錐角很小時比較有用，無法在錐角或螺旋間距增加的情況下提供可接受的重建精度。Balogh等[8]提出了一種基于統計的迭代環去除算法，有效地檢測和校正投影數據中的錯誤像素值，而不會在重建圖像中造成新的偽影，可以去除缺陷檢測單元產生的環偽影。Han等[9]提出了一種在圓軌道錐束CT系統中減少錐束偽影的雙能方法，可以有效恢復扭曲的骨性結構和強度值，尤其在大錐角區域。Tang等[7]基于圓柱形探測器，提出了一種三維加權(3D, three-dimensional)螺旋FDK重建算法，可以在中等錐角和大螺距下顯著提高重建的精度。Kong等[10]基于平板探測器，將錐形束數據轉換為錐形平行幾何體，提出了一種3D加權螺旋PFDK算法，在錐角和螺距較大時可以提供良好的重建結果。Magkos等[11]調整CT軌跡的直接迭代重建算法來估計錐束偽影并對其進行補償，產生了質量更高的重建體積。

而風車偽影是由z軸方向采樣不足引起的，常出現在z軸解剖結構和密度變化很大的部位，在其重建圖像上圍繞中心點呈黑/白圖案組成風車狀的偽影，尤其是在其他參數保持一定的情況下，重建圖像中風車偽影的嚴重程度與螺距的大小成正比，準直器層越厚風車偽影就越嚴重[5,11]。針對螺旋CT存在的風車偽影問題，可以通過提高z軸方向的采樣率得到改善，人們提出了各種方法來減少風車偽影或提高縱向分辨率，這些方法可分為三類：正弦圖上采樣、超分辨率重建和圖像后處理。李嫣嫣等[12]提出一種結合黃金角變密度螺旋采樣、并行成像和基于同倫l0范數最小化的壓縮感知的圖像重建的三維動態磁共振成像方法，可以提高欠采樣率，實現較高的空間分辨率和時間分辨率。減少風車偽影最常用的方法之一是改進CT系統硬件中的縱向采樣，如德國西門子公司的z軸動態飛焦技術(z-FFS, z-Flying Focal Spot)利用焦斑在縱向上的周期性運動將同時獲取的切片數量增加一倍，以提高縱向分辨率和減少風車偽影[13]。Brown等[14]提出了一種基于圖像域總方差最小化的風車偽影抑制方法，該方法能夠去除風車偽影，同時保持圖像中解剖結構的分辨率。Shi等[15]采用改進的低秩和全變分正則化從低分辨率圖像中初步恢復細節和邊緣，然后將高分辨率圖像與受風車偽影影響的原始重建圖像相結合，可以獲得更準確的圖像。黃銅港等[16]提出一種基于雙域濾波與距離變換的算法對風車偽影圖像進行校正，能夠在校正風車偽影的同時較好地保留原始螺旋CT圖像中的細節信息。

上述方法在一定程度上解決了螺旋CT重建中錐束偽影和風車偽影的影響，但是都有一定的局限性，例如，部分投影數據由于消除數據冗余而被浪費，噪聲特性或劑量效率降低；空間分辨率下降，因為必須沿z方向至少執行一次插值；由于拉東平面或Tam窗的參與導致不可接受的時間分辨率。但高螺距螺旋掃描在臨床實踐中是非常必要的，比如，選擇高螺旋距掃描可以減少病人輻射劑量，可以提高重建圖像質量等等。因此本文提出使用三維加權螺旋FDK算法[6]和改進的雙域濾波算法[17-19]對螺旋CT圖像進行重建和校正，雙域濾波可以在保留更多細節的同時保留邊緣等高對比度特征，因此能有效解決上述問題并在較大的錐角和螺距下保持更好的圖像質量。

1 理論與方法

1.1 基于3D加權的錐束偽影去除

螺旋CT中，隨著錐角和螺距的增大，重建圖像質量變差，這是由于重建時共軛射線不一致引起的。X射線在機架內旋轉180°后，可以獲得相同方向的X射線，此射線就被稱為共軛射線，且兩條射線的投影會存在差異，這種差異我們通常稱為共軛射線不一致。如圖1所示，給定一個要重建的點P(x,y,z)和一條穿過它的射線SP(即直射射線)，那么就存在一條也穿過點P并與之對應的共軛射線S′P，其中α和αc分別對應直射射線和共軛射線的錐角。由于直射射線的錐角不等于共軛射線的錐角，因此如果重建的平面離中心平面越遠，重建圖像的質量就會越差。Patch(2004)提出了一種3D權重，用于抑制原生錐束幾何體中錐束重建的偽影，其對一對共軛射線中錐角較小的射線給予有利權重，并對錐角較大的射線給予不利權重，便可以提高螺旋重建精度。基于這一情況，下面給出了一種基于平板探測器的3D加權螺旋FDK重建算法[10]。

圖1 螺旋CT投影數據中一對共軛射線的示意圖

首先給出二維加權(2D, two-dimensional)FDK算法，在2D加權的基礎上考慮錐角的影響，再給出一種3D加權算法，從而得到一種3D加權FDK算法：

(1)

(2)

(3)

然后將預加權和過濾后的投影進行反投影來重建圖像，在圖1所示的錐形平行幾何中，3D加權螺旋FDK重建算法的解析表達式為:

(4)

(5)

其中:式(5)是3D加權函數，k是一個隨螺距變化的參數。并且加權函數w3d(θ,l,α)滿足:

w3D(θ,l,α)+w3D(θc,-l,αc)=1

(6)

2D加權函數可以抑制由射線冗余引起的偽影，而3D加權函數可以抑制由錐角和螺距過大引起的偽影。

1.2 基于雙域濾波的風車偽影去除

雙域濾波包括雙邊濾波和短時傅里葉變換(STFT, short-time fourier transform)兩個部分，其主要思想是將圖像分為兩個層次，包括基礎層和細節層，其中基礎層具有高對比度的特點，可以采用空間域的算法進行處理，細節層具有低對比度的特點，能采用變換域的算法進行處理。它是在空間域使用雙邊濾波，可以保留邊緣等高對比度特征，在變換域使用帶小波收縮的STFT，可以很好的保留細節，將空間域和變換域結合到一起比單獨使用其中一種效果更好[17]。而螺旋CT的螺距值越大，螺旋掃描一周所涉及的用于重建的探測器排數就越多，重建圖像中的風車偽影也會變嚴重。首先輸入一幅含有風車偽影的圖像f，同時將中間去噪圖像作為導引圖像來繼續后面的迭代，用于調整STFT的系數，并且隨著迭代的增加使得噪聲越來越小，不斷逼近原始圖像。雙域去噪采用三步迭代進行，三次迭代以后，引導圖像逐漸逼近原始干凈圖像。去噪算法的模型可以定義為:

f=g+η

(7)

式(7)中f表示含風車偽影的圖像，g表示原始圖像，η表示方差為σ2的加性噪聲。雙域濾波算法的目標是從噪聲污染圖像f來估計原始圖像g。以下是對雙域濾波算法的雙邊濾波、STFT和小波收縮的詳細介紹。

1.2.1 雙邊濾波

● 誤食外形較為光整的固體物品：如紐扣、果核、小硬幣等。若是哽在喉部，可嘗試催吐，若是吞入腹中，可等待與糞便中排出，若是大量吞入腹中引起消化道不適，很緊急就醫。

雙邊濾波是一種非線性濾波器，僅是在局部鄰域樣本的加權平均值，其中權重是基于中心樣本和相鄰樣本之間的空間和輻射距離計算的，通過這種方式，邊緣得到很好的保留，而噪聲被平均化。它需要對導引圖像和含風車偽影的圖像同時進行雙邊濾波，定義了一個正方形窗口Np，Np以像素p為中心，r為半徑，導引圖像和含風車偽影圖像像素點p雙邊濾波后的高對比度值為：

(8)

(9)

(10)

雙邊濾波的邊緣保護特性是因為使用了范圍核(以及空間核)，用于控制邊緣附近的擴散，并且空間核和距離核都是高斯函數。其中強烈的不連續性通常發生在邊緣附近，這由范圍核處理，然后用于抑制空間的擴散。另一方面，范圍核在強度平滑變化的區域中變得不起作用，然后空間核進行處理。雙邊濾波器作為標準的擴散濾波器，空間核和范圍核一起在均勻區域中執行平滑，并同時保留邊緣。

1.2.2 STFT

通過式(8)、(9)得到了圖像的基礎層，然后將圖像信息從空間域轉換到頻域，即將原圖像分別減去導引圖像和含噪圖像的高對比度去噪圖像來得到各自的低對比度圖像，然后進行STFT，為小波收縮做準備。以時間信號為例，短時傅里葉變換就是對信號先乘以僅在一段時間(圖像信號是空間)上不為零的窗函數，然后對加窗的函數進行傅里葉變換。STFT所得系數Gp,v和Sp,v的自變量為頻率v，定義域為頻率窗口Vp，該窗大小與Np相同。分別得到它們的傅里葉系數Gp,v，Sp,v及方差為:

(11)

(12)

(13)

STFT步驟加強了中心部分的細節特征，它每次迭代噪聲都會減少，而只會引入很少的偏差，在進行離散傅里葉變換時需要執行加窗操作，是為了避免圖像邊緣出現新的偽影。雙邊濾波核采用了高斯函數，STFT過程就等同于對低對比度圖像進行了Gabor變換。

1.2.3 小波收縮

(14)

其中:γf表示小波收縮因子調整參數，通過對Gp,v構造的小波收縮系數Kp,v，對Sp,v收縮得到低對比度值，將細節信號的信息反變換回空間域，即:

(15)

其中:Mp為Sp,v對應的頻率分布矩陣，|Mp|為Mp中的元素個數。那么，對于雙域濾波從含噪圖像f還原未知圖像g，便可以合理估計為：

(16)

雙邊濾波核掩蓋了可能導致小波收縮振鈴的高對比度邊緣，而局部傅里葉變換檢測并保留了雙邊濾波核可能會模糊的低對比度重復結構，實現起來較為容易。

1.2.4 改進的雙邊濾波核函數

雙域濾波算法經過多次迭代會造成時間過長，并且將原始含風車偽影圖像作為引導圖像來進行下一步迭代，當風車偽影過于嚴重時，需要增強雙域濾波的強度，進而會造成圖像細節的消失和邊緣部分變得模糊。

式(10)中的d(p,q)指的是像素間的歐氏距離，用來度量鄰域間的相似性。當偽影過于嚴重時，對于結構和細節信息比較豐富和復雜的區域，由于歐氏距離較大就會導致相似度低的鄰域被賦予較大的權重。此時，歐氏距離就不能較好的反映鄰域間的相似性，達不到抑制風車偽影同時保留細節信息的目的。本文將歐氏距離改為漢明距離來度量鄰域間的相似性，公式如下：

(17)

其中:L(p,q)指的是像素p、q鄰域之間的漢明距離。兩個向量x、y∈{0,1}n之間的漢明距離定義為:

(18)

其中:x=(x1,x2,…,xn)T，y=(y1,y2,…,yn)T，⊕為模2加運算[20]。

改進后的雙域濾波可以選取更小的σ值進而能有效避免原始雙域濾波的局限性。

2 實驗結果

2.1 實驗設置

本研究的主要目的是去除錐束偽影和風車偽影，將通過仿真實驗驗證提出算法的有效性。文中所有算法都在MATLAB和C++的混合模式下實現，接口在MATLAB中，大規模計算全部都在C++中。實驗中，CT系統的掃描方式為螺旋掃描，射線源到旋轉中心的距離為200 mm，探測器為平板探測器，探測器到旋轉中心的距離為50 mm，探測器單元尺寸為1 mm×4 mm，探測器陣列大小為350×20，掃描的圈數為3圈，每圈均勻采樣，并且每圈采集的投影數為360，重建圖像尺寸為256×256×256。實驗分別對相同材料和不同材料的植入物進行了研究，并同時采用原始雙域濾波算法對偽影進行校正，用于與本文所提算法進行對比。本文所有實驗中雙域濾波均進行3次迭代。圖2是對不同材料植入物編號的示意圖。

圖2 不同材料植入物編號

圖2的模型植入物參數見表1。

表1中，A為植入物的灰度值，O為植入物的中心坐標，a，b，c分別為植入物在X，Y，Z軸方向上的半軸長系數。

表1 模型植入物參數

2.2 實驗分析

圖3、4、5為同種材料的植入物在螺距為1時第107層、105層和104層的重建圖像與校正圖像，σ為雙域濾波的參數值。第一列是3D加權螺旋FDK算法重建圖像，由于重建圖像的錐角依次增大，掃描對象結構變窄，圖像沿旋轉軸的欠采樣而遭受到更嚴重的風車偽影；第二列是螺旋FDK算法的重建圖像，由于其忽略了直射射線和共軛射線之間的錐角差異，因此圖像中還出現了較為明顯的錐束偽影，并且隨著錐角的增大，圖像中的錐束偽影也隨之增加，嚴重影響了CT圖像的質量；第三列是經過3D加權和原始雙域濾波的校正圖像，由于圖5(a)的風車偽影過于嚴重，圖5(c)的去噪圖像需要增加雙域濾波的強度，雖然可以去除一部分風車偽影，卻造成了邊緣和細節信息的丟失，并且偽影仍有殘留；第四列是本文算法處理后的重建圖像，圖5(d)不僅可以去除大部分的錐束偽影和風車偽影，并且保留了更多的細節和邊緣信息，效果要比圖5(c)好。可以看到，圖5(c)和圖5(d)的σ值分別為0.019 8和0.000 15，圖5(d)的雙域濾波強度遠遠小于圖5(c)，因此可以避免圖像變得模糊。

圖3 同種材料107層切片螺距為1的重建圖像與校正圖像

圖4 同種材料105層切片螺距為1的重建圖像與校正圖像

圖5 同種材料104層切片螺距為1的重建圖像與校正圖像

圖6、7、8為同種材料植入物105層切片在螺距為0.75、1和1.25下的重建圖像與校正圖像。第一列是采用3D加權螺旋FDK算法的重建圖像，圖像沿著高對比邊緣附近出現了螺旋狀圖案，并且螺距越大，重建出來的圖像質量越差，風車葉片的數量也越來越多；第二列是螺旋FDK算法得到的重建圖像，圖像中存在風車偽影和錐束偽影，并且隨著螺距的增大，植入物之間的黑白帶狀偽影也變得更明顯；第三列是原始雙域濾波對第一列重建圖像去噪后的圖像，錐束偽影和風車偽影基本去除，但結構的邊緣變得模糊，較小的植入物也消失不見；第四列是采用本文算法處理后的校正圖像，去除了圖像中的錐束偽影和風車偽影，保留了更多的細節特征，并且結構邊緣也更加清晰。

圖6 同種材料螺距為0.75的重建圖像與校正圖像

圖7 同種材料螺距為1的重建圖像與校正圖像

圖8 同種材料螺距為1.25的重建圖像與校正圖像

圖9為不同材料植入物不同螺距下第105層的螺旋CT重建圖像與校正圖像，植入物的參數如圖2和表1所示，其第一、二和三行分別是螺距為0.75、1和1.25的重建圖像。由圖9第一列可以看出，由于植入物的灰度值不同，風車偽影的嚴重程度也不相同，并且植入物與背景材料的灰度值差距越大，風車偽影也越嚴重。本文算法采用的σ值遠遠小于原始雙域濾波的σ值，因此可以在去除了絕大部分的錐束偽影和風車偽影的同時一定程度上保留圖像的邊緣信息，由此可看出本文算法在整體圖像質量與細節特征的恢復上有著明顯的優勢。

圖9 不同材料不同螺距下的重建圖像與校正圖像

由上述實驗結果可以看出，重建平面如果離中心平面越遠，即錐角越來越大時，重建出來的圖像質量就越差，錐束偽影也就越嚴重，而3D加權螺旋FDK算法能夠有效減少錐束偽影，通過調節參數k的值來調整3D加權函數的權重，并且隨著錐角的增大，參數k的值也隨之增大。對于風車偽影，由圖3、4、5可以看出，隨著錐角的增大，風車偽影越來越嚴重，重建圖像的質量變差；由圖6、7、8可以看出，隨著螺距的增大，風車偽影葉片的數量越來越多。原始雙域濾波算法可以有效去除風車偽影，但當風車偽影過于嚴重，雙域濾波的強度過大時，雙邊濾波中核函數的歐氏距離較大便可能會導致相似度低的鄰域被賦予較大的權值，最終造成邊緣模糊和細節丟失等問題，改進后的雙濾波算法用漢明距離來度量鄰域間的相似性，可以選取更小的σ值，有效避免原始雙域濾波的局限性，可以在抑制風車偽影同時保留細節信息。

3 結束語

本文提出了一種基于三維加權和改進的雙域濾波的螺旋CT偽影校正算法，提高了螺旋CT重建圖像中偽影校正的能力。首先使用3D加權螺旋FDK算法去除錐束偽影，再使用改進的雙域濾波算法對圖像中的風車偽影進行校正。3D加權螺旋FDK算法可以抑制錐角和螺距增大對重建結果的影響，以此來去除錐束偽影，而雙域濾波將空域和變換域結合起來，與單獨使用雙邊濾波或者STFT小波收縮相比效果要更好，可以在很大程度上去除風車偽影。實驗結果表明，本文算法有效去除了螺旋CT重建圖像中大部分的偽影，并且可以保留更多的細節特征和邊緣信息。