洪賢勇，喬志偉

(中北大學(xué)計(jì)算機(jī)與控制工程學(xué)院，山西太原030051)

計(jì)算機(jī)斷層成像技術(shù)(Computed Tomography，CT)是現(xiàn)今較好的無損檢測技術(shù)之一，已廣泛應(yīng)用于人體組織成像和工業(yè)無損檢測中［1－2］。CT算法主要包括解析法和迭代法［3］兩種，由于解析法具有重建時(shí)間快的特點(diǎn)而處于主導(dǎo)地位，解析法中具有代表性的是濾波反投影(FBP)算法和反投影濾波(BPF)算法。

FBP算法出現(xiàn)于20世紀(jì)60年代，并很快從平行束掃描模式推廣到圓軌道錐束掃描模式。然而，1984年提出的FDK算法只是一種近似算法，此后近二十年，研究者們一直在尋找精確解析三維FBP算法。2002年到2004年，Katsevich提出了螺旋錐束FBP算法以及一些改進(jìn)算法［4－5］，這一算法具有公式簡單、易于實(shí)現(xiàn)等優(yōu)點(diǎn)，是精確三維錐束CT的一個(gè)重大突破，但是這是一種全局重建算法，在處理小物體時(shí)具有優(yōu)勢，當(dāng)處理比較大的物體時(shí)，并且只需要重建物體局部的圖像時(shí)，這種算法就出現(xiàn)了局限性。2004年，潘曉川課題組提出了一種基于PI線的BPF算法［6－7］，由于這一算法是通過求導(dǎo)、反投影、濾波(有限希爾伯特變換)來重建圖像，因而在處理數(shù)據(jù)截?cái)鄷r(shí)具有很好的效果［8－9］，而且該算法能夠?qū)崿F(xiàn)局部重建［10－11］。

為了實(shí)現(xiàn)精確的局部重建，本文在三維螺旋錐束BPF算法的基礎(chǔ)上，設(shè)計(jì)了一種利用微分和有限希爾伯特變換的局部特性來重建局部圖像的平行束BPF算法，這一算法能夠?qū)崿F(xiàn)精確的局部重建，而且有效解決了重建過程中出現(xiàn)的亮條偽像問題［12］，在工程中具有很好的實(shí)用性。

1 BPF重建算法的原理

BPF算法是由芝加哥大學(xué)潘曉川研究組提出的，其最初形式是針對螺旋錐束三維CT掃描幾何的。

后來，該組將這種方法推廣到了任意掃描幾何軌跡;又將該算法改寫成了扇形束二維CT的形式。

本節(jié)將分析該算法的原理、關(guān)鍵技術(shù)、實(shí)現(xiàn)公式及有限希爾伯特變換的實(shí)現(xiàn)步驟。

1.1 算法實(shí)現(xiàn)的框架圖

根據(jù)螺旋錐束三維CT的基于PI線的BPF算法，將重建的基本步驟和思想移植到平行束二維CT，可以得到圖1所示的重建步驟。

圖1 BPF算法的流程圖

但是該步驟只是一種重建思想的移植，所以1.2節(jié)，給出了算法的推導(dǎo)。

1.2 平行束BPF算法

1)斜變?yōu)V波可以分解為微分和希爾伯特變換［13］。

傅里葉變換的兩個(gè)性質(zhì):

(1)在傅里葉域中乘以j2πω相當(dāng)于在空域(s域)中求導(dǎo)。

由以上兩個(gè)性質(zhì)可知，經(jīng)典的FBP算法可以分解成由微分、希爾伯特變換和反投影三部分組成，改變它們之間的順序可以實(shí)現(xiàn)不同的算法，本文算法來源于上面的基本思想，但是具體的實(shí)現(xiàn)公式是基于潘曉川的BPF重建公式并假定平行束掃描幾何是螺旋錐束掃描幾何的特殊形式而改寫得到的。

2)本文算法的實(shí)現(xiàn)公式［14］

(1)對每個(gè)角度下得到的平行投影數(shù)據(jù)取ROI區(qū)域投影并求導(dǎo)，得到局部微分投影。

式中:p(t，θ)是在角度θ下的投影;p'(t，θ)是投影信號(hào)的微分。

(2)對這一局部微分投影反投影到ROI區(qū)域，反投影對應(yīng)的角度范圍為0～π。

式中:t=xcosθ+ysinθ;(x，yπ)是平面上的一點(diǎn)，x是該點(diǎn)在x軸上的坐標(biāo)，yπ是該點(diǎn)在y軸上的坐標(biāo)，也是PI坐標(biāo)軸的坐標(biāo)。注意，在這里PI線是平行于y軸的。

(3)對反投影圖像沿著PI線(如圖2)進(jìn)行有限Hilbert變換，得到物體的ROI重建圖。

圖2 局部平行PI示意圖

式中:l，s分別為PI線的兩個(gè)端點(diǎn)，如圖2所示;C(0)為沿著PI線方向的投影，也就是當(dāng)旋轉(zhuǎn)角度為0時(shí)的投影，物體的平行旋轉(zhuǎn)投影圖見圖3。

圖3 實(shí)驗(yàn)用仿真模型投影圖

1.3 有限希爾伯特變換的實(shí)現(xiàn)步驟

在整個(gè)算法中，有限希爾伯特變換的實(shí)現(xiàn)方式是最關(guān)鍵的一個(gè)步驟。如果處理不好，將使得圖像的兩邊有亮條狀偽像。

本文在實(shí)現(xiàn)該變換的時(shí)候，采用了加權(quán)希爾伯特變換的方式，完全避免了這種偽像。

具體的實(shí)現(xiàn)步驟如下:

1)第一步，對PI線上的反投影圖像做加權(quán)，得到

2)第二步，對每個(gè)PI線上的反投影加權(quán)圖像，做一個(gè)無限希爾伯特變換，并取值［l，s］;

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

2.1 實(shí)驗(yàn)仿真模型說明

在1個(gè)由7個(gè)小圓柱內(nèi)插到1個(gè)大圓柱的某斷層中利用反投影濾波算法進(jìn)行重建，實(shí)驗(yàn)所用仿真模型見圖3。探測器在每個(gè)角度分別采集201個(gè)點(diǎn)(采樣間隔為0.05 cm)，共采集180個(gè)角度的投影數(shù)據(jù)(角度間隔為1°)，采用像素驅(qū)動(dòng)方式進(jìn)行反投影，卷積運(yùn)算的方式實(shí)現(xiàn)濾波，重建出的圖像為以旋轉(zhuǎn)中心為中心、大小為201×201的整體重建圖。實(shí)驗(yàn)?zāi)Ｐ偷膮?shù)見表1。

表1 實(shí)驗(yàn)?zāi)Ｐ偷膮?shù)說明

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

首先用FBP算法和BPF算法做完全重建(整體投影正弦圖見圖4);然后以右邊半圓作為ROI區(qū)域，取僅僅覆蓋ROI區(qū)域的局部投影(見圖5)，用FBP和BPF兩種算法重建。

圖4 投影正弦圖

圖5 右半圓投影正弦圖

2.3 結(jié)果分析

由于FBP算法和BPF算法重建的整體圖都是精確重建，因而引入歸一化均方距離判距d和歸一化平均絕對距離判距r來分析FBP算法和BPF算法的精度［15］(比較結(jié)果見表2)

式中:tu，v，ru，v分別表示原始模型和重建后的圖像的第u行、v列的像素密度;為測試模型密度的值;圖像的像素為N×N。由式子分析可知d越大標(biāo)明與原始圖像的區(qū)別越大，也就是說圖像越不準(zhǔn)確，而r也類似的說明與原始圖像的誤差，不過它側(cè)重反映出現(xiàn)多數(shù)點(diǎn)存在小誤差的情況。

表2 本文全局算法與經(jīng)典FBP全局算法誤差比較

分析圖6和圖7，以及表2中的數(shù)據(jù)，可以看出本文所對應(yīng)的全局重建算法可以達(dá)到精確重建的目的，而且其精度與FBP的精度相當(dāng)。

圖6 BPF算法整體圖

圖7 FBP算法整體圖

分析圖8和圖9，兩者所取的局部微分投影一致，但是所使用的算法不同，前者使用本文提出的BPF局部重建算法，后者則是經(jīng)典的FBP重建算法。圖9重建的圖像出現(xiàn)了數(shù)據(jù)的丟失，而圖8則顯示了精確的局部重建圖像，這是由于本文算法具有局部重建的特性，而FBP算法不具備這一特性，由此可以看出本文算法能夠?qū)崿F(xiàn)精確的局部重建。

圖8 BPF右半圓圖

圖9 FBP右半圓圖

另外要特別指出的是，圖8沒有出現(xiàn)先前文獻(xiàn)中所出現(xiàn)的頂部和底部的亮條偽像，其根本原因是采用了1.3節(jié)所描述的有限希爾伯特變換實(shí)現(xiàn)算法。

3 小結(jié)

本文的算法是根據(jù)潘曉川BPF算法的基礎(chǔ)上設(shè)計(jì)的，通過對ROI區(qū)域投影數(shù)據(jù)求導(dǎo)，然后反投影到ROI區(qū)域，最后沿著PI線進(jìn)行有限希爾伯特濾波得到局部重建圖像，仿真實(shí)驗(yàn)的結(jié)果說明，本文算法所對應(yīng)的全局重建是一種精度很高的重建算法;當(dāng)所取的局部投影區(qū)域和重建區(qū)域一致時(shí)，可以精確重建出局部圖像，也就是在投影覆蓋ROI而且僅僅覆蓋ROI時(shí)能夠?qū)崿F(xiàn)對這一局部區(qū)域的精確重建，說明本文的算法是一種精確的局部重建算法。

當(dāng)遇到大物體的時(shí)候，一般只需要分析物體某一部分或者檢查某一部分是否出現(xiàn)問題，此時(shí)采用這種方法，一方面可以減少射線劑量，另一方面可以快速重建?？梢夿PF算法是一種快速有效的ROI重建方法。

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

［2］謝強(qiáng).計(jì)算機(jī)斷層成像技術(shù)［M］.北京:科學(xué)出版社，2006.

［3］ JIANGM，WANG G.Convergence studies on iterative algorithms for image reconstruction［J］.IEEE Trans.Medical Imaging，2003，22(5):569－579.

［4］ KATSEVICH A.A general scheme for constructing inversion algorithm for cone beam CT［J］.International Journal of Mathematics and Mathematical Science，2003(21):1305－1321.

［5］ YU H，WANG G.Studies on artifacts of the Katsevich algorithm of spiral cone－beam CT［J］.Proceeding of SPIE，2004(35):540－549.

［6］ ZOUY，PANX.Exact image reconstruction on PI－lines fromminimum data in helical cone beam CT［J］.Phys.Med.Bio.，2004，49(6):941－959.

［7］ ZOUY，PAN X.Image reconstruction on PI－linesby useof filtered backprojection in helical cone－beam CT［J］.Phys.Med.Bio.，2004，49(12):2717－2731.

［8］ ZOU Y，PAN X，SIDKY E.Image reconstruction in region－of－interest from truncated projections in a reduced fan－beam scan［J］.Phys.Med.Bio.，2005，50(1):13－27.

［9］ YE Y，YU H，WEI Y，et al.A general local reconstruction approach based on truncated［EB/OL］.［2013－06－10］.http://www.hindawi.com/journals/ijbi/2007/063634/abs/.

［10］ PAN X，ZOU Y，XIA D.Image reconstruction in peripheral and central region－of－interestand data redundancy［J］.Med.Phys.，2005，32(3):673－684.

［11］ ZOU Y，PAN X，SIDKY E.Theory and algorithm for image reconstruction on chords and within regions of interest［EB/OL］.［2013－06－10］.http://www.ncbi.nlm.nih.gov/pubmed/16304723.

［12］谷建偉，張麗，陳志強(qiáng)，等.基于平行PI線的平行束CT重建［J］.清華大學(xué)學(xué)報(bào)，2007，47(3):393－395.

［13］喬志偉.能消除CT角誤差的斜變?yōu)V波器長度選擇方法［J］.計(jì)算機(jī)工程與設(shè)計(jì)，2012，33(4):1486－1490.

［14］李亮，陳志強(qiáng)，張麗，等.潘曉川教授的反投影濾波(BPF)新型重建算法介紹［J］.CT 理論與應(yīng)用研究，2006，15(3):68－73.

［15］李錚.CT重建算法的理論及綜合應(yīng)用的研究［D］.沈陽:東北大學(xué)，2005.