基于Wiener濾波的反投影圖像重建算法

馬 敏 趙 亮 姬晶晶

（中國民航大學(xué)航空自動化學(xué)院 天津 300300）

計算機層析成像(CT)是通過對物體進行不同角度的射線投影測量而獲取物體橫截面信息的成像技術(shù)。CT技術(shù)的核心是投影重建圖像，其實質(zhì)是由掃描所得的投影數(shù)據(jù)反求出成像平面上每個點的衰減系數(shù)值。圖像重建方法大致分為[1]：變換法和級數(shù)展開法。濾波反投影算法是目前應(yīng)用最廣泛的基于變換法的圖像重建算法[2,3]，具有重建速度快、空間和密度分辨率高等優(yōu)點，缺點是對投影數(shù)據(jù)的完備性要求較高，否則會產(chǎn)生偽影。本文分析了濾波反投影重建算法產(chǎn)生偽影的原因，提出一種新的圖像重建方法，將 傳 統(tǒng)濾波反投影算法與Wiener濾波相結(jié)合。實驗結(jié)果表明，該方法較好地克服了圖像偽影現(xiàn)象，圖像質(zhì)量更好。

1 重建圖像偽影

1.1 偽影概述[4]

CT圖像偽影又稱偽像，是指圖像重建過程中，不同類型的圖像互相干擾和其他各種非隨機干擾在圖像上的表現(xiàn)，使部分影像根本不反映被測物體的對應(yīng)區(qū)域。CT投影是整個CT過程的輸入，投影的好壞直接影響重建圖像的質(zhì)量。若投影數(shù)據(jù)先天不足，則無論何種算法都不能得到完美的重建圖像，得到的只是帶有偽影的圖像，且投影數(shù)據(jù)截斷越大，偽影越多。

造成投影質(zhì)量差的原因主要為[5]幾何誤差、探測器缺陷、射線源不理想以及各種噪聲。幾何誤差因掃描方式而異，對二維CT圖像中心旋轉(zhuǎn)中心不重合，會造成“拖尾”偽影；探測器缺陷和探測器陣列元素響應(yīng)不一致，會造成環(huán)形偽影；射線源的多色性，會導(dǎo)致射線穿透試件時的硬化效應(yīng)，從而造成杯狀偽影；射線源能量歧離也會產(chǎn)生條紋狀偽影；噪聲主要包括隨機噪聲、散射噪聲等。在射線穿過物體時，由于像素對射線采樣的不足，通道信號會產(chǎn)生混疊失真，出現(xiàn)偽影。而濾波反投影重建的本質(zhì)是把取自有限物體空間的射線投影均勻地反投影到射線所及的無限空間的各點之上，因此產(chǎn)生的是星狀偽影。

1.2 偽影原理

本文以反投影重建算法為例來闡述偽影產(chǎn)生機理[6,7]，該算法可等效為以原圖像為輸入，重建后圖像為輸出的成像系統(tǒng)(圖1)。

圖1 反投影重建的等效成像系統(tǒng)Fig.1 Equivalent imaging system of back projection reconstruction.

該系統(tǒng)的點擴展函數(shù)(PSF)可如下求得。設(shè)位于坐標原點的點源δ(x,y)為x-y斷面中唯一像點，并設(shè)掃描方式為平移/旋轉(zhuǎn)。待重建圖像為f(x,y)，射線方向由旋轉(zhuǎn)坐標軸xr與固定坐標軸x的夾角確定，射線透射物體后的投影為pφ(xr)。

掃描坐標系統(tǒng)見圖2，虛線代表射線投影方向，其方向始終與yr平行，當φ角從0°旋轉(zhuǎn)到90°時，得到一系列投影pφi(xr)xr=rcos(θ–φ)。

圖2 平移/旋轉(zhuǎn)掃描方式的坐標系統(tǒng)Fig.2 Coordinate system of tran s lation/rotation scanning mode.

設(shè)輸入圖像為點源，利用 反 投 影重建算法，重建圖像即為該系統(tǒng)的點擴展函數(shù)[5]。

式中，φ0為α(φ)=0的唯一解，令α(φ)=rcos(θ–φ)，則α(φ0)=rcos(θ–φ0)=0。

因r≠0，故sin(θ–φ0)=1。所以

可見，相應(yīng)于反投影重建算法的系統(tǒng)，它的點擴展函數(shù)不是δ-函數(shù)，雖在r=0處能反映原圖是點源的情況，但在r≠0處，像素值不為 0 ，這是此種算法存在偽影的本質(zhì)。

2 Wiener濾波概述

Wiener濾波屬于反卷積算法，是從噪聲中提取信號的一種濾波方法，由Wiener[8]提出，在一、二維信號處理和圖像復(fù)原領(lǐng)域得到廣泛應(yīng)用。在濾波反投影算法(Filtered Back Projection, FBP)對圖像進行重建時，濾波函數(shù)的選擇對消除反投影所產(chǎn)生的星形偽影及圖像重建質(zhì)量影響很大，目前常用R-L濾波函數(shù)和S-L濾波函數(shù)。R-L濾波函數(shù)優(yōu)點是形式簡單、實用，重建圖像輪廓清晰，缺點為重建圖像有明顯的振蕩現(xiàn)象，且當投影數(shù)據(jù)含有噪音時，重建質(zhì)量比較差。S-L濾波函數(shù)優(yōu)點是對投影數(shù)據(jù)的高頻成分有抑制作用，重建圖像的震蕩響應(yīng)小，缺點是圖像分辨率不高。相比而言，Wiener濾波法頻域形式簡單，計算效率高，復(fù)原效果良好，且抗噪性能優(yōu)良，在重建圖像的精度和速度方面都有明顯的優(yōu)勢。

以一線性系統(tǒng)為例來闡述Wiener濾波原理[9]，若它的單位樣本響應(yīng)為h(n)，Wiener濾波器的輸入-輸出關(guān)系如圖3所示。

圖3 Wiener濾波原理示意圖Fig.3 Schematics of the principle of multistage Wiener filter.

由圖3，x(n)為隨機輸入信號，且

其中，s(n)表示信號，v(n)表示噪聲，則輸出y(n)為

y(n)為s(n)的估計值，用 ?(n)表示，即

用e(n)表示真實值與估計值之間的誤差，即

顯然，e(n)是一隨機變量，可能為正或負。因此，選擇用它的均方誤差來表達誤差，所謂均方誤差最小即它的平方的統(tǒng)計期望最小：

3 Wiener?Hopf方程的求解[10]

按最小均方誤差準則確定 Wiener濾波器的沖激響應(yīng)h(n)，令ξ(n)對h(j)的導(dǎo)數(shù)等于零，即得式中，Rxs(m)=E(x(n)s(n+m))是s(n)與x(n)的互相關(guān)函數(shù)，Rxx(m) =E(x(n)x(n+m))是x(n)的自相關(guān)函數(shù)。式(8)所列為 W iener濾波器的標準方程或Wiener-Hopf方程。若已知Rxs(m)和Rxx(m)，則解此方程即可求出Wiener濾波器的沖激響應(yīng)。

設(shè)濾波器沖激響應(yīng)序列的長度為N[11]，沖激響應(yīng)矢量為

濾波器輸入數(shù)據(jù)矢量為

則濾波器的輸出為

Wiener-Hopf 方程也可寫成

其中

式(13)為s(n)與x(n)的互相關(guān)函數(shù)，是N維列矢量；R=E?x(n)xT(n)? 是x(n)的自相關(guān)函數(shù)，是N階方陣。利用求逆矩陣的方法直接求解式( 1 2 )，得

式中，opt表示“最佳”，即FIR Wiener濾波器的沖激響應(yīng)。

利用Wiener濾波器以平滑圖像邊緣，它對于圖像變化越明顯的區(qū)域所起的作用越強。實驗結(jié)果顯示，在程序中加入Wiener濾波后得到圖像更清晰。

4 仿真及實驗結(jié)果

4.1 數(shù)值仿真

為驗證該法的有效性，進行了數(shù)值仿真實驗。本文重建圖像分辨率為 560×420。對三種氣液兩相流流型進行仿真。程序在 Intel(R) Core(TM)2 Duo CPU、1 G內(nèi)存的PC上運行。操作系統(tǒng)為Windows XP系統(tǒng)，程序運行平臺為MAT LAB7.0。由投影數(shù)據(jù)分別對泡狀流、核心流和層流采用不同圖像重建算法進行圖像重建，其中濾波反投影重建算法采用的是R-L濾波函數(shù)，兩種方法迭代次數(shù)均為20。仿真重建圖像如圖4所示。

圖4 FBP與改進的FBP算法重建的圖像Fig.4 Images reconstructed using the FBP and improved FBP algorithm.

4.2 實驗結(jié)果

所用多相流CT實驗裝置的射線源為1.11×1010Bq(300 mCi)241Am源，γ射線能量為59.5 keV，探測器陣列由17個CdZnTe探測器組成，單個探測器的尺寸為3 mm×7 mm×3 mm。將直徑為29、22.4、13.4 mm的有機玻璃棒插入充滿水的流體管道中，模擬三個泡狀體的泡狀流；將Φ60 mm有機玻璃管套在Φ82 mm流體管內(nèi)，將水充入內(nèi)管道模擬核心流；用Φ82 mm有機玻璃管模擬流體管道，在該管道中充入一半水模擬層流。從7個角度采集實驗數(shù)據(jù)，在 M atlab7.0環(huán)境下分別用濾波反投影重建算法和加Wiener濾波的反投影重建算法作圖像重建。

為比較改進前后算法的仿真結(jié)果，用空泡份額及相對誤差對重建圖像的質(zhì)量進行評估?？张莘蓊~是指氣液兩相流中，氣相所占截面積與總流通截面積之比。以流體管道中含有三個氣泡為例計算空泡份額，并與實驗?zāi)Ｐ椭械膶嶋H比例進行比較。根據(jù)已知實驗?zāi)M的流體管道管徑和三個模擬泡狀氣體管徑數(shù)據(jù)，可算出三個模擬泡狀體管體的截面占流體管道截面的比例為Pa=0.1251、Pb=0.0746、Pc=0.0267。用圖像處理技術(shù)算出的各模擬泡狀管截面所成像所占流體管道截面成像的比例用P'a、P'b、P'c表示，空泡份額及相對誤差比較結(jié)果見表1。

濾波反投影重建算法與基于 Wiener濾波的反投影重建算法均采用20次迭代，對泡狀流、核心流和層流進行圖像重建，重建時間見表2。

表1 空泡份額及相對誤差比較Table 1 Comp a r ison of void fra c tio n and relative e rrors.

表2 FBP與改uter time of the FBP算法的重建 oved FBP and impr algo進的 時間Table 2 Comp FBP ns).(20次迭代)rithm (20 iteratio

從實驗結(jié)果對比可見，本文所用改進的重建算法的重建效果明顯好于傳統(tǒng)的濾波反投影的重建效果。在近似相等的時間內(nèi)，濾波反投影重建圖像雖然也可以重建圖像，但表面有許多細條紋和亮紋，有明顯的 偽 影，而基于維納濾波的反投影重建算法重建的圖像則清晰得多，表層均勻，無雜質(zhì)。由表2和表3，加入維納濾波后，重建圖像明顯偽影減少，更接近原始圖像，因此有一定的可行性和應(yīng)用性。

5 結(jié)語

本文提出了一種應(yīng)用于CT成像系統(tǒng)圖像重建算法的改進方案，針對濾波反投影算法產(chǎn)生偽影的問題，通過采取將Wiener濾波技術(shù)融入傳統(tǒng)算法的方法提高重建圖像的質(zhì)量。此方法對有泡狀流、核心流和層流等典型流型的模擬管道進行圖像重建，實驗結(jié)果表明，方法較好地克服了偽影，成像的邊界更加清晰，獲得較好的圖像重建效果。

1 Gilbert B K , Welch B, Holmes D R. Rapid execution of fan beam image reconstruction algorithms using efficient computational techniques and special-purpose[J]. IEEE Trans Biomedical Engineering, 1981, 28(2): 98–115

2 劉杰, 施寅, 阮秋琦. CT快速圖像重建算法研究[J]. 中國醫(yī)學(xué)物理學(xué)雜志, 2003, 3: 15–17 LIU Jie, SHI Yin, RUAN Qiuqi. The study of fast image reconstruction[J]. Chinese Journal of medical Physics,2003, 3: 15–17

3 Lucas G P, Albusaidi K. An axial scanning system for investigating vertical gas-liquid flows in the slug and bubbly-slug transition regimes[J]. Flow Measurement and Instrumentation, 1998, 9: 9–35

4 Luenberger D G. Linear and nonlinear programming[J].MA: Addison-Lesley, 1989, 5: 19–25

5 莊天戈. CT原理與算法[M]. 上海: 上海交通大學(xué)出版社, 1992: 34–35,77 ZHUANG Tiange. Principles andAlgorithms of CT[M].Shanghai: Shanghai Jiaotong UniversityPress, 1992:34–35,77

6 Fossa M. Design and performance of a conductance probe for measuring the liquid fractionin two-phase gas-liquid flows[J]. Flow Measurement and Instrumentation, 1998, 9:103–109

7 Thorn R, Hamme E A, Johanen G A. Recent developments inhere-phase flow measurement[J].Measurement Science Technology, 1997, 8(7): 691–701

8 Klaus Muelle. Anti-aliased three-dimentional cone-beam reconstruction of low-constrast objects with algebraic methods[J]. IEEE Transactiol Imagining, 1999, 5(2):116–134

9Alexander D P, Zayed M R. Adaptive filtering primer with MATLAB[J]. Chinese Journal of medical Physics,2006, 10: 50–75

10 Johansen G A, Froystein T, Hjertaker B T. A dual sensor flow imaging topographic system[J]. Measurement Sci Tech, 1996, 7(2): 297–307

11 Jiang Hsieh. Computer tomography principle, design,artifacts and recent advances[J]. SPIE of America, 2003,5(10): 10–12