洪振厚，周 彬，郭金川

深圳大學(xué)光電工程學(xué)院，光電子器件與系統(tǒng)教育部/廣東省重點實驗室，廣東深圳518060

【光電工程 / Optoelectronic Engineering】

錐束X射線CT圖像重建的新型濾波函數(shù)

洪振厚，周 彬，郭金川

深圳大學(xué)光電工程學(xué)院，光電子器件與系統(tǒng)教育部/廣東省重點實驗室，廣東深圳518060

針對X射線計算機(jī)斷層成像(computed tomography，CT)圖像重建FDK算法中，采用通常的濾波函數(shù)會導(dǎo)致明顯的Gibbs現(xiàn)象，影響重建圖像的質(zhì)量．基于混合濾波和加權(quán)平均理論，設(shè)計了一種新型的NEW-MS-L濾波函數(shù)．先將S-L濾波函數(shù)加權(quán)平均為M3S-L濾波函數(shù)，再與NEW濾波函數(shù)疊加,得到NEW-MS-L濾波函數(shù)．通過數(shù)值仿真，分別對比NEW、R-L-S-L、R-L-NEW、R-L-MS-L和NEW-MS-L濾波函數(shù)的重建圖像效果，結(jié)果表明，NEW-MS-L濾波函數(shù)能夠在保持較高圖像分辨率的情況下，更有效地抑制Gibbs現(xiàn)象，使重建圖像效果更佳．

光學(xué)工程；計算機(jī)斷層成像；FDK算法；濾波函數(shù)；Gibbs現(xiàn)象；混合濾波；加權(quán)平均

隨著計算機(jī)斷層成像(computed tomography, CT)在醫(yī)療診斷、工業(yè)無損探傷和食品安全檢測等領(lǐng)域廣泛的應(yīng)用[1-5]，其對圖片質(zhì)量的要求也越來越高，進(jìn)而對軟硬件特別是算法的要求也越來越高．在眾多算法中，F(xiàn)DK算法[6]應(yīng)用最為普遍．FDK算法是由Feldkamp、Davis和Kress提出的一種基于圓軌道掃描的近似重建算法，其本質(zhì)是濾波反投影(filtered back projection，F(xiàn)BP)，它對圖像質(zhì)量的影響非常明顯．傳統(tǒng)的濾波函數(shù)能保持較高的空間分辨率，但同時引起明顯的Gibbs現(xiàn)象(即有明顯的振蕩效應(yīng))，致使重建效果較差．因此，在保持高圖像分辨率的情況下設(shè)計新的濾波函數(shù)，降低Gibbs現(xiàn)象就尤為重要．為此，研究人員相繼提出NEW濾波函數(shù)[7]、R-L-S-L、R-L-NEW和R-L-MS-L等混合濾波函數(shù)[8-10]來消除Gibbs現(xiàn)象．雖然這些濾波函數(shù)都具有保持高空間分辨率的同時減小Gibbs現(xiàn)象的作用，但無法完全消除Gibbs現(xiàn)象，主要原因是與之結(jié)合的R-L濾波函數(shù)的近旁瓣突出，遠(yuǎn)旁瓣的幅度和寬度較大[9-10]，而NEW濾波函數(shù)和R-L濾波函數(shù)一樣，由理想濾波器推導(dǎo)出，所以其近旁瓣也較為突出[7].本研究基于FDK算法，設(shè)計了一種新型濾波函數(shù)進(jìn)行圖像重建，在保持較高空間分辨率的情況下大幅減小Gibbs現(xiàn)象，提高了密度分辨率，進(jìn)而改進(jìn)了重建圖像質(zhì)量．

1 FDK算法介紹

FDK算法因其高效、簡便以及易于圖像處理器(graphics processing unit，GPU)加速，至今仍被大量使用．FDK算法本質(zhì)上是FBP，把得到的錐束投影數(shù)據(jù)進(jìn)行濾波，然后利用扇形束近似重建而得．FDK算法的計算公式[11]如式(1)和式(2)．

對投影數(shù)據(jù)P的加權(quán)濾波為

(1)

其中，R為光源到物體中心的距離；β為源繞中心旋轉(zhuǎn)軸z軸的旋轉(zhuǎn)角度；a為旋轉(zhuǎn)中心的橫坐標(biāo)；b為旋轉(zhuǎn)中心的縱坐標(biāo)；P(β,a,b)為投影數(shù)據(jù)；G(a)為濾波函數(shù)；符號*表示卷積運算．

反投影重建結(jié)果為

(2)

其中，U(x,y,β)=R+xcosβ+ysinβ為加權(quán)函數(shù)，這里x和y為探測器坐標(biāo)；P′(β,a,b)為加權(quán)濾波后的投影數(shù)據(jù)；其他變量定義如式(1)．

影響重建結(jié)果好壞的關(guān)鍵在于濾波函數(shù)的選擇、重建過程中插值的方式，以及錐角的大小等．從式(2)可直觀地看出，濾波函數(shù)能直接影響到反投影重建的結(jié)果．所以，合適的濾波器是實現(xiàn)高質(zhì)量重建圖像的關(guān)鍵因素之一．

2 新型濾波函數(shù)的構(gòu)建

2.1 選擇要混合和加權(quán)平均的濾波函數(shù)

要選擇合適的濾波函數(shù)進(jìn)行混合和加權(quán)平均，必須要先判斷濾波函數(shù)的優(yōu)劣特性．判斷一個濾波函數(shù)對重建結(jié)果的影響主要有：主瓣、近旁瓣和遠(yuǎn)旁瓣．主瓣高而窄，說明空間分辨率好；遠(yuǎn)旁瓣的幅度和幅值越小，說明Gibbs現(xiàn)象越小、密度分辨率越好[12-13]．范惠榮等[14]研究表明，NEW濾波函數(shù)能夠保持空間分辨率的同時減小Gibbs現(xiàn)象，而S-L濾波函數(shù)[11]可通過這3點加權(quán)平均使遠(yuǎn)旁瓣的幅度和幅值大幅減小，其近旁瓣收斂相比R-L[11]、S-L和NEW濾波函數(shù)更快，更能有效地抑制Gibbs現(xiàn)象，但因其主瓣變矮，空間分辨率會變得很差，如圖1．其中，n為采樣點；h[n]為采樣點n所對應(yīng)的濾波函數(shù)值．圖2中除R-L濾波函數(shù)外，剩下3種濾波函數(shù)的遠(yuǎn)旁瓣幾乎重疊，表明M3S-L濾波函數(shù)[10]的遠(yuǎn)旁瓣的幅度和幅值非常小，與R-L濾波函數(shù)相比無明顯振蕩，說明Gibbs現(xiàn)象非常小．因此，可將NEW濾波函數(shù)和加權(quán)平均后的M3S-L濾波函數(shù)混合疊加，得到新的NEW-MS-L濾波函數(shù)，大幅減小Gibbs現(xiàn)象，且能保持與單獨使用NEW濾波函數(shù)后相當(dāng)?shù)目臻g分辨性能．

圖1 R-L、S-L、NEW和M3S-L濾波函數(shù)的空域主瓣分布Fig.1 (Color online) Main lobe distributions in spatial domain of R-L,S-L,NEW, and M3S-L filter functions

圖2 R-L、S-L、NEW和M3S-L濾波函數(shù)的遠(yuǎn)旁瓣分布Fig.2 Far lobe distributions in spatial domain of R-L, S-L, NEW, and M3S-L filter functions

2.2 構(gòu)建新型濾波函數(shù)NEW-MS-L

新型濾波函數(shù)是根據(jù)加權(quán)平均和混合濾波的思想構(gòu)建的．首先需將S-L濾波函數(shù)進(jìn)行加權(quán)平均．

S-L濾波函數(shù)的離散形式[11]為

(3)

其中，n=0,±1,±2,…, 為采樣點(后面各式含義相同)；τ為探測器單元的大小，一般設(shè)為1．

基于文獻(xiàn)[15]的研究結(jié)果，本研究對S-L濾波函數(shù)進(jìn)行3點加權(quán)平均，記為M3S-L，該函數(shù)的離散形式[10]為

hM3S-L[n]= 0.6hS-L[n]+0.2hS-L[n-1]+

0.2hS-L[n+1]

(4)

將M3S-L濾波函數(shù)和NEW濾波函數(shù)進(jìn)行混合濾波，得到新的濾波函數(shù)NEW-MS-L．NEW濾波函數(shù)[7]和NEW-MS-L濾波函數(shù)的離散形式分別為

(5)

hNEW-MS-L[n]=k1×hM3S-L[n]+

k2×hNEW[n]

(6)

其中，k1和k2為權(quán)重，k1+k2=1． 本研究取k1=0.5；k2=0.5．

NEW-MS-L、NEW和M3SL濾波函數(shù)在空域的主瓣分布對比如圖3．

圖3 NEW、M3S-L和NEW-MS-L濾波函數(shù)的空域主瓣分布Fig.3 (Color online) Main lobe distributions in spatial domain of NEW,M3S-L, and NEW-MS-L filter functions

從圖3可見，新型濾波函數(shù)的hNEW-MS-L[n]主瓣高度與NEW濾波函數(shù)的相當(dāng)，近旁瓣幅度比NEW濾波函數(shù)的小，且收斂更快．遠(yuǎn)旁瓣的幅度和幅值很小，表明NEW-MS-L濾波函數(shù)可實現(xiàn)與NEW濾波函數(shù)相當(dāng)?shù)目臻g分辨率，且能更有效抑制Gibbs現(xiàn)象，進(jìn)一步提高密度分辨率，獲得更優(yōu)的圖像質(zhì)量．

3 仿真結(jié)果與分析

基于FDK算法，分別用NEW、R-L-S-L、R-L-N、R-L-MS-L和NEW-MS-L濾波函數(shù)對尺寸為256×256×256像素的Shepp-Logan三維頭模型[16]進(jìn)行重建，然后觀察中心面處的截面及其縱向第128行的灰度曲線圖，結(jié)果如圖4．

圖5(a)至圖5(f)分別為采用NEW、R-L-S-L、R-L-NEW、NEW-MS-L和R-L-MS-L濾波函數(shù)重建的縱向第128行的灰度曲線圖．灰度曲線圖的波動大小反應(yīng)與原始模型的差異，波動越大說明差異大，表現(xiàn)為圖像越粗糙．由圖5可見，NEW-MS-L較其他濾波函數(shù)能夠更好的還原原始模型．

圖4 原始模型圖像與用NEW-MS-L濾波函數(shù)的重建圖像Fig.4 Original model image and reconstruction image by NEW-MS-L filter function

圖5 原始模型及使用不同濾波函數(shù)重建結(jié)果的灰度曲線Fig.5 The profile pictures of original model and reconstructed results by different filter functions

僅通過灰度曲線圖的對比來說明NEW-MS-L濾波函數(shù)的優(yōu)越性是不夠的，還需通過計算NEW、R-L-S-L、R-L-N、R-L-MS-L和NEW-MS-L濾波函數(shù)的歸一化均方誤差d和歸一化平均絕對誤差r， 如式(7)和式(8)，進(jìn)行定量對比.

(7)

(8)

d反映少數(shù)點的大誤差，d=0表示重建后圖像完全再現(xiàn)了原始模型圖像；d越小表示兩者的偏差越?。畆反映多數(shù)點的小誤差，r=0表示重建圖像與頭模型原始圖像無誤差；r越小說明誤差越?。?/p>

表1為分別采用NEW、R-L-S-L、R-L-NEW、R-L-MS-L和NEW-MS-L濾波函數(shù)重建圖像后的d值和r值對比．

表1 重建結(jié)果的誤差值對比

由表1可見，采用NEW-MS-L濾波函數(shù)后的d值與采用NEW濾波函數(shù)后的d值相當(dāng)，且所得r值較其他濾波函數(shù)的都小，說明采用NEW-MS-L濾波函數(shù)所得的重建效果更佳．

上述結(jié)果還可通過調(diào)整k1和k2值來進(jìn)一步優(yōu)化[10]．d和r隨k1的變化如表2，由表2可知，當(dāng)k1=0.8時，d值最??；但當(dāng)k1=0.7、 0.8和0.9時，d值雖然相差無幾，但r值隨k1的增加而變大，且變化較明顯．所以，本研究權(quán)衡d和r的值，選取k1=0.7．

表2 誤差值d與r隨k1的變化

當(dāng)k1=0.7，k2=0.3時，可得到d=0.315 5，r=0.737 3． 雖然r值增大了，但d值進(jìn)一步減?。畬EW-MS-L濾波函數(shù)的d值和r值與表1的其他濾波函數(shù)的d值和r值進(jìn)行對比，結(jié)果反映新濾波函數(shù)性能比其他濾波函數(shù)要好．

結(jié) 語

基于混合濾波和加權(quán)平均理論，構(gòu)建新型濾波函數(shù)NEW-MS-L．其空域主瓣分布圖具有主瓣高而窄、近旁瓣小、遠(yuǎn)旁瓣衰減快等特征，能很好地抑制Gibbs現(xiàn)象，保持較高的空間分辨率．仿真結(jié)果表明，通過對比NEW、R-L-S-L、R-L-NEW、R-L-MS-L濾波函數(shù)重建結(jié)果與灰度曲線圖，并分析歸一化均方誤差和歸一化平均絕對誤差，發(fā)現(xiàn)NEW-MS-L濾波函數(shù)可在保持較高空間分辨率的情況下更有效地抑制Gibbs現(xiàn)象，密度分辨率更高，重建結(jié)果也更平滑．此外，根據(jù)實際情況適當(dāng)調(diào)整k1和k2， 能使圖像效果達(dá)到符合預(yù)期的要求．

NEW-MS-L濾波函數(shù)在本次仿真實驗中主要針對傳統(tǒng)基于吸收效應(yīng)的CT，今后可進(jìn)一步探討NEW-MS-L濾波函數(shù)在相襯CT[17-22]中的應(yīng)用．

/ References：

[1] Egan C K,Jacques S D M,Wilson M D, et al.3D chemical imaging in the laboratory by hyperspectral X-ray computed tomography[J].Scientific Reports,2015,5: 15979.

[2] Bieberle M,Barthel F.Combined phase distribution and particle velocity measurement in spout fluidized beds by ultrafast X-ray computed tomography[J].Chemical Engineering Journal,2016,285:218-227.

[3] Edlund R,Skorpil M,Lapidus G,et al.Cone-beam CT in diagnosis of scaphoid fractures[J].Skeletal Radiology,2016,45(2):197-204.

[4] Senyshyn A,Mühlbauer M J,Dolotko O,et al.Homogeneity of lithium distribution in cylinder-type Li-ion batteries[J].Scientific Reports,2015,5:18380.

[5] Plessis A D, Olawuyi B J, Boshoff W P, et al. Simple and fast porosity analysis of concrete using X-ray computed tomography[J]. Materials and Structures, 2016, 49(1): 553-562.

[6] Feldkamp L A, Davis L C, Kress J W. Practical cone-beam algorithm[J]. Journal of the Optical Society of America A, 1984, 1(6): 612-619.

[7] Wei Yuchuan, Wang Ge, Hsieh J. An intuitive discussion on the ideal ramp filter in computed tomography (I)[J]. Computers & Mathematics with Applications, 2005, 49(5/6): 731-740.

[8] 謝 強(qiáng)．計算機(jī)斷層成像技術(shù):原理、設(shè)計、偽像和進(jìn)展[M]．張朝宗，郭志平，王宗賢，譯．北京：科學(xué)出版社，2006． Hsieh J.Computed tomography: principle, design, artifacts and recent advances[M]. Zhang Chaozong,Guo Zhiping, Wang Zongxian, trans. Beijing: Science Press,2006.(in Chinese)

[9] 胡永梅，全 紅，徐利明，等．新型混合濾波器在FDK算法中的應(yīng)用[J]．中國醫(yī)學(xué)影像技術(shù)，2011，27(2)：397-400． Hu Yongmei, Quan Hong, Xu Liming, et al. Application of a new mixed filter in FDK algorithm [J]. Chinese Journal of Medical Imaging Technology, 2011, 27(2): 397-400.(in Chinese)

[10] 王曉鵬，王明泉，侯慧玲．基于R-L-MS-L濾波函數(shù)的CT圖像重建[J]．電視技術(shù)，2014，38(7)：26-28． Wang Xiaopeng,Wang Mingquan,Hou Huilin. CT image reconstruction based on R-L-MS-L filter function[J].Video Engineering, 2014,38(7): 26-28.(in Chinese)

[11] 閆 鑌，李 磊．CT圖像重建算法[M]．北京：科學(xué)出版社，2014． Yan Bin,Li Lei.CT image reconstructinalgorithm[M].Beijing:Science Press,2014.(in Chinese)

[12] 鄭君里．信號與系統(tǒng)[M]．北京：高等教育出版社，2002． Zhen Junli. Signals and systems[M].Beijing:Higher Education Press,2002.(in Chinese)

[13] 楊位欽，顧 嵐．時間序列分析與動態(tài)數(shù)據(jù)建模[M]．北京：北京工業(yè)學(xué)院出版社，1986． Yang Weiqin,Gu Lan. Time series analysis and dynamic data modeling [M].Beijing:Beijing Institute of Tchnology Press,1986.(in Chinese)

[14] 范惠榮，徐茂林，邱 鈞，等．關(guān)于CBP算法的一種新型濾波函數(shù)和它的性質(zhì)[J]．電子學(xué)報，2004，32(2)：232-235． Fan Huirong,Xu Maolin,Qiu Jun,et al. A new filter function in CBP and its properties[J].Acta Electronica Sinica,2004,32(2):232-235.(in Chinese)

[15] 劉 曉，楊朝文．RL濾波函數(shù)的改進(jìn)對卷積反投影圖像重建的影響[J]．四川大學(xué)學(xué)報自然科學(xué)版， 2004，41(1)：112-117． Liu Xiao,Yang Chaowen. The influence of modified-RL filter function on the quality of image reconstruction using convolution backprojection algorithm [J].Journal of Sichuan University Natural Science Edition,2004,41(1):112-117.(in Chinese)

[16] 劉 澤，孫豐榮，李艷玲，等．基于3D Shepp-Logan頭部模型的三維醫(yī)學(xué)圖像重建仿真[J]．生物醫(yī)學(xué)工程學(xué)雜志，2006，23(5)：938-943． Liu Ze,Sun Fengrong,Li Yanling, et al.The three-dimension medical image reconstruction simulation on 3D Shepp-Logan head phantom [J].Journal of Biomedical Engineering,2006,23(5):938-943.(in Chinese)

[17] Lei Yaohu,Huang Jianheng,Liu Xin,et al.Improvement of visibility of moir fringe in X-ray differential phase-contrast imaging[J].Journal of Shenzhen University Science and Engineering,2016, 33(5):506-510.

[18] 劉 鑫，郭金川．部分相干光源微分相襯雙圖像相位恢復(fù)[J]．深圳大學(xué)學(xué)報理工版，2014，31(2)：169-173． Liu Xin,Guo Jinchuan.Dual-images phase retrieval in differential phase-contrast imaging with partial coherence source[J].Journal of Shenzhen University Science and Engineering,2014,31(2):169-173.(in Chinese)

[19] Jerjen I, Revol V, Kottler C, et al. Phase contrast cone beam tomography with an X-ray grating interferometer [C]// International Conference on Advanced Phase Measurement Methods in Optics & Imaging.Ascona, Switzerland: American Institute of Physics, 2010, 1236(1):227-231.

[20] Huang Zhifeng, Kang Kejun, Li Zheng, et al. Direct computed tomographic reconstruction for directional-derivative projections of computed tomography of diffraction enhanced imaging[J]. Applied Physics Letters, 2006, 89(4): 041124.

[21] McCann M T, Nilchian M, Stampanoni M, et al. Fast 3D reconstruction method for differential phase contrast X-ray CT[J]. Optics Express, 2016, 24(13): 14564-14581.

[22] Fu Jian, Hu Xinhua, Velroyen A, et al. 3D algebraic iterative reconstruction for cone-beam X-ray differential phase-contrast computed tomography[J]. PloS One, 2015, 10(3): e0117502.

【中文責(zé)編：英 子；英文責(zé)編：木 南】

2016-11-29；Accepted：2017-03-09

Professor Guo Jinchuan. E-mail: jcguo@szu.edu.cn

A new filter function for the image reconstruction of cone beam X-ray CT

Hong Zhenhou, Zhou Bin, and Guo Jinchuan

College of Optoelectronic Engineering, Key Laboratory of Optoelectronic Devices and Systems of Ministry of Education and Guangdong Province, Shenzhen University, Shenzhen 518060, Guangdong Province, P.R.China

The conventional filter function of the FDK algorithm in X-ray computed tomography (CT) can result in the obvious Gibbs phenomenon, which has great effect on the CT image. A NEW-MS-L filter function is deduced based on the rationale of mixed filter theory and weighted average theory. First, the S-L filter function is weighted and averaged to form M3S-L filter function.Then, as a superposition of the M3S-L filter function and the NEW filter function, the NEW-MS-L filter function is obtained. Compared with the other filter functions such as NEW, R-L-S-L, R-L-NEW, and R-L-MS-L, the NEW-MS-L can suppress the Gibbs phenomenon while maintaining high resolution. The results show that the NEW-MS-L filter function can provide X-ray CT with better reconstruction images.

optical engineering; computed tomography; FDK algorithm; filter function; Gibbs phenomenon; mixed filter theorem; weighted average theorem

：Hong Zhenhou,Zhou Bin,Guo Jinchuan.A new filter function for the image reconstruction of cone beam X-ray CT [J]. Journal of Shenzhen University Science and Engineering, 2017, 34(3): 284-289.(in Chinese)

國家自然科學(xué)基金資助項目(11674232，11074172)

洪振厚(1991—)，男，深圳大學(xué)碩士研究生．研究方向：X光成像．E-mail：paul_hongzh@163.com

TP 391.41

A

10.3724/SP.J.1249.2017.03284

Foundation：National Natural Science Foundation of China(11674232, 11074172)

引 文：洪振厚，周 彬，郭金川．錐束X射線CT圖像重建的新型濾波函數(shù)[J]. 深圳大學(xué)學(xué)報理工版，2017，34(3)：284-289.