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

(1.中北大學(xué) 電子測(cè)試技術(shù)重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051；2.中北大學(xué) 儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051)

?

基于極坐標(biāo)的扇束CT周期性?xún)?yōu)化算法

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

(1.中北大學(xué) 電子測(cè)試技術(shù)重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051；2.中北大學(xué) 儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室，山西 太原 030051)

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

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

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

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

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

1　扇束CBP重建原理

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

(1)

(2)

(3)

(4)

(5)

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

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

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

2　極坐標(biāo)下進(jìn)行重建

(6)

(7)

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

3　三角函數(shù)周期性?xún)?yōu)化

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

三角函數(shù)中30°，45°，60°和90°是幾個(gè)特殊的角度，三角函數(shù)值很容易計(jì)算得出.通過(guò)簡(jiǎn)單的三角函數(shù)的誘導(dǎo)公式，可以得出120°,135°,150°,180°,210°,225°,240°,270°,300°,315°和330°的三角函數(shù)值，這些角度中除了45°及其誘導(dǎo)產(chǎn)生的角度之外，其余角度都是呈現(xiàn)出的等差關(guān)系，這些角度的正余弦值一定程度上表現(xiàn)出了周期性，給PCBP算法優(yōu)化提供了很好的條件.分別讓?duì)?θ=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)系只是差了一個(gè)正負(fù)號(hào)，而且A,B和A1,B1以及A2,B2具有式(12)所示的關(guān)系，因此只要求解出[0°,30°)范圍內(nèi)的正余弦值就可以求得所有掃描角度下β-θ的正余弦值，極坐標(biāo)三角函數(shù)周期性?xún)?yōu)化算法(PPCBP)算法的反投影過(guò)程流程如圖4所示.

(12)

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

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

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

(13)

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

圖5　雙線(xiàn)性插值示意圖Fig.5　The schematic of bilinear interpolation

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

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

圖7　各種算法相對(duì)CBP重建時(shí)間加速比Fig.7　The speedup time compared with CBP

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

4.1重建速度比較

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

表1　各算法重建速度對(duì)比表

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

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

(14)

(15)

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

圖9　各算法第145列像素大小對(duì)比圖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可知，無(wú)論是從歸一化均方距離還是歸一化平均絕對(duì)距離上PTCBP算法跟PCBP,OPTCBP和SPCBP算法重建出的圖像在質(zhì)量上沒(méi)有明顯的差別，因此進(jìn)一步的說(shuō)明PTCBP算法在重建圖像時(shí)并沒(méi)有引入誤差，PTCBP算法在保證圖像質(zhì)量上有良好的適用性.

5　結(jié)　論

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

[1]曾更生.醫(yī)學(xué)圖像重建[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].上海：上海交通大學(xué)出版社，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 理論與應(yīng)用研究.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]楊民,路宏年,黃朝志.用查找表和極坐標(biāo)反投影法實(shí)現(xiàn)計(jì)算機(jī)斷層掃描快速重建[J].兵工學(xué)報(bào)，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]張鑾，桂志國(guó).扇束等距CT極坐標(biāo)反投影重建算法的速度優(yōu)化[J].無(wú)損檢測(cè)，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極坐標(biāo)反投影重建算法的對(duì)稱(chēng)優(yōu)化[J].計(jì)算應(yīng)用.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

國(guó)家自然科學(xué)基金資助項(xiàng)目(61271357)；國(guó)家重大科學(xué)儀器設(shè)備開(kāi)發(fā)專(zhuān)項(xiàng)資助項(xiàng)目(2014YQ240445)；山西省自然科學(xué)基金資助項(xiàng)目(2015011046)

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

TP391.4

Adoi：10.3969/j.issn.1671-7449.2016.05.005