斜變?yōu)V波器時(shí)域加窗對(duì)圖像重建的影響

喬志偉

（中北大學(xué)電子與計(jì)算機(jī)科學(xué)技術(shù)學(xué)院，山西太原 030051）

斜變?yōu)V波器時(shí)域加窗對(duì)圖像重建的影響

喬志偉

（中北大學(xué)電子與計(jì)算機(jī)科學(xué)技術(shù)學(xué)院，山西太原 030051）

提出了斜變?yōu)V波器的單位脈沖響應(yīng)的時(shí)域加窗方法。斜變?yōu)V波器的單位脈沖響應(yīng)應(yīng)取奇數(shù)個(gè)點(diǎn)并偶對(duì)稱，以保證其嚴(yán)格零相位。可以用來(lái)截取單位脈沖響應(yīng)的窗函數(shù)有矩形窗、三角窗、漢寧窗、海明窗、布萊克曼窗等。仿真實(shí)驗(yàn)表明，矩形窗精度最好，三角窗因過(guò)度地壓低了第一旁瓣，精度最差，其余窗對(duì)應(yīng)的CT圖像的精度介于兩者之間。在對(duì)斜變?yōu)V波器的單位脈沖響應(yīng)截?cái)鄷r(shí)應(yīng)該使用矩形窗。

計(jì)量學(xué)；斜變?yōu)V波器；窗函數(shù)；圖像重建；濾波反投影算法

1 引 言

CT（Computed Tomography）即計(jì)算機(jī)斷層成像技術(shù)是最好的無(wú)損檢測(cè)手段之一，已經(jīng)廣泛應(yīng)用到醫(yī)學(xué)、工業(yè)、地震學(xué)和考古學(xué)等領(lǐng)域［1，2］。X射線工業(yè)CT算法分為解析式算法和代數(shù)式算法2種，其中解析式算法居于優(yōu)勢(shì)地位［3，4］。在解析式重建算法中，又分為二維CT算法和三維CT算法，二維算法的經(jīng)典代表是濾波反投影（filtered back projection，F(xiàn)BP）算法，而三維算法的經(jīng)典代表是FDK算法［5，6］。從本質(zhì)上說(shuō)，由于FDK算法就是一種近似的濾波反投影算法，因此要研究影響此類算法的各種因素時(shí)，往往以FBP算法為研究對(duì)象。

在濾波反投影算法中，斜變?yōu)V波器的設(shè)計(jì)方法為【7，8】：

（1）通過(guò)對(duì)理想斜變?yōu)V波器ω加窗形成實(shí)用濾波器的頻率響應(yīng)H（ω）。加矩形窗的叫R-L濾波器；加sinc窗的叫S-L濾波器。

（2）由H（ω）計(jì)算得到單位沖激響應(yīng)h（t）。

（3）根據(jù)沖激響應(yīng)不變法，計(jì)算得到單位脈沖響應(yīng)h（n）。

（4）根據(jù)投影信號(hào)和濾波投影信號(hào)均為有限長(zhǎng)度的特點(diǎn)，截取h（n）為有限長(zhǎng)度的單位脈沖響應(yīng)。

在經(jīng)典的FIR濾波器的窗函數(shù)法設(shè)計(jì)中，要用不同于矩形窗的更加平滑的窗去截取有限長(zhǎng)度，以減少頻譜泄露。對(duì)于上述第（4）步，各種文獻(xiàn)均未見(jiàn)論及加窗的問(wèn)題，即相當(dāng)于均用矩形窗截?cái)嗟姆椒?。而根?jù)信號(hào)處理的基本理論可知，加三角窗、漢寧窗、海明窗等可以有效減少頻譜泄露。因此，有必要研究時(shí)域加窗對(duì)圖像重建效果的影響。本文以S-L濾波器為例研究此問(wèn)題。

2 斜變?yōu)V波器單位脈沖響應(yīng)長(zhǎng)度確定

根據(jù)離散時(shí)間信號(hào)卷積理論，2個(gè)有限長(zhǎng)的信號(hào)做卷積，則卷積和信號(hào)的長(zhǎng)度為2個(gè)信號(hào)的長(zhǎng)度之和減1。

設(shè)離散的投影信號(hào)p（n）的取值范圍為［－N，N］，S-L濾波器的單位脈沖響應(yīng)h（n）的取值范圍為［－∞，∞］，則離散的濾波投影信號(hào)g（n）的取值范圍為［－∞，∞］。但在實(shí)際的CT工程實(shí)踐中，濾波投影信號(hào)沒(méi)有必要取無(wú)限長(zhǎng)，而只需在［－N，N］這個(gè)范圍內(nèi)取值［9，10］。

根據(jù)卷積的定義，可知g（n）在2個(gè)邊界點(diǎn)－N和N的值為

可見(jiàn)，h（n）的取值范圍為［－2N，2N］，共4N＋1個(gè)點(diǎn)，就可以求得所需范圍內(nèi)的濾波投影信號(hào)g（n）。

總之，若投影信號(hào)和濾波投影信號(hào)的長(zhǎng)度均取M點(diǎn)，則斜變?yōu)V波器的單位脈沖響應(yīng)h（n）只需要取2M－1點(diǎn)。

需要注意的是，截取后的斜變?yōu)V波器h（n）必須是嚴(yán)格零相位的，因此，在截取2M－1個(gè)點(diǎn)時(shí)，必須以0為對(duì)稱中心，即取值范圍應(yīng)該為［－（M－1），M－1］。

3 時(shí)域加窗方法

如果直接取［－（M－1），M－1］范圍的h（n）作為濾波器的單位脈沖響應(yīng)，那么這種方式相當(dāng)于在無(wú)限長(zhǎng)的h（n）之上加了1個(gè)以0為中心的奇數(shù)個(gè)點(diǎn)的矩形窗q（n）。

矩形窗是最簡(jiǎn)單的一種截取無(wú)限長(zhǎng)單位脈沖響應(yīng)為有限長(zhǎng)單位脈沖響應(yīng)的窗函數(shù)，因?yàn)槭峭蝗唤財(cái)嗟?，?huì)引起強(qiáng)烈的頻譜泄露。為了減小頻譜泄露，就需要加一些不突然截?cái)嗟拇昂瘮?shù)。這些窗函數(shù)有三角窗、漢寧窗、海明窗、布萊克曼窗等，其性能參數(shù)見(jiàn)表1。－80

表1 6種窗函數(shù)的基本參數(shù)比較

顯然，任何一種窗函數(shù)均不可能做到過(guò)渡帶最窄，同時(shí)阻帶衰減最大。過(guò)渡帶寬度和阻帶衰減程度是2個(gè)矛盾的參數(shù)，在工程中一般是根據(jù)對(duì)過(guò)渡帶寬和阻帶衰減的要求來(lái)選擇窗函數(shù)。

可以認(rèn)為表1中的后5種窗函數(shù)均為對(duì)矩形窗函數(shù)的優(yōu)化和校正，其中漢寧窗、海明窗、布萊克曼窗其實(shí)均為余弦類的窗函數(shù)，只是各自性能有差異，而凱澤窗的波形也與海明窗近似。

在FIR濾波器的設(shè)計(jì)理論中，因?yàn)镕IR濾波器的單位脈沖響應(yīng)是因果序列，窗函數(shù)是加在［0，N－1］范圍內(nèi)的。而在對(duì)斜變?yōu)V波器的單位脈沖響應(yīng)截取時(shí)，信號(hào)要偶對(duì)稱，且選取奇數(shù)個(gè)點(diǎn)，故窗函數(shù)的表示方法與FIR濾波器設(shè)計(jì)理論中的形式不同。

為了表示方便，現(xiàn)假設(shè)要對(duì)h（n）在［－N，N］之間截取，則相應(yīng)的窗函數(shù)為

（1）三角窗

三角窗變化比較劇烈，呈線性變化，而漢寧窗、海明窗和布萊克曼窗因均為升余弦窗，其變化較平滑。

4 仿真實(shí)驗(yàn)及結(jié)果分析

4.1 仿真模型及實(shí)驗(yàn)條件

采用濾波反投影算法對(duì)一個(gè)由4個(gè)細(xì)圓柱插到1個(gè)粗圓柱的模型的某斷層進(jìn)行重建。探測(cè)器分別在每個(gè)角度采集101點(diǎn)（采樣間隔為0.1 cm），共采集180個(gè)角度的投影（角度步進(jìn)間隔為1°），斜變?yōu)V波器采用S-L濾波器；采用卷積運(yùn)算方式濾波，像素驅(qū)動(dòng)方式為反投影。重建的圖像為以旋轉(zhuǎn)中心為圖像中心的101×101的圖像。

模型及其采集示意圖見(jiàn)圖1，模型參數(shù)見(jiàn)表2。

圖1 實(shí)驗(yàn)用仿真模型及0°采集示意圖

表2 模型參數(shù)

對(duì)S-L濾波器的h（n）采用5種不同的窗函數(shù)截?cái)?，研究不同窗函?shù)對(duì)重建精度的影響。

對(duì)重建圖像的質(zhì)量評(píng)估采用歸一化均方距離判據(jù)d和歸一化平均絕對(duì)距離判據(jù)r這2個(gè)參數(shù)來(lái)評(píng)估。誤差的公式分別為

4.2 實(shí)驗(yàn)結(jié)果

窗函數(shù)對(duì)重建精度的影響見(jiàn)表3，重建的圖像見(jiàn)圖2。

表3 S-L濾波器時(shí)域不同窗函數(shù)對(duì)重建精度影響

圖2 時(shí)域加窗重建的圖像對(duì)比

4.3 結(jié)果分析

由表3可以看出，按照精度從高到低的順序?yàn)椋壕匦未?、海明窗、漢寧窗、布萊克曼窗、三角窗。

由圖2可以看出，只有三角窗對(duì)應(yīng)的CT圖像失真明顯，與其它幾個(gè)區(qū)別較大，而其它幾個(gè)重建圖像的區(qū)別用肉眼不好分辨。三角窗的失真主要表現(xiàn)在中心圓變亮，同時(shí)中心圓周圍的區(qū)域也明顯變亮。

表3和圖2反映出來(lái)的規(guī)律一致，能得出共同的結(jié)論。只是這個(gè)結(jié)論應(yīng)該引起特別的注意，因?yàn)檫@里的加窗性能與FIR濾波器設(shè)計(jì)方法中的加窗性能不一致。

在FIR濾波器的設(shè)計(jì)中，為了減小吉布斯效應(yīng)，應(yīng)該使用過(guò)渡平滑的除矩形窗以外的其它濾波器。在犧牲過(guò)渡帶寬不大的情況下，阻帶性能從好到差的順序?yàn)椴既R克曼窗、海明窗、漢寧窗、三角窗、矩形窗。一般在工程中認(rèn)為，這個(gè)順序是一個(gè)性能由高到低的順序。

然而，旁瓣衰減迅速（阻帶衰減性能好）的窗函數(shù)應(yīng)用到圖像重建中，圖像精度并沒(méi)有提高。

采用過(guò)渡平滑的窗函數(shù)的優(yōu)點(diǎn)是可以使得吉布斯效應(yīng)減小，然而它的一個(gè)顯著的缺點(diǎn)是窗函數(shù)的形狀改變了信號(hào)時(shí)域值的真實(shí)表示。平滑可以減小吉布斯效應(yīng)，從而提高精度；另一方面，加窗引起了信號(hào)時(shí)域的失真，從而降低精度，這是相互制約的一對(duì)矛盾。

當(dāng)因?yàn)榧哟耙鸬木冉档统潭瘸^(guò)了平滑引起的精度提高的程度，總的加窗精度就要比不加窗（加矩形窗）的精度要低，加窗就沒(méi)有意義了。

由表3可以看出，矩形窗的精度最高，三角窗的精度最低，其它升余弦類的窗函數(shù)的重建精度介于其間。這說(shuō)明在對(duì)斜變?yōu)V波器的單位脈沖響應(yīng)時(shí)域加窗時(shí)，加窗引起的精度降低程度超過(guò)了平滑過(guò)渡引起的精度提高程度?？偟膩?lái)說(shuō)，加窗反而引起了精度的降低。

這種現(xiàn)象從h（n）的特征也可以解釋。h（n）的第一旁瓣特別重要，它的大小對(duì)重建圖像的精度影響重大。而無(wú)論是三角窗還是海明窗等，均對(duì)第一旁瓣做了較大的衰減，從而引起較大的精度損失。圖3給出了加矩形窗的S-L濾波器的單位脈沖響應(yīng)與加三角窗的S-L濾波器的單位脈沖響應(yīng)的對(duì)比。

圖3n?。郏?00，100］h（n）在［2，5］范圍的波形

由圖3可以看出，加三角窗后的單位脈沖響應(yīng)絕對(duì)值變小，在［2，5］范圍內(nèi)三角窗的波形整體提升到了矩形窗對(duì)應(yīng)的波形上面，這相當(dāng)于對(duì)每一個(gè)旁瓣的幅值均做了壓縮。而這種壓縮必然導(dǎo)致較大的誤差，所以三角窗形式的斜變?yōu)V波器重建出來(lái)的CT圖像誤差較大。而其它的余弦類的窗函數(shù)也均有不同程度地改變了旁瓣幅值，導(dǎo)致了精度的降低。

可見(jiàn)，在斜變?yōu)V波器的設(shè)計(jì)中，精度最好的加窗方式是矩形窗，沒(méi)有必要加起平滑作用的窗函數(shù)。

5 結(jié) 論

在用FIR濾波器的窗函數(shù)法設(shè)計(jì)濾波器時(shí)，根據(jù)過(guò)渡帶寬和阻帶衰減速度的要求，可以適當(dāng)選擇三角窗、漢寧窗、海明窗等不是突然截?cái)嗟拇昂瘮?shù)，以減少吉布斯效應(yīng)。

但是，在用窗函數(shù)截取無(wú)限長(zhǎng)的單位脈沖響應(yīng)為有限長(zhǎng)序列時(shí)，只能用矩形窗，而不能用三角窗和升余弦類窗。其原因在于，斜變?yōu)V波器的單位脈沖響應(yīng)的第一旁瓣對(duì)重建精度有顯著影響，非矩形窗使得第一旁瓣的幅值被壓縮，從而增大了誤差。

總之，在設(shè)計(jì)濾波反投影算法的斜變?yōu)V波器的單位脈沖響應(yīng)時(shí)，只能用矩形窗。

［1］ 王召巴.任意角度入射的三維CT投影方法［J］.計(jì)量學(xué)報(bào)，2002，23（1）：1－5.

［2］ 張朝宗，郭志平，張朋，等.工業(yè)CT技術(shù)和原理［M］.北京：科學(xué)出版社，2009.

［3］ 張順利，張定華，趙歆波.一種射束與像素的快速遍歷和求交算法［J］.中國(guó)圖象圖形學(xué)報(bào)，2009，14（10）：1961－1965.

［4］ 王東明，盧虹冰，張軍英.基于統(tǒng)計(jì)特性的小波噪聲抑制在低劑量CT中的應(yīng)用［J］.中國(guó)圖象圖形學(xué)報(bào)，2008，13（05）：876－881.

［5］ 馬建華，陳武凡.不含旋轉(zhuǎn)角度微分的螺旋錐束CT重建［J］.中國(guó)圖象圖形學(xué)報(bào)，2008，13（04）：647－653.

［6］ 張全紅，路宏年，楊民.錐束工業(yè)CT中Feldkamp重建算法的快速實(shí)現(xiàn)［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2006，27（6）：931－933.

［7］ 喬志偉，魏學(xué)業(yè)，韓焱.解析法圖像重建中的插值技術(shù)研究［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2009，30（9）：2213－2216.

［8］ 曾更生.醫(yī)學(xué)圖像重建［M］.北京：高等教育出版社，2010.

［9］ 李保磊，傅健，魏東波，等.工業(yè)計(jì)算機(jī)斷層成像系統(tǒng)轉(zhuǎn)臺(tái)旋轉(zhuǎn)中心的確定［J］.航空動(dòng)力學(xué)報(bào)，2009，24（07）：1544－1548.

［10］ 李保磊，傅健，黃巧珍，等.一種基于正弦圖的工業(yè)CT系統(tǒng)轉(zhuǎn)臺(tái)旋轉(zhuǎn)中心自動(dòng)確定方法［J］.航空學(xué)報(bào)，2009，30（07）：1341－1345.

Im pact of Time Domain W indow-adding for the Ramp Filter on Image Reconstruction

QIAO Zhi-wei
（School of Electronics and Computer Science and Technology，North University of China，Taiyuan，Shanxi030051，China）

The time domain window-adding method for the unit impulse response of the ramp filter is proposed.the unit impulse response of the ramp filter should be symmetrical and its length should be odd number to ensure that it is strict zero-phase.There are many windows to truncate the unit impulse response，for example，rectangle window，triangle window，Hanning window，Hamming windows，Blackman window etc.The simulation experiment demonstrates that the rectangle window is the best and the triangle window is the worst，for the triangle window depresses the first side petal overly.The precision of the CT images reconstructed by the otherswindows is in-between.The rectanglewindow should be used when the unit impulse response of the ramp filter is been truncated.

Metrology；Ramp filter；Window function；Image reconstruction；Filtered back projection algorithm

TB973

A

1000-1158（2014）03-0268-04

10．3969／j．issn．1000-1158．2014．03．15

2012-01-16；

2013-08-30

國(guó)家自然科學(xué)基金（61071193，61171179）

喬志偉（1977－），男，山西洪洞人，中北大學(xué)副教授，博士，主要研究方向?yàn)樾畔@取及處理技術(shù)。zhiweiqiaook＠nuc.edu.cn