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

馬 敏 趙 亮 姬晶晶

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

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

1 重建圖像偽影

1.1 偽影概述[4]

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

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

1.2 偽影原理

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

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

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

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

圖2 平移/旋轉掃描方式的坐標系統Fig.2 Coordinate system of tran s lation/rotation scanning mode.

設輸入圖像為點源，利用 反 投 影重建算法，重建圖像即為該系統的點擴展函數[5]。

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

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

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

2 Wiener濾波概述

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

以一線性系統為例來闡述Wiener濾波原理[9]，若它的單位樣本響應為h(n)，Wiener濾波器的輸入-輸出關系如圖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)是一隨機變量，可能為正或負。因此，選擇用它的均方誤差來表達誤差，所謂均方誤差最小即它的平方的統計期望最小：

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

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

設濾波器沖激響應序列的長度為N[11]，沖激響應矢量為

濾波器輸入數據矢量為

則濾波器的輸出為

Wiener-Hopf 方程也可寫成

其中

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

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

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

4 仿真及實驗結果

4.1 數值仿真

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

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

4.2 實驗結果

所用多相流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流體管內，將水充入內管道模擬核心流；用Φ82 mm有機玻璃管模擬流體管道，在該管道中充入一半水模擬層流。從7個角度采集實驗數據，在 M atlab7.0環境下分別用濾波反投影重建算法和加Wiener濾波的反投影重建算法作圖像重建。

為比較改進前后算法的仿真結果，用空泡份額及相對誤差對重建圖像的質量進行評估。空泡份額是指氣液兩相流中，氣相所占截面積與總流通截面積之比。以流體管道中含有三個氣泡為例計算空泡份額，并與實驗模型中的實際比例進行比較。根據已知實驗模擬的流體管道管徑和三個模擬泡狀氣體管徑數據，可算出三個模擬泡狀體管體的截面占流體管道截面的比例為Pa=0.1251、Pb=0.0746、Pc=0.0267。用圖像處理技術算出的各模擬泡狀管截面所成像所占流體管道截面成像的比例用P'a、P'b、P'c表示，空泡份額及相對誤差比較結果見表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

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

5 結語

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

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]. 中國醫學物理學雜志, 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]. 上海: 上海交通大學出版社, 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