Kaiser窗在濾波反投影圖像重建中的應(yīng)用研究?

于尚民

（青島大學(xué)電子信息學(xué)院 青島 266071）

1 引言

重建算法是CT（Computed Tomography，計(jì)算機(jī)斷層成像）圖像重建技術(shù)的核心。迭代法和解析法是目前最常見的兩大類重建算法。濾波反投影（Filter Back Projection，F(xiàn)BP）算法是一種最為常用的解析圖像重建算法［1］。濾波反投影算法重建速度快、重建圖像精度高，因而它是一種應(yīng)用十分廣泛的CT圖像重建算法。在直接反投影算法中，重建后的圖像是由所有方向下的投影數(shù)據(jù)累加而形成的。在累加過程中，無限空間中原本像素值為零的點(diǎn)的像素值被均勻地回抹為有限物體空間的投影值，即其值變?yōu)榉橇恪＿@種非零的像素值導(dǎo)致了明顯的“星”狀偽影。由投影定理可知，在反投影重建之前將投影數(shù)據(jù)用濾波函數(shù)進(jìn)行濾波可消除“星”狀偽影［2］。因此，濾波反投影算法的關(guān)鍵在于濾波函數(shù)的設(shè)計(jì)。窗函數(shù)法常用于濾波函數(shù)的設(shè)計(jì)。R-L（Ram-Lak）濾波函數(shù)和 S-L（Shepp-Lo?gan）濾波函數(shù)是濾波反投影算法中典型的兩種基于傳統(tǒng)窗函數(shù)的濾波函數(shù)。其中R-L濾波函數(shù)是一種簡(jiǎn)單易行的濾波函數(shù)。它的重建圖像質(zhì)量較高，但卻伴隨著嚴(yán)重的吉布斯（Gibbs）現(xiàn)象［3］。Kai?ser窗是一種可由參數(shù)調(diào)節(jié)的窗函數(shù)，有著廣泛的應(yīng)用空間。本文將Kaiser窗濾波函數(shù)與基于傳統(tǒng)窗函數(shù)的R-L濾波函數(shù)和S-L濾波函數(shù)進(jìn)行對(duì)比，經(jīng)實(shí)驗(yàn)驗(yàn)證，Kaiser窗濾波函的重建圖像精度更高、抗噪能力更強(qiáng)。

2 濾波反投影算法實(shí)現(xiàn)原理

根據(jù)投影定理，在反投影重建之前將投影數(shù)據(jù)進(jìn)行濾波，對(duì)濾波后的數(shù)據(jù)再進(jìn)行反投影運(yùn)算，可消除“星”狀偽影。其具體步驟如下：

1）設(shè) Sθ(ω)為某角度下的投影數(shù)據(jù) pθ(t)的一維傅里葉變換；

2）用一維權(quán)重因子 ||ρ乘以Sθ(ω)；

3）對(duì)第2）步所得結(jié)果作一維傅里葉逆變換，記為Qθ(t)；

4）對(duì) 0°～180°內(nèi)的全部修正過的 Qθ(t)作直接反投影運(yùn)算，得到重建圖像 f(x，y)。

3 濾波器的設(shè)計(jì)思想

濾波器設(shè)計(jì)的好壞會(huì)對(duì)重建圖像的精度產(chǎn)生不可忽視的影響。FBP算法理論上要求濾波器的系統(tǒng)函數(shù)滿足HR(ρ)= ||ρ，但現(xiàn)實(shí)中的濾波函數(shù)的頻帶寬度是無窮大的。在實(shí)際應(yīng)用中，由Paley-Wiener準(zhǔn)則可知，這種頻帶寬度無窮大的濾波器是無法實(shí)現(xiàn)的。但由于實(shí)際CT設(shè)備的分辨率是有限的，這使其收集到的投影數(shù)據(jù)是一種頻譜能量分布于低頻區(qū)域的密度分布［4］。因此可以忽略超出某個(gè)截止頻率的投影數(shù)據(jù)，即

式（4）中的W(ρ)稱為窗函數(shù)。W(ρ)的好壞對(duì)重建圖像精度起著極其重要的作用。因此要求窗函數(shù)W(ρ)應(yīng)遵守下列原則：

1）主瓣寬度要窄，提高識(shí)別精度，以增大圖像分辨率；

2）最大旁瓣幅值盡可能的要小，增大阻帶衰減、削弱Gibbs效應(yīng)。

3.1 常用濾波器

3.1.1 R-L濾波器

R-L濾波函數(shù)頻率響應(yīng)為

其中：

與濾波函數(shù)HR-L(ρ)相對(duì)應(yīng)的卷積函數(shù)沖激響應(yīng)hR-L(R)為

其中B=1/(2d)，sinc(x)=sin(x)/x。

其離散式（采樣序列）為

R-L濾波函數(shù)由于沒有復(fù)雜的三角函數(shù)計(jì)算，因而運(yùn)算速度較快。其所得重建圖像較為清晰，但卻伴隨有嚴(yán)重的Gibbs現(xiàn)象，并對(duì)含有噪聲的投影數(shù)據(jù)重建效果較差。

3.1.2 S-L濾波器

S-L濾波函數(shù)緩解了R-L濾波函數(shù)的震蕩，其頻率響應(yīng)為

其離散表達(dá)式為

S-L濾波函數(shù)通過對(duì)投影數(shù)據(jù)中的高頻成分的抑制，來減小Gibbs現(xiàn)象。同時(shí)S-L濾波器有著較強(qiáng)的抗噪能力［5］。但在低頻段，其重建圖像的質(zhì)量不及R-L濾波函數(shù)，這是由于hS-L(ρ)在低頻段偏離 ||ρ的原因。

除此之外，常用濾波函數(shù)還有三角窗濾波函數(shù)、Hanning窗濾波函數(shù)、Kaiser窗濾波函數(shù)等［6］。

3.2 圖像評(píng)價(jià)指標(biāo)

多數(shù)情況下肉眼難以分辨圖像質(zhì)量，因而需要引入量化評(píng)判判據(jù)，其表達(dá)式如下：

1）歸一化均方距離d：

其中，N*N為圖像的像素，tˉ為實(shí)際物體模型密度的平均值，重建圖像和物體模型中第u行、v列的像素密度分別用ru，v和tu，v表示。d的值越大，表示CT重建圖像和實(shí)際物體的偏差越大，圖像越不清晰。d值受少數(shù)偏差大的像素點(diǎn)的影響較大。

2）歸一化平均絕對(duì)距離r，即：

歸一化平均絕對(duì)距離r與歸一化均方距離d相比受少數(shù)偏差大的像素點(diǎn)影響較小。歸一化平均絕對(duì)距離r反映了大多數(shù)像素點(diǎn)存在誤差的情形。

依據(jù)歸一化均方距離d和歸一化平均絕對(duì)距離r兩項(xiàng)指標(biāo)，對(duì)矩形窗、S-L窗、Hanning窗、Ham?ming窗這四種窗函數(shù)的重建圖像誤差進(jìn)行比較，如表1所示。

表1 窗函數(shù)濾波器的誤差比較

由表1可知，加Hanning窗重建圖像的歸一化均方距離d和歸一化平均絕對(duì)距離r要比其他三種窗函數(shù)要小，重建圖像精度較高。

4 Kaiser濾波函數(shù)

4.1 Kaiser窗函數(shù)性質(zhì)分析

Kaiser窗函數(shù)是一種常用的窗函數(shù)，其表達(dá)式為

其中，Kaiser濾波函數(shù)時(shí)域內(nèi)的主瓣寬度和旁瓣高低可以通過參數(shù)β來調(diào)節(jié)，β值越大，主瓣的寬度越大，旁瓣幅值越小，I0(x)為零階第一類修正貝塞爾（Bessel）函數(shù)，其級(jí)數(shù)形式為

表2是Kaiser窗函數(shù)在不同β值下的性能。

表2 Kaiser窗函數(shù)在不同β值下的性能

各窗函數(shù)（長(zhǎng)度為64）的主瓣寬度和第一旁瓣衰減值如表3所列。

表3 主瓣寬度與第一旁瓣相對(duì)主瓣衰減值

由和表3可以看出，當(dāng)β=6.5、長(zhǎng)度為64時(shí)，Kaiser窗的主瓣寬度和Hanning窗主瓣寬度一樣，但旁瓣幅值又比Hanning窗和矩形窗的小。根據(jù)濾波函數(shù)的設(shè)計(jì)原理和圖像重建誤差的影響因素可知，主瓣寬度越小，識(shí)別精度越高，圖像重建的精度越高；旁瓣幅值越小，Gibbs效應(yīng)越弱。因此，可以發(fā)現(xiàn)Kaiser窗濾波函數(shù)的空間分辨率、密度分辨率均較高，其圖像重建的精度較高。

4.2 Kaiser窗函數(shù)實(shí)驗(yàn)結(jié)果檢驗(yàn)與分析

實(shí)驗(yàn)中采用了二維平行束FBP算法進(jìn)行計(jì)算機(jī)實(shí)驗(yàn)仿真來驗(yàn)證Kaiser窗濾波函數(shù)的圖像重建精度。所用的投影數(shù)據(jù)由Shepp-Logan模型生成，其大小為512*512，重建圖像的大小為514*514，三種濾波函數(shù)的重建圖像如圖1所示。

圖1 三種濾波函數(shù)的重建結(jié)果

由表4可知，Kaiser窗（β=6.5）濾波函數(shù)重建的圖像，其歸一化絕對(duì)距離r和歸一化均方距離d分別為0.1372和0.32537比R-L濾波函數(shù)和S-L濾波函數(shù)要小，因而圖像重建精度較高。

表4 三種濾波函數(shù)重建的d和r值

4.3 Kaiser窗在噪聲環(huán)境下的重建結(jié)果分析

在實(shí)際CT數(shù)據(jù)的采集中，投影數(shù)據(jù)常常混雜有噪聲。因此，本文引入了泊松噪聲來對(duì)Kaiser窗濾波函數(shù)的重建圖像精度進(jìn)行分析。

改變參數(shù)λ，就可以改變產(chǎn)生泊松噪聲的大小。實(shí)驗(yàn)中選取大小為512*512的Shepp-Logan模型的投影數(shù)據(jù)。噪聲測(cè)試中λ的取值分別為0.1，1，2。其重建圖像如圖2～圖4所示。

圖2 λ=0.1時(shí)Kaiser窗濾波函數(shù)的重建結(jié)果

圖3 λ=1時(shí)Kaiser窗濾波函數(shù)的重建結(jié)果

圖4 λ=2時(shí)Kaiser窗濾波函數(shù)的重建結(jié)果

由圖2～圖4可以看出，Kaiser窗濾波函數(shù)在泊松噪聲下仍然可以獲得高精度的重建圖像。

圖5是對(duì)大小為512*512的Shepp-Logan模型的投影數(shù)據(jù)，以參數(shù)β為自變量，分別以歸一化均方距離d和歸一化絕對(duì)距離r為因變量，分析在λ=0.1、λ=1和λ=2的噪聲情況下Kaiser窗濾波函數(shù)的β值對(duì)圖像重建精度的影響曲線圖。

圖5 評(píng)價(jià)函數(shù)變化曲線

由圖5可以看出，在λ=0.1、λ=1和λ=2時(shí)，隨著β的增大歸一化平均絕對(duì)距離r和歸一化均方距離d均呈先減小后增大的變化趨勢(shì)。這是由于，當(dāng)β增大時(shí)，Kaiser窗的旁瓣幅值減小，Gibbs效應(yīng)減弱，同時(shí)主瓣寬度增加，識(shí)別精度降低。在曲線的前半段，旁瓣幅值減小的影響大于主瓣寬度增加的影響，故d和r減??；在曲線后半段，主瓣寬度增加的影響大于旁瓣幅值減小的影響，故d和r增大。對(duì)λ=1的d和r變化曲線分別進(jìn)行擬合可得：

對(duì)式（16）、（17）求一階導(dǎo)數(shù)可知，當(dāng)β=12.2658時(shí)，d有最小值0.3265；當(dāng)β=22.3863時(shí)，r有最小值0.1629，這與實(shí)驗(yàn)結(jié)果較為吻合。

經(jīng)實(shí)驗(yàn)驗(yàn)證，一般當(dāng)β≈0時(shí)，Kaiser窗濾波函數(shù)近似于R-L濾波函數(shù)；當(dāng)β≈4時(shí)，Kaiser窗濾波函數(shù)近似于S-L濾波函數(shù)；當(dāng)β≈5.15時(shí)，Kaiser窗濾波函數(shù)近似于Hamming窗濾波函數(shù)。當(dāng)β值取得過大時(shí)，會(huì)使重建圖像出現(xiàn)模糊。在實(shí)際應(yīng)用中參數(shù)β的取值取決于多個(gè)方面的因素，例如CT投影數(shù)據(jù)的采集設(shè)備、投影數(shù)據(jù)中混有的噪聲種類、原始物體模型的種類以及評(píng)價(jià)函數(shù)等。

5 結(jié)語

本文通過對(duì)Kaiser窗濾波函數(shù)重建圖像的分析，計(jì)算了其歸一化絕對(duì)距離r和歸一化均方距離d，證明了Kaiser窗濾波函數(shù)的圖像重建精度比基于傳統(tǒng)窗函數(shù)的R-L、S-L濾波函數(shù)的更高。此外本文還分析了泊松噪聲對(duì)Kaiser窗濾波函數(shù)重建圖像的影響，驗(yàn)證了Kaiser窗濾波函數(shù)的抗噪能力，并通過曲線擬合的方法來確定最優(yōu)的β值，證實(shí)了Kaiser窗濾波函數(shù)的優(yōu)越性。