基于極坐標的扇束CT周期性優(yōu)化算法

劉　輝，張　權(quán)，劉　祎，桂志國,2

(1.中北大學 電子測試技術(shù)重點實驗室，山西 太原 030051；2.中北大學 儀器科學與動態(tài)測試教育部重點實驗室，山西 太原 030051)

?

基于極坐標的扇束CT周期性優(yōu)化算法

劉輝1，張權(quán)1，劉祎1，桂志國1,2

(1.中北大學 電子測試技術(shù)重點實驗室，山西 太原 030051；2.中北大學 儀器科學與動態(tài)測試教育部重點實驗室，山西 太原 030051)

為了提高傳統(tǒng)扇束計算機斷層掃描(Computed Tomography,CT)卷積反投影重建算法(Conventional Back-Projection,CBP )的重建速度，提出了一種快速的扇束CT重建算法.該算法結(jié)合極坐標特性，利用三角函數(shù)一定程度上的周期性，對傳統(tǒng)CBP重建算法進行了改進，減少了正余弦函數(shù)的計算量和重建所需要的外圍循環(huán)次數(shù).同樣利用該特性對極坐標重建后的圖像坐標轉(zhuǎn)換過程進行改善，從而對CBP算法進行了改進，降低了重建時間.

計算機斷層掃描；卷積反投影；三角函數(shù)特性；極坐標

隨著計算機斷層成像技術(shù)(Computed Tomography，CT)的不斷發(fā)展，CT重建技術(shù)越來越多地運用到醫(yī)學診斷[1]、工業(yè)無損檢測等行業(yè)中[2].CT重建主要有解析算法和迭代算法[3-5].與迭代算法相比，解析重建法以其重建速度快、重建圖像質(zhì)量高的特點，已經(jīng)被廣泛應用到工程上.CBP重建算法是解析算法中最具代表性的重建算法之一[6-8].雖然CBP重建算法重建所需數(shù)據(jù)量較少，但是運算數(shù)據(jù)量依然很大，還需要進一步提高其重建速度，才能更好地運用到工業(yè)上.

加權(quán)、卷積濾波和反投影是CBP算法的3個主要部分，其中反投影過程就占據(jù)了98.36%的重建時間[9]，因此反投影重建過程的改進是主要研究方向.文獻[10]運用極坐標特性，提出了極坐標反投影重建算法(Polar Coordinates Back-Projection,PCBP)，但極坐標向笛卡爾坐標的轉(zhuǎn)換過程中要使用插值運算，重建速度提高有限；文獻[11]利用正余弦函數(shù)特性，將16幅投影同時進行反投影重建，提出了對極坐標重建的優(yōu)化算法(Optimize Polar Conventional Back-Projection OPCBP)，該算法將重建速度提高了4倍；文獻[12]利用了三角函數(shù)的對稱性，提出了極坐標反投影對稱優(yōu)化算法(Symmetry Polar Conventional Back-Projection SPCBP)，同時對64幅預處理后的數(shù)據(jù)進行重建，將重建速度提高8倍.

本文利用三角函數(shù)一定程度上的周期性和極坐標反投影重建的特點，將144幅反投影數(shù)據(jù)同時進行反投影重建，使外圍循環(huán)減少至傳統(tǒng)算法的1/144，同時大大減少了三角函數(shù)的運算量，從而實現(xiàn)對傳統(tǒng)算法的優(yōu)化.本文稱之為極坐標反投影周期性優(yōu)化算法(Periodicity Polar-coordinate Conventional Back-Projection，PPCBP).

1　扇束CBP重建原理

扇束CT掃描示意圖如圖1所示，圖中a(r,θ)是待重建圖像切片的射線衰減系數(shù)分布函數(shù)；β為射線源投射角度；函數(shù)p(β,s)是該角度下的投影數(shù)據(jù)；D為射線源到物體重建中心的距離，重建公式為

(1)

(2)

(3)

(4)

(5)

由式(1)～式(5)可知計算機實現(xiàn)CBP反投影重建時，首先通過式(3)給投影數(shù)據(jù)進行加權(quán)；其次通過式(2)對加權(quán)過的數(shù)據(jù)進行卷積濾波；最后根據(jù)式(4)和式(5)確定的射線的反投影地址s和加權(quán)因子U對濾波后的數(shù)據(jù)進行累加反投影重建.假設重建的圖像的大小為N×N，投影角度為[0°,360°),投影視圖數(shù)為M，該算法的反投影過程流程如圖2所示，從圖2中可以看出CBP反投影重建中運最復雜的就是對反投影地址s和加權(quán)因子U求解.由式(4)和式(5)可知，每次求s和U需要9次乘法(含除法)和6次正余弦函數(shù)運算，反投影過程總共需要次乘法和次正余弦運算，因此主要考慮對s和U求解過程的簡化.

圖1　扇束掃描示意圖Fig.1　The scanning schematic of fan-beam

圖2　CBP反投影重建過程流程圖Fig.2　The flowchart of CBP back-projection reconstruction

2　極坐標下進行重建

(6)

(7)

圖3　PCBP反投影重建流程圖Fig.3　The flowchart of PCBP algorithm

3　三角函數(shù)周期性優(yōu)化

3.1反投影過程的優(yōu)化

三角函數(shù)中30°，45°，60°和90°是幾個特殊的角度，三角函數(shù)值很容易計算得出.通過簡單的三角函數(shù)的誘導公式，可以得出120°,135°,150°,180°,210°,225°,240°,270°,300°,315°和330°的三角函數(shù)值，這些角度中除了45°及其誘導產(chǎn)生的角度之外，其余角度都是呈現(xiàn)出的等差關(guān)系，這些角度的正余弦值一定程度上表現(xiàn)出了周期性，給PCBP算法優(yōu)化提供了很好的條件.分別讓β-θ=30m(0≤m≤11,m∈Z)，討論其正余弦值并研究各角度下的余弦值之間的關(guān)系

(8)

(9)

(10)

(11)

圖4　PPCBP反投影重建流程圖Fig.4　The flowchart of PPCBP algorithm

式(8)～式(10)之間關(guān)系只是差了一個正負號，而且A,B和A1,B1以及A2,B2具有式(12)所示的關(guān)系，因此只要求解出[0°,30°)范圍內(nèi)的正余弦值就可以求得所有掃描角度下β-θ的正余弦值，極坐標三角函數(shù)周期性優(yōu)化算法(PPCBP)算法的反投影過程流程如圖4所示.

(12)

通過以上分析，在一個投影角度β下通過三角函數(shù)的變換，只需要計算兩次正余弦值就可以計算出θ+30m(0≤m≤11,m∈z)下的s和U的值，因此一個循環(huán)可以完成12個投影視角下的144個投影數(shù)據(jù)參與重建，外層循環(huán)減少為極坐標下的1/144，同時將三角函數(shù)計算量降低為極坐標下的1/72，大大提高了重建速度.所以PPCBP算法無論在循環(huán)，乘法或者是三角函數(shù)運算上都得到了很大的提高，在理論上證實了PPCBP的優(yōu)越性.

3.2重建圖像坐標轉(zhuǎn)換的優(yōu)化

本文在圖像的轉(zhuǎn)換過程中采用雙線性插值方法，如圖5所示，雙線性插值通過臨近目標點G的A,B,C和D 4點經(jīng)過3次線性插值從而求得G點像素值的過程，其原理如式(13)所示，其中P(x)代表了x點的像素值，a=BE/AB,b=DF/CD，c=FG/EF.

(13)

針對采用雙線性插值增加的計算量，本文使用正切函數(shù)對稱性對其簡化.圖6中所標記的點在極坐標中具有相同的值，且其正切值之間具有正負和互為倒數(shù)的關(guān)系.因此，如果知道極坐標下任意一個點的坐標，就能同時知道極坐標下其他7個點的坐標轉(zhuǎn)換，使外圍循環(huán)減少至改變前的1/8，同時使正余切函數(shù)值的計算量也減少至改變前的1/8，提高了坐標轉(zhuǎn)換的速度.

圖5　雙線性插值示意圖Fig.5　The schematic of bilinear interpolation

圖6　坐標轉(zhuǎn)換中像素點的選取示意圖Fig.6　Pixels chosen schematic during coordinate transformation

4　實驗結(jié)果

圖7　各種算法相對CBP重建時間加速比Fig.7　The speedup time compared with CBP

為了驗證本文算法的可行性，將PTCBP算法和CBP算法、PCBP算法、OPCBP算法和SPCBP算法的重建結(jié)果進行對比.實現(xiàn)本文算法的軟件設備為Windows 7 32 b操作系統(tǒng)和VC++6.0編譯環(huán)境；硬件設備為英特爾 Pentium(奔騰)雙核 E5400@2.70 GHz處理器，DDR3 2G金士頓內(nèi)存.投影數(shù)據(jù)參數(shù)如下：射線源到旋轉(zhuǎn)中心的距離為800 mm，物體中心到探測器的距離為150 mm，射線為ISOVOLT 450，探測器為paxscan2520，探元大小0.127 mm，電壓110 kV，電流2.5 mA，投影數(shù)據(jù)大小為570×460，投影視圖數(shù)為360，重建圖像大小為512×512.

4.1重建速度比較

各算法的重建時間如表1所示，加速比如圖7所示：從重建的時間和加速比上可以看出，PTCBP算法的速度是CBP算法的將近12.5倍，是PCBP算法的近10倍，是SPCBP的近1.5倍，是OPCBP的2.1倍.因此，本文所提出的PPCBP算法大大提高了CBP算法的重建速度.由于使用的是三角函數(shù)一定程度上表現(xiàn)出來的周期性原理，使得重建過程中所占用的內(nèi)存也大大減少，在計算機編程中充分展示了PPCBP算法的優(yōu)越性.

表1　各算法重建速度對比表

4.2重建圖像質(zhì)量比較

各種算法圖像如圖8所示，從圖中我們很難看出圖像之間的差別，圖9展示的是各算法重建出圖像的第145列像素大小之間的對比，從中我們發(fā)現(xiàn)各圖像之間差別很小，說明了利用PPCBP算法重建出的圖像與其他重建算法重建出的圖像在質(zhì)量上沒有明顯的差別.為了量化本文算法重建出的圖像與其他算法重建圖像之間的差別，本文采用歸一化均方距離 (Normalized Mean Square Distance,NMSD)和歸一化平均絕對距離(Normal Average Absolute Distance,NAAD)對圖像質(zhì)量進行評估，其公式分別為

(14)

(15)

圖8　各種算法重建圖像質(zhì)量對比圖Fig.8　The quality comparison chart of the images constructed by different algorithms

圖9　各算法第145列像素大小對比圖Fig.9　The comparison of pixels in the 145th row but from different images reconstructed by vary algorithms.

算法CBPPCBPPTCBPOPTCBPSPCBPNMSD-0.3881030.3881080.3881040.405673NAAD-0.1984530.1984590.1984540.217913

由表2可知，無論是從歸一化均方距離還是歸一化平均絕對距離上PTCBP算法跟PCBP,OPTCBP和SPCBP算法重建出的圖像在質(zhì)量上沒有明顯的差別，因此進一步的說明PTCBP算法在重建圖像時并沒有引入誤差，PTCBP算法在保證圖像質(zhì)量上有良好的適用性.

5　結(jié)　論

本文提出的PTCBP算法結(jié)合三角函數(shù)一定程度上的周期性以及對稱性分別對反投影過程和坐標轉(zhuǎn)換過程進行了改進.在保證了重建圖像質(zhì)量的前提下，很大程度上提高了傳統(tǒng)CBP算法的重建速度.

[1]曾更生.醫(yī)學圖像重建[M].北京：高等教育出版社,2010.

[2]Franco L,Tahoces P G,Martínez-Mera J A.Visualization software for CT：fan/cone beam and metrology applications [J].Procedia Engineering,2013,63(63)：779-785.

[3]莊天戈.CT原理與算法[M].上海：上海交通大學出版社，1992.

[4]Rit S,Oliva M V,Brousmiche S,et al.The Reconstruction Toolkit (RTK),an open-source cone-beam CT reconstruction toolkit based on the Insight Toolkit (ITK)[C].Journal of Physics：Conference Series.2014,489：012079.

[5]Wang A S,Stayman J W,Otake Y,et al.Soft-tissue imaging with C-arm cone-beam CT using statistical reconstruction[J].Physics in Medicine and Biology.2014,59(4)：1005.

[6]Zhang F,Du J,Jiang H,et al.Iterative geometric calibration in circular cone-beam computed tomography[J].Optik-International Journal for Light and Electron Optics,2014,125(11)：2509-2514.

[7]Zheng Fang,Xin Gu,Xing Liang Zhang,et al.Research on the relationship between projection number and image noise level in FBP reconstruction [J].Advanced Materials Research.2013,718-720：2324-2328.

[8]Xin Y,Yang J,Xie Z Q,et al.An overlapping semantic community detection algorithm base on the ARTs multiple sampling models[J].Expert Systems with Applications,2014,42(7)：3420-3432.

[9]劉曉平.扇束卷積反投影法的程序優(yōu)化[J].CT 理論與應用研究.1996,5(1)：35-37.

Liu Xiaoping.Program optimization of swctor-beam convolution backprojection[J].Computerized Tomography Theory and Applications,1996,5(1)：35-37.(in Chinese)

[10]楊民,路宏年,黃朝志.用查找表和極坐標反投影法實現(xiàn)計算機斷層掃描快速重建[J].兵工學報，2004,35(4)：476-479.

Yang Min,Lu Hongnian,Huang Chaozhi.Fast computed tomography reconstruction with look-up table and polar coordonate back-projection[J].Acta Armamentarii,2004,35(4)：476-479.(in Chinese)

[11]張鑾，桂志國.扇束等距CT極坐標反投影重建算法的速度優(yōu)化[J].無損檢測，2010,32(2)：95-98.

Zhang Luan,Gui Zhiguo.Speed optimization of polar coordinate back-projection rectonstruction algorithm for fan beam collinear equispaced CT[J].Nondestructive Testing,2010,32(2)：95-98.(in Chinese)

[12]張晶，張權(quán)，劉祎，等.扇束CT極坐標反投影重建算法的對稱優(yōu)化[J].計算應用.2014,34(6)：1711-1714.

Zhang Jing,Zhang Quan,Liu Yi,et al.Symmetry optimization of polar coordinate back-projection reconstruction algorithm for fan beam CT[J].Journal of Computer Applications,2014,34(6)：1711-1714.(in Chinese)

Periodic Optimization Algorithm Based on Fan-Beam CT Polar Coordinates

LIU Hui1,ZHANG Quan1,LIU Yi1,GUI Zhiguo1,2

(1.State Key Laboratory of Electronic Measurement Technology (North University),Taiyuan 030051,China；2.Instrument Science & Dynamic Measurement,Ministry of Education Key laboratory of (North University),Taiyuan 030051,China.)

To improve the reconstruction speed of traditional fan-beam CT image reconstruction based on fan-beam conventional back projection (CBP),a fast polar coordinate back projection reconstruction algorithm was proposed.Combining with the features of polar coordinate,the proposed algorithm applied the trigonometric periodicity in some degree,and improved traditional CBP to reduce the amount of computation of trigonometric function and the times of external cycle in reconstruction process.Similarly,the mentioned properties were also used in the process of the conversion between polar coordinate and Cartesian coordinate to improve the CBP algorithm and the reconstruction time gets reduced.

computed tomography; conventional back-projection; trigonometric function characteristics; polar coordinate

1671-7449(2016)05-0394-06

2016-03-01

國家自然科學基金資助項目(61271357)；國家重大科學儀器設備開發(fā)專項資助項目(2014YQ240445)；山西省自然科學基金資助項目(2015011046)

劉輝(1989-)，男，碩士生，主要從事圖像重建的研究.

TP391.4

Adoi：10.3969/j.issn.1671-7449.2016.05.005