楊宏成，高 欣，張 濤1，

(1.中國科學(xué)院長春光學(xué)精密機(jī)械與物理研究所，吉林長春130000;2.中國科學(xué)院大學(xué)，北京100049;3.中國科學(xué)院蘇州生物醫(yī)學(xué)工程技術(shù)研究所，江蘇蘇州215163)

數(shù)字減影血管造影(digital subtraction angiography，DSA)是一種通過計算機(jī)把血管造影片上的骨與軟組織的影像消除，僅在影像片上突出血管的X射線技術(shù).近年來隨著平板探測器(flat panel detectors，F(xiàn)PD)的引入，錐束DSA可獲取血管3D影像信息.錐束DSA采用短掃描方式，可有效減少掃描時間、提高時間分辨率、降低射線劑量.但受掃描方式限制，錐束DSA很難獲得類似CT的圖像質(zhì)量.針對短掃描方式的錐束重建算法，目前主要有兩大類:一類是基于FDK的算法［1］，一類是基于Grangeat的算法［2］.與基于 Grangeat的算法相比，基于 FDK的算法重建速度快，空間分辨率高.但由于圓軌跡的短掃描方式不滿足完全精確重建條件，基于FDK的短掃描重建算法存在錐束偽影.

目前改善錐束偽影的方法主要有以下兩種:① 采用“圓弧 + 圓弧”［3］、“ 圓弧 + 線段”［4］、“圓弧+螺旋”［5］等掃描方式以滿足完備重建條件，利用精確重建算法來進(jìn)行圖像重建;② 采用改進(jìn)的FDK算法，很多學(xué)者提出了許多FDK類改進(jìn)算法，如 CW-FDK［6］、CB-FBP［7］、ACE［8-9］等.這些改進(jìn)方法的中心思想是通過一些操作來減少錐角增大時不精確重建形成的偽影，這些算法對錐角較小(小于±5°)情況有較高的重建質(zhì)量.第1種方法要求增加軌跡精度高，增加了DSA機(jī)械設(shè)計難度;第2種方法當(dāng)錐角大于5°時，改善效果有限.

文中擬提出一種針對DSA圓軌跡短掃描的FDK錐束重建算法.根據(jù)Radon空間數(shù)據(jù)缺失特點設(shè)計一種全新的反投影權(quán)重函數(shù)，將其引入圓軌跡短掃描錐束重建中.先分析錐束偽影形成的原因，提出一種反投影權(quán)重，并將其作為約束條件引入到Parker-FDK算法中實現(xiàn)抑制偽影的目的.采用無噪聲投影數(shù)據(jù)、有噪聲投影數(shù)據(jù)和自行研制的錐束DSA設(shè)備上采集的數(shù)據(jù)對新算法進(jìn)行試驗，列出該算法與Parker-FDK算法試驗對比結(jié)果，并最終對算法的的性能給出總結(jié).

1 算法及理論分析

1.1 FDK算法

FDK算法是 L.A.Feldkamp等人［10］把二維扇束重建公式推廣到三維空間的錐束重建算法.圖1表示錐束圓形掃描軌跡.

圖1 圓軌跡短掃描錐束投影示意圖

FDK重建公式可以表示為

式中:g(λ;u，v)為從平面探測器上獲得的投影數(shù)據(jù);λ為旋轉(zhuǎn)角度;(u，v)為探測器系統(tǒng)下坐標(biāo)值;γ為半錐角;D為射線源到探測器中心的垂直距離;h(u)為濾波函數(shù);?表示一維卷積.

Liu Ying等［11］將 Parker窗函數(shù)引進(jìn)到 FDK 算法中，提出了針對［0，π +2γmax］角度范圍內(nèi)圓軌道掃描的Parker-FDK重建算法.該方法利用類似FDK的方法進(jìn)行重建.其中，γmax表示最大半扇角.重建過程包括加權(quán)、數(shù)據(jù)冗余補(bǔ)償(Parker窗加權(quán))、濾波和反投影，公式表示如下:

式中:w(λ，γ)為 Parker窗函數(shù)，表示為

圖2為圓軌跡掃描Radon數(shù)據(jù)缺失示意圖.

圖2 圓軌跡掃描Radon數(shù)據(jù)缺失示意圖

由圖2可見，圓軌跡掃描方式只能采集到部分Radon數(shù)據(jù)，在旋轉(zhuǎn)軸z軸方向上存在著較大的陰影區(qū)域，整個區(qū)域類似面包圈結(jié)構(gòu).而FDK類算法將這些陰影區(qū)域填充為0，這就是FDK類算法在大錐角區(qū)域產(chǎn)生錐束偽影的主要原因.故對陰影區(qū)域進(jìn)行補(bǔ)償是減少錐束偽影的一個思路.

1.2 反投影權(quán)重錐束重建算法

首先引入共軛射線和共軛射線權(quán)重.圖3為共軛射線示意圖.

圖3 共軛射線示意圖

圖3中，重建平面內(nèi)任意一點P，掃描軌跡上的某一點S，兩者之間的連線在掃描軌跡平面內(nèi)的垂直投影的延長線與掃描軌跡的交點為S'，則S'P稱為SP的共軛射線，S'則是S關(guān)于P點的共軛點.以α表示射線SP的錐角，α'表示共軛射線S'P的錐角，一般情況下，α≠α'.考慮共軛射線的連續(xù)性，對直接射線和其共軛射線采用不同因子進(jìn)行加權(quán)處理，記作共軛射線權(quán)重.對于一組共軛射線來說，錐角小的那根射線投影的影響因子大一些.故可設(shè)共軛射線權(quán)重為

實際上，F(xiàn)DK算法是二維扇束重建公式在三維空間推廣，并沒有考慮離掃描軌道距離z不同的重建平面對反投影的影響.為此文中在FDK算法的反投影過程中引入基于距離權(quán)重的補(bǔ)償因子，記作反投影權(quán)重(backprojection weight，bpw)，對 Radon 缺失數(shù)據(jù)進(jìn)行補(bǔ)償，從而改善錐束偽影.

綜上所述，反投影權(quán)重應(yīng)該是錐角α和距離z的函數(shù)，記作wbpw(α，z)，可表示為

式中:k為(0，1)區(qū)間的松弛因子;zmax為軸向最大重建距離.x=-25 mm處的20°和200°時的反投影權(quán)重如圖4所示.可以看出反投影權(quán)重隨著|z|的增大而變大.

圖4 x=-25 mm時20°和200°的反投影權(quán)重

由此文中提出了基于FDK的反投影權(quán)重錐束重建算法(back projection weighted FDK，BPWFDK).在公式2中引入反投影權(quán)重，反投影權(quán)重重建算法可表示為

式中:h(u)為濾波函數(shù);w(λ，γ)為Parker權(quán)重函數(shù).應(yīng)該注意的是:短掃描方式下，部分直接射線的共軛射線不在掃描軌跡上，如圖5所示S'.這種情況，直接射線的反投影權(quán)重設(shè)置為1.

圖5 短掃描方式下的共軛射線示意圖

2 試驗與結(jié)果

2.1 無噪聲投影重建

文中基于三維Shepp-Logan頭模型模擬投影數(shù)據(jù)，對Parker-FDK和BPW-FDK算法進(jìn)行了比較.X射線源距離重建中心、探測器中心的垂直距離分別為400mm和600 mm，最大錐角和扇角都為18°，探測器陣列數(shù)為512×512，尺寸為0.39 mm×0.39 mm.圖6給出了理論模型以及Parker-FDK和BPW-FDK算法的重建結(jié)果.切片圖像大小為256×256，窗口顯示范圍為［1，1.05］.

圖6 x=-25 mm時的無噪聲重建結(jié)果切片圖像比較

從圖6可以看出:在投影錐角較大的區(qū)域，Parker-FDK算法重建圖像出現(xiàn)了明顯的錐束偽影;而BPW-FDK算法精確地重建出原物體.圖7給出了圖6重建的Sheep-Logan圖像在x=-25 mm，y=0 mm和z=35 mm，x=0 mm處剖面灰度值比較.圖7采用密度曲線的方法定量地顯示上述重建結(jié)果，可以看出與Parker-FDK算法相比，BPW-FDK算法所獲得的重建結(jié)果更好地擬合了模型的密度曲線.

圖7 圖6重建的Sheep-Logan圖像的剖面灰度值比較圖

表1給出了無噪聲情況Shepp-Logan頭模型之間的兩個誤差測度:歸一化均方距離判據(jù)d、歸一化平均絕對距離判據(jù)r作為重建質(zhì)量好壞的兩個標(biāo)準(zhǔn).d，r定義如下:

式中:tu，v和ru，v分別為模型和重建后圖像中第u行、v列的像素密度;ˉt為模型密度的平均值.從表1可以看出，BPW-FDK算法歸一化所獲得結(jié)果的均方距離判據(jù)d和歸一化平均絕對距離判據(jù)r均比Parker-FDK算法所獲得的結(jié)果小，表明誤差更小，重建圖像質(zhì)量更好.

表1 無噪聲重建誤差比較

2.2 有噪聲投影重建

為了檢測該算法的抗噪能力，在投影數(shù)據(jù)中加入0.1%的高斯隨機(jī)噪聲.圖8是加入噪聲后，用Parker-FDK和BPW-FDK算法重建的x=-25 mm的切片圖像.仿真參數(shù)同無噪聲投影情況下相同，切片圖像窗口顯示范圍為［1，1.05］.重建結(jié)果表明當(dāng)存在噪聲時，BPW-FDK算法仍能夠取得較好的結(jié)果.變化松弛因子k仍可以彌補(bǔ)當(dāng)錐角較大時出現(xiàn)的錐束偽影問題.

圖8 x=-25 mm時的有噪聲重建結(jié)果切片圖像比較

圖9給出了圖8重建的Sheep-Logan圖像在x=-25 mm，y=0 mm和z=35 mm，x=0 mm處剖面灰度值比較.從圖中可以看出當(dāng)存在噪聲時，與Parker-FDK算法相比，BPW-FDK算法所獲得的重建結(jié)果仍更逼近模型的密度曲線.

圖9 圖8重建的Sheep-Logan圖像的剖面灰度值比較圖

表2給出了有噪聲情況Shepp-Logan頭模型之間的兩個誤差測度.從表2可以看出，BPW-FDK算法歸一化均方距離判據(jù)d、歸一化平均絕對距離判據(jù)r均比Parker-FDK算法所獲得的結(jié)果小，表明重建圖像質(zhì)量更好，證明算法對有噪聲數(shù)據(jù)的重建有效.

表2 有噪聲重建誤差比較

2.3 DSA數(shù)據(jù)重建

為了更好地驗證BPW-FDK算法的性能，采用自行研制的DSA數(shù)據(jù)進(jìn)行重建試驗.DSA成像系統(tǒng)結(jié)構(gòu)如下:X射線源距離重建中心、探測器中心的垂直距離分別為756.4 mm和1 060 mm，掃描角度為202°，最大錐角和扇角都為12°，探測器陣列數(shù)為768×1 024，探測器單元尺寸為0.388 mm×0.388 mm.分別使用 Parker-FDK、BPW-FDK算法對其重建.圖10分別給出Parker-FDK和BPW-FDK算法重建的x=0mm處的切片圖像，切片圖像大小為256×256.重建物體為頭骨模型，右下角為面顱，具有豐富的細(xì)節(jié)信息.可以看出，與Parker-FDK算法相比，BPW-FDK算法在右下角處通過對遠(yuǎn)離掃描平面進(jìn)行補(bǔ)償，可重建出更多細(xì)節(jié)，標(biāo)明算法對實際DSA數(shù)據(jù)的重建是有效.

圖10 DSA數(shù)據(jù)重建結(jié)果比較

BPW-FDK算法保持了濾波反投影結(jié)構(gòu)，而濾波反投影算法中加權(quán)反投影運(yùn)算量占總運(yùn)算量的90%以上.BPW-FDK算法中反投影權(quán)重的引入并沒有增加反投影復(fù)雜度，與Parker-FDK算法相比，計算效率并沒有降低.

3 結(jié)論

文中分析了圓軌跡錐束掃描陰影區(qū)域形成的原因及對圖像重建的影響，提出了基于FDK的反投影權(quán)重錐束重建算法，結(jié)論如下:

1)BPW-FDK算法保持了濾波反投影結(jié)構(gòu)，反投影復(fù)雜度并沒有增加，保持了Parker-FDK算法重建速度快的優(yōu)點.

2)提出了基于距離權(quán)重的陰影區(qū)域補(bǔ)償方法，設(shè)計了一種新型的權(quán)重函數(shù)，在反投影過程中將基于距離的反投影權(quán)重作用于濾波后的投影數(shù)據(jù)，重建出圖像值.

3)仿真試驗表明，與Parker-FDK算法相比，BPW-FDK算法重建圖像的歸一化均方距離判據(jù)和歸一化平均絕對距離判據(jù)均降低了5%.

4)文中提出的BPW-FDK算法能夠較好地重建出錐束短掃描下低對比度物體，減少錐束偽影，提高圖像重建質(zhì)量.

References)

［1］ Li Liang，Xing Yuxiang，Chen Zhiqiang，et al.A curve-filtered FDK(C-FDK)reconstruction algorithm for circular cone-beam CT ［J］.Journal of X-Ray Science and Technology，2011，19(3):355-371.

［2］ Zhu L，Starman J，F(xiàn)ahrig R.An efficient estimation method for reducing the axial intensity drop in circular cone-beam CT［J］.International Journal of Biomedical Imaging，2008，doi:10.1155/2008/242841.

［3］ Zambelli J，Nett B E，Leng S，et al.Novel C-arm based cone-beam CT using a source trajectory of two concentric arcs［C］∥Proceedings of SPIE.Bellingham:SPIE，doi:10.1117/12.713813 2007.

［4］ Zamyatin A A，Katsevich A，Chiang B S.Exact image reconstruction for a circle and line trajectory with a gantry tilt［J］.Physics in Medicine and Biology，2008，doi:10.1088/0031-9155/53/23/N02.

［5］ Schomberg H，van de Haar P，Baaten W.Cone-beam CT using a C-arm system as front end and a spherical spiral as source trajectory［C］∥Proceedings of SPIE.Bellingham:SPIE，doi:10.1117/12.811545.

［6］ 陳 煉，吳志芳，王欽婭.錐束 CT的分區(qū)短掃描FDK重建算法［J］.清華大學(xué)學(xué)報:自然科學(xué)版，2009，49(6):860-863.Chen Lian，Wu Zhifang，Wang Qinya.Segmentation short scan FDK reconstruction algorithm for cone-beam CT［J］.Journal of Tsinghua University:Science and Technology，2009，49(6):860-863.(in Chinese)

［7］ Tang X Y，Hsieh J，Hagiwara A，et al.A three-dimensional weighted cone beam filtered backprojection(CBFBP)algorithm for image reconstruction in volumetric CT under a circular source trajectory ［J］.Physics in Medicine and Biology，2005，doi:10.1088/0031-9155/50/16/016.

［8］ Zhu L，Yoon S，F(xiàn)ahrig R.A short-scan reconstruction for cone-beam CT using shift-invariant FBP and equal weighting［J］.Medical Physics，2007，34(11):4422-4438.

［9］ Nett B E，Zhuang T L，Leng S，et al.Arc based conebeam reconstruction algorithm using an equal weighting scheme［J］.Journal of X-Ray Science and Technology，2007，15(1):19-48.

［10］ Feldkamp L A，Davis L C，Kress JW.Practical conebeam algorithm［J］.Journal of the Optical Society of A-merica A，1984，1(6):612-619.

［11］ Liu Ying，Liu Hong，Wang Ying，et al.Half-scan cone-beam CT fluoroscopy with multiple X-ray sources［J］.Medical Physics，2001，28(7):1466-1471.