[摘要]針對錐束CT錐角增大而產(chǎn)生的錐束偽影問題，比較不同三維重建算法的重建效果。本文采用三維加權(quán)重建算法來重建高質(zhì)量的小麥切片圖像。該算法結(jié)構(gòu)簡單，易于實(shí)現(xiàn)，不需要將原錐形束投影數(shù)據(jù)重新排列成錐形平行束數(shù)據(jù)，只需在反投影過程中增加三維加權(quán)函數(shù)，以減小軸向密度下降，提高FDK重建算法的精度。使用改進(jìn)Shepp-Logan頭模型和小麥數(shù)據(jù)來評估兩個(gè)算法的性能，結(jié)果表明這種算法在降低軸向失真方面取得了較好的效果，提高了FDK算法的重建質(zhì)量。

[關(guān)鍵詞]錐束CT;FDK重建;反投影

中圖分類號：TP751.1

文獻(xiàn)標(biāo)識碼：A

DOI：10.16465/j.gste.cn431232ts.20181227

近年來，計(jì)算機(jī)斷層掃描技術(shù)（CT）已廣泛應(yīng)用于醫(yī)學(xué)成像和工業(yè)無損檢測方面[1]。其中對于標(biāo)準(zhǔn)醫(yī)療CT系統(tǒng)和錐束形掃描儀，圓軌跡的錐束掃描模式發(fā)展迅速。目前，在網(wǎng)形掃捕的錐束CT系統(tǒng)中，主要的重建算法有迭代重建算法和解析類重建算法[2]。其中，解析重建算法以Radon變換為理論基礎(chǔ)，重建速度快，易于實(shí)現(xiàn)，能夠保證很好的重建精度，所以在實(shí)際系統(tǒng)中應(yīng)用廣泛。

FDK重建算法作為典型的解析圖像重建算法，圓軌跡的掃描方式不滿足數(shù)據(jù)充分（Data SufficiencyCondition，DSC）條件[3]，該算法中錐角的增加會導(dǎo)致錐束偽影，使圖像失真。目前，有許多種方法可以改善錐束偽影，主要可分為4類：

（1）通過改進(jìn)圓形軌跡來滿足DSC條件，如采用“圓弧+圓弧”“圓弧+線段”“網(wǎng)弧+螺旋”等掃捕軌跡l41。

（2）通過將錐束投影重排成平行投影后進(jìn)行重建，如T-FDK、CW-FDK、CB-FBP、C-FDK[5]。

（3）基于Grangeat公式對缺失的數(shù)據(jù)進(jìn)行估計(jì)[6]。

（4）利用反投影三維加權(quán)補(bǔ)償缺失數(shù)據(jù)，如X-FDK、FAFDK等[7]。

第1類方法需要高精度的軌跡和大量的計(jì)算，增加了機(jī)械設(shè)計(jì)難度，實(shí)際應(yīng)用中可行性低;第2類方法需進(jìn)行數(shù)據(jù)重排，會降低空間分辨率，造成損失，僅在錐角很小時(shí)比較有用;第3類方法在灰度下降比FDK更低的情況下，效率更高，但需要大量的計(jì)算;第4類方法具有簡單和直觀的結(jié)構(gòu)，能減少錐形束偽影。

針對FDK算法的錐角效應(yīng)，本文采用一種加權(quán)FDK算法，在反投影階段加入權(quán)重函數(shù)，來實(shí)現(xiàn)抑制偽影，達(dá)到提高重建圖像質(zhì)量的目的。通過對頭模型和實(shí)際小麥數(shù)據(jù)進(jìn)行實(shí)驗(yàn)，對實(shí)驗(yàn)結(jié)果進(jìn)行對比，以驗(yàn)證算法的有效性和實(shí)用性。

1 算法及理論分析

1.1 FDK算法及缺失數(shù)據(jù)

FDK算法由Feldkamp等于1984年提出，是迄今為止圓軌道錐束重建中最實(shí)用的一種算法，主要由于其一維移不變?yōu)V波和反投影的結(jié)構(gòu)[8-9]。FDK算法假定使用一個(gè)有限尺寸的等距平板探測器，傳統(tǒng)的錐束掃描結(jié)構(gòu)見圖1。

射線源S圍繞Z軸做網(wǎng)周旋轉(zhuǎn)，旋轉(zhuǎn)半徑也就是射線源到樣本中心的距離為R，當(dāng)待重建點(diǎn)在中心平面soy'時(shí)，這是等距扇形束的情況，因此可采用其重建公式進(jìn)行重建。然而，當(dāng)重建點(diǎn)不在中心平面時(shí)，如圖1中重建點(diǎn)P位于傾斜平面So'Q上，點(diǎn)Q（a，b）是點(diǎn)P在虛擬探測器上對應(yīng)的投影點(diǎn)，上面提到三維錐束FDK重建算法是二維扇束重建算法的一種擴(kuò)展，因此，對于傾斜平面So'Q應(yīng)用扇束公式重建之前，對應(yīng)的旋轉(zhuǎn)半徑和旋轉(zhuǎn)角度修正：

R'=R2+b2

dβ'=R/√R2+b2dβ

（1）

將R'和dβ'分別代替二維扇形束重建公式中的R和dβ，可得到三維錐束重建公式。FDK重建算法的具體實(shí)現(xiàn)可歸納為3個(gè)步驟：

第一步，通過加權(quán)因子來修正二維投影數(shù)據(jù)p（β，a，b）： p'（β，a，b）=R/√R2+a2+b2p（β，a，b） （2）

加權(quán)因子R/√R2+a2+b2的幾何意義是任何一條射線與中心射線之間夾角的余弦。

第二步，沿行方向?qū)πＵ蟮耐队皵?shù)據(jù)進(jìn)行濾波： p（β，a，b）=p'（β，a，b）×h（a）

（3） 第三步，加權(quán)反投影： fp（x，y，z）=∫02πR2/υρ（β，a，b）dβ

（4） 式中，U=R+xcosβ+ysinβ、a=R-xsinβ+ycosβ/R+xcosβ+ysinβ、b=Rz/R+xcosβ+ysinβ。

作為近似重建算法，F(xiàn)DK算法在小錐角的情況下，重建效果顯著。然而在大錐角情況下，由于圓形掃描模型白有的Radon空間數(shù)據(jù)缺失，在距中心平面有大的縱向距離的周邊區(qū)域中，會產(chǎn)生灰度降低的問題，影響重建質(zhì)量[10]。

通過將錐形束數(shù)據(jù)轉(zhuǎn)換成Radon數(shù)據(jù)，在三維Radon空間中填充類似環(huán)形的區(qū)域。由于存在空區(qū)域（圓環(huán)孔），造成三維逆Radon變換數(shù)據(jù)缺失，這是由于圓形軌跡未能滿足三維重建的數(shù)據(jù)充分條件。圖2是Radon域中的垂直平面二維圖，平面中的測量數(shù)據(jù)由兩個(gè)圓盤組成，而外部數(shù)據(jù)（如陰影區(qū)域所示）丟失。FDK算法在軌道平面上是精確的，而離軌道平面越遠(yuǎn)，Radon數(shù)據(jù)缺失越嚴(yán)重（FDK算法假設(shè)未測量數(shù)據(jù)為零），灰度和強(qiáng)度下降越明顯，重建誤差越大[11-12]。對于FDK中缺失數(shù)據(jù)的問題，研究表明軸向的偏差可以描述為帽式函數(shù)（在其一階導(dǎo)數(shù)是遞增函數(shù)），并且利用典型的經(jīng)驗(yàn)函數(shù)如余弦函數(shù)和高斯函數(shù)來補(bǔ)償軸向的缺失數(shù)據(jù)是可行的。此時(shí)，反投影過程中的權(quán)重相當(dāng)于帽形函數(shù)，而不是FDK中0.5的固定權(quán)重。

1.2 三維加權(quán)錐束重建算法

由于錐束重建的衰減效應(yīng)可以用余弦函數(shù)建模，可以用余弦的倒數(shù)表示權(quán)重函數(shù)，對于軸向（橫跨z平面）的衰減效應(yīng)，權(quán)重函數(shù)應(yīng)包含變量z，則Weight-FDK算法中提出采用參數(shù)C1和Cc2進(jìn)行調(diào)節(jié)縱向距離變量z和徑向距離變量r，從而調(diào)整權(quán)重函數(shù)的形狀，其數(shù)學(xué)表達(dá)式為： f（x，y，z）=∫02πR2/U×WW-FDK3D（z，c1 ，C2）

xp（β，a，b）dβ

（5）

各個(gè)變量的含義與方程式（4）相同，其中，用于補(bǔ)償衰減退化的經(jīng)驗(yàn)權(quán)重公式為：

WW-FDK3D（Z， C1， c2）=1/cos[c1z/（R-c2r）]

（6）

與FDK不同，經(jīng)驗(yàn)權(quán)重在非中心平面上總是不相等且大于0.5，其應(yīng)該滿足不等式：

W3DW-FDK（z，c1，C2）+ W3DW-FDK（z，cl，C2）≥1（7）

其中r=√x2+y2 +z2，對于帽形函數(shù)的描述，參數(shù)Cc1和c2的選擇應(yīng)該滿足以下條件[13]：

1≤c1z÷（R-c2r）lt;π/2

（8）

通過設(shè)置參數(shù)c1和c2，包括如下特殊情況：

W3DW-FDK（Z，c1，o）=W3DW-FDK（z，c1）1/cos（c1z÷R）

（9）

W3DW-FDK（o，c1，c2）=1

（10）

W3W-FDK（z，o，c2）=1

（11）

式（9）僅表示縱向衰減效應(yīng)的權(quán)重，式（10）代表對中心平面沒有補(bǔ)償，式（11）表示對整個(gè)數(shù)據(jù)體沒有補(bǔ)償。

該算法只需在反投影階段考慮軸向距離的影響而引入三維權(quán)重函數(shù)，通過人為經(jīng)驗(yàn)調(diào)節(jié)兩個(gè)參數(shù)來改變權(quán)重的大小，以彌補(bǔ)數(shù)據(jù)缺失，從而改善錐束偽影。該算法操作簡單，易于實(shí)現(xiàn)，并且在中心平面相當(dāng)于FDK算法。

1.3 評價(jià)指標(biāo)

為了進(jìn)一步比較不同算法的重建質(zhì)量，這里引入峰值信噪比[14]、均方誤差[15]、結(jié)構(gòu)相似度作為圖像重建質(zhì)量的客觀評價(jià)指標(biāo)，表達(dá)式為：

（1）峰值信噪比（peak signal to noise ratio）：

PSNR=10lg2552/1/MNM∑i=1N∑j=1[X（x，j）-Y（i，j）]2

（12）

X和Y分別表示待評價(jià)圖像和參考圖像，圖像大小都是M×N。該指標(biāo)數(shù)值越大，表示重建圖像失真越小，質(zhì)量越好。

（2）均方根誤差（roommean squared error）：

RMSE=√1/MNM∑i=1N∑j=1[X（i，j）-Y（i，j）]2

RMSE表征測量值與真實(shí)值之間的偏差，該指標(biāo)越小，表示圖像質(zhì)量越好。

（3）結(jié)構(gòu)相似度（structure similarity）：

SSIM（x，y）=（2μxμy+c1）（2σxy+c2）/（μx2+μy2+c1）（σx2+y2+c2）

（14）

μx、μy、σx、σy、σxy分別表示圖像X、y的均值、方差和協(xié)方差。c1、c2為常數(shù)，是衡量兩幅圖像之間相似性的新指標(biāo)，該指標(biāo)的值越大越好，最大值為1。

2 結(jié)果與分析

2.1 頭模型仿真重建結(jié)果

這里以三維Shepp—Logan改進(jìn)頭模型的投影數(shù)據(jù)為對象，比較FDK算法和Weight-FDK算法的重建效果。其中，網(wǎng)軌跡掃描半徑是300mm，從光源到探測器的距離為600mm，投影數(shù)據(jù)是從360。圓周上均勻采集360個(gè)投影。

為了在Weight-FDK算法中獲得最佳參數(shù)，首先選擇c2=0并調(diào)整c1以校正沿x=y=z的強(qiáng)度下降，確定了c1，并調(diào)整C2來校正殘差作為徑向距離的函數(shù)。x=Omm處頭模型如圖3（a）所示，F(xiàn)DK算法和Weight-FDK算法在y-z平面的重建效果如圖3（b）和圖3（c）所示，切片圖像大小為256×256mm。

由圖3可見，在投影錐角相對較大的區(qū)域產(chǎn)生了錐角效應(yīng)，而Weight-FDK算法可以精確地重建三維頭模型。圖3（a）、3（b）、3（c）對應(yīng)的Profile曲線對比如圖3（d）所示，由圖3可見，Weight-FDK算法重建結(jié)果和原模型更接近。為了更加客觀地評價(jià)不同算法的重建質(zhì)量，這里通過RMSE、PSNR、SSIM 三個(gè)參數(shù)進(jìn)行評價(jià)，其評價(jià)結(jié)果如表1所示。據(jù)表1可知，Weight-FDK算法重建結(jié)果的PSNR和SSIM比FDK算法的分別大1.2217和0.0078，而RMSE比FDK算法小2.6477，表明Weight- FDK算法的重建結(jié)果的誤差更小，重建圖像更清晰。

2.2 實(shí)際數(shù)據(jù)重建結(jié)果

為了驗(yàn)證Weight-FDK算法的性能，將該算法應(yīng)用于實(shí)際的小麥數(shù)據(jù)重建中，光源到探測器的距離為282.9mm，到重建物體的距離為30.6mm，投影數(shù)據(jù)是在360°的圓周上均勻采集1000個(gè)投影。分別采用FDK算法和Weight-FDK算法進(jìn)行重建，得到兩個(gè)算法在x=Omm處的重建切片圖像如圖4所示。重建物體為小麥實(shí)體，可以看出與FDK算法相比，Weight-FDK算法重建結(jié)果在兩端對比度更高、圖像更清晰，表明算法對實(shí)際數(shù)據(jù)是有效的，可以明顯改善小麥圖像質(zhì)量。圖4在x=Omm，y=14mm處的局部Profile曲線如圖5所示。由圖5司知，Weight-FDK算法小麥重建結(jié)果在坐標(biāo)為[70，140]和[360，460]，小麥區(qū)域的灰度值較FDK算法重建圖像有明顯增大，重建效果有了明顯的提升，表明該算法具有很大的實(shí)際應(yīng)用價(jià)值。

3 結(jié)論

本文分析了網(wǎng)軌跡錐束掃捕陰影區(qū)域產(chǎn)生的原因，并采用Weight-FDK算法提高了圖像的重建質(zhì)量，該算法不用對錐束投影數(shù)據(jù)進(jìn)行重排，保證了重建圖像的分辨率，易于實(shí)現(xiàn)，且在中心平面等效于FDK算法。通過對比實(shí)驗(yàn)，對于改進(jìn)三維Shepp-Logan頭模型和實(shí)際小麥投影數(shù)據(jù)，不難發(fā)現(xiàn)Weight-FDK算法能夠明顯改善FDK算法的錐束偽影，可以顯著提高圖像的重建質(zhì)量，且該技術(shù)可擴(kuò)展到其他谷物投影數(shù)據(jù)的重建。

參考文獻(xiàn)

[1]何昌保，馬秀麗，余長明，基于Random Walks算法的心臟雙源CT左心房分割[J]電子測量技術(shù)，2016 （5）：75-79.

[2] TUY， HEANG K. An inversion formula for cone-beam reconstruction[J].SIAM Joumal on Applied Mathematics，1983（3）：546-552

[3]鄭喜艷.基于協(xié)同進(jìn)化遺傳算法的PET圖像重建研究[D].鄭州：鄭州大學(xué)，2016，

[4]李增光，李磊，韓玉，等.錐柬CT圓加直線軌跡反投影濾波重建算法[J]光學(xué)學(xué)報(bào)，2016 （9）：157-164.

[5]楊宏成，鬲欣，張濤.基于FDK反投影權(quán)重的錐柬DSA重建算法[J]江蘇大學(xué)學(xué)報(bào)（自然科學(xué)版），2013 （2）：190-195.

[6]王革，SEUNG WOOKLEE.論錐束CT掃描Grangeat-型Katsevich-型的算法（英文）[J].CT理論與應(yīng)用研究，2003 （2）：45-55.

[7]楊亞飛，張定華，黃魁東，等.基于橢球包圍盒的錐束CT三維加權(quán)重建算法[J]儀器儀表學(xué)報(bào)，2016 （11）：2563-2571.

[8]張晶，張權(quán)，劉袆，等.錐束CT中FDK算法的快速對稱優(yōu)化與GPU實(shí)現(xiàn)[J].計(jì)算機(jī)工程與設(shè)計(jì)，2014 （8）：2813-2816

[9]張劍，陳志強(qiáng).三維錐形束CT威像FDK重建算法發(fā)展綜述[J].中國體視學(xué)與圖像分析，2005 （2）：116-121.

[10] HE CH B， MAXL，YU CH M. Study ofleft atrium segmentationin dual source CT image with random walks algorithms[J]Electronic Measurement Technology，2016（5）：75 -79

[11] HUANG K D， ZHANG H， SHI Y K，et al. Scatter correctionmethod for cone-beam CT based on interlacing-slit scan[J].Chinese Physics B，2014（9）：515-521

[12] ZHU L STARMAN J， FAHRIG R An efficient estimation methodfor reducing the axial intensity drop in circular cone-beam CT[J]International Joumal of Biomedical imaging，2014（1）：242 841

[13] GRIMMER R， OELHAFENM， ELSTROM U，et al. Cone-beam CT image reconstruction with extended z range[J].MedicalPhysics，2009（7）：3363-3370

[14]李英，楊秋翔，雷海衛(wèi)，等.基于深度圖像繪制技術(shù)的Criminisi算法的改進(jìn)[J].計(jì)算機(jī)工程與設(shè)計(jì)，2017 （5）：1287-1291.

[15]張小英蘋果近紅外光譜的預(yù)處理[J]紅外，2016 （5）：43-48