惠苗

（三明學(xué)院信息工程學(xué)院，福建三明365004）

一種含噪投影的螺旋錐束CTKatsevich重建算法

惠苗

（三明學(xué)院信息工程學(xué)院，福建三明365004）

在實(shí)際工程應(yīng)用中，投影體數(shù)據(jù)含有大量噪聲.Katsevich算法作為一種精確的重建算法，它可以得到質(zhì)量較高的重建圖像，但它需要對投影數(shù)據(jù)旋轉(zhuǎn)角和二維探測器坐標(biāo)做求導(dǎo)運(yùn)算，求導(dǎo)會放大噪聲，甚至使重建圖像失真.結(jié)合Katsevich算法中對投影數(shù)據(jù)的旋轉(zhuǎn)角避開求導(dǎo)和探測器坐標(biāo)求導(dǎo)采用五點(diǎn)差分方法代替直接求導(dǎo)，可以降低對旋轉(zhuǎn)角度的采樣率，減少對病人的射線輻射和射線成本，減小求導(dǎo)誤差和噪聲帶來的負(fù)面影響，從而提高圖像質(zhì)量.

螺旋錐束CT；圖像重建；Katsevich算法；投影數(shù)據(jù)求導(dǎo)

CT作為一種先進(jìn)的無損檢測技術(shù)，已廣泛應(yīng)用于醫(yī)學(xué)、考古、航空航天、安全檢測、農(nóng)業(yè)、設(shè)備制造等相關(guān)領(lǐng)域.其中，CT的掃描方式和重建算法是CT研究的關(guān)鍵技術(shù)和研究熱點(diǎn).在掃描方式上，目前主要應(yīng)用的有圓形掃描和螺旋掃描兩種，對于長物體采用螺旋掃描方式有明顯的優(yōu)勢，螺旋掃描時(shí)間短、射線利用率高、軸向分辨率高，而且它的掃描軌跡滿足三維精確重建的完全條件.在重建算法上，Katsevich算法是一種精確的重建算法，該算法重建圖像的質(zhì)量較高[1-3]，但該算法在重建的過程中，需要對投影數(shù)據(jù)求導(dǎo)，求導(dǎo)過程會放大噪聲，影響圖像的邊緣信息，甚至使重建圖像失真.針對這種問題，本文結(jié)合了避開對投影數(shù)據(jù)旋轉(zhuǎn)角求導(dǎo)和對探測器求導(dǎo)采用五點(diǎn)差分方法求導(dǎo).

1 螺旋錐束掃描方式

射線源掃描軌跡描述為：

式中R為射線源到物體中心的距離，也稱為旋轉(zhuǎn)半徑，h為螺距，s為旋轉(zhuǎn)角度.

感興趣區(qū)域重建坐標(biāo)點(diǎn)定義為X=(x,y,z),待重建圖像在點(diǎn) X的函數(shù)值表示為 f(X),區(qū)域x2+y2＜R稱為視場FOV（fieldofview）.在旋轉(zhuǎn)角s確定下，X點(diǎn)在探測器上投影的位置表示為(u,v)，投影數(shù)據(jù)值表示為g(s,u,v).

定義3個(gè)兩兩垂直的單位向量d1,d2,d3，用來表示平板探測器上錐束投影測量的局部坐標(biāo)系，這里 d1=(-sins,coss,0),d2=(0,0,1),d3=(-coss,-sins,0).

如圖1所示，射線源y(s)在平板探測器上的投影點(diǎn)是平板探測器的坐標(biāo)原點(diǎn)O,物體上任意點(diǎn)X投影為平板探測器上的坐標(biāo)記為(u,v).

圖1 平板探測器上錐束投影測量的局部坐標(biāo)系圖

2 Katsevich重建算法

根據(jù)Katsevich螺旋錐束CT的算法框架[4]，視場內(nèi)的點(diǎn)X處的函數(shù)灰度值 f(X)可通過點(diǎn)X的通過PI線段確定的螺旋軌道區(qū)間上的投影數(shù)據(jù)進(jìn)行濾波運(yùn)算，再經(jīng)過加權(quán)反投影得到.重建公式可寫為：

這里濾波數(shù)據(jù)可由下式確定：

式中，D為射線源到探測器平面中心的距離；h表示 Hilbert變 換 核 函 數(shù) ，定 義 為 ：h(t)=表示對投影數(shù)據(jù)g(s,u,v)求導(dǎo)運(yùn)算.

螺旋錐束Katsevich算法的實(shí)現(xiàn)過程分為投影數(shù)據(jù)求導(dǎo)、濾波、反投影重建三部分.

（I）投影數(shù)據(jù)求導(dǎo)：

（II）求導(dǎo)后的數(shù)據(jù)進(jìn)行希爾伯特變換：

（III）反投影重建：

其中sb(X)、st(X)表示由X點(diǎn)所確定的PI線的兩個(gè)端點(diǎn).

3 Katsevich重建算法回避投影角度求導(dǎo)

在式（4）的計(jì)算中，需要沿著旋轉(zhuǎn)軌道對投影數(shù)據(jù)g(s,u,v)進(jìn)行求導(dǎo)運(yùn)算，然而在數(shù)值計(jì)算中求導(dǎo)數(shù)運(yùn)算是通過差分運(yùn)算來逼近的，這使得求導(dǎo)運(yùn)算依賴于旋轉(zhuǎn)角度的采樣率和探測器單元的采樣率.其中螺旋旋轉(zhuǎn)角的采樣率遠(yuǎn)低于探測器單元的采樣率，而且螺旋旋轉(zhuǎn)角的采樣率越高，在醫(yī)學(xué)上病人射線輻射就越高，射線成本就越高，運(yùn)算量也很大.所以，考慮Katsevich的重建公式中避開對旋轉(zhuǎn)角度的求導(dǎo)，這將在醫(yī)學(xué)臨床有極其重要的意義[5-6].式(9)給出了避開旋轉(zhuǎn)角求導(dǎo)的重建公式.

式中：

實(shí)際的螺旋錐束成像系統(tǒng)中，隨著探測器的發(fā)展和采樣率逐步提高，使得投影數(shù)據(jù)對探測器變量u和v的數(shù)值求導(dǎo)運(yùn)算精度遠(yuǎn)高于投影數(shù)據(jù)對旋轉(zhuǎn)角的求導(dǎo)運(yùn)算，因此重建公式通過數(shù)學(xué)上的變形，避開了旋轉(zhuǎn)角度的求導(dǎo)運(yùn)算，這具有重要的實(shí)際意義.

4 投影數(shù)據(jù)的五點(diǎn)方法求導(dǎo)

螺旋錐束Katsevich算法的實(shí)現(xiàn)過程，需要對投影數(shù)據(jù)求導(dǎo)，有多項(xiàng)式擬合求導(dǎo)[7]，求導(dǎo)過程較復(fù)雜，本文采用五點(diǎn)差分方式求導(dǎo)，可以有效增強(qiáng)圖像邊緣.下面公式是投影數(shù)據(jù)對探測器u和v求導(dǎo)，其中定義旋轉(zhuǎn)角度采樣間隔為Δs,采樣數(shù)為K,每個(gè)探測器單元的寬度為Δu,探測的寬度數(shù)量為M,每個(gè)探測器單元的高度為Δv,探測器的高度數(shù)量為N,探測器的寬為M·Δu,探測器高為N·Δv,投影數(shù)據(jù)離散形式為 g(sk,um,vn)，這里 0≤k≤K-1，0≤m≤M-1，0≤n≤N-1.投影數(shù)據(jù)g(sk,um,vn)避開對旋轉(zhuǎn)角度的求導(dǎo)，對探測器橫向單元求導(dǎo)如下：

當(dāng)m=0時(shí)：

當(dāng)m=1時(shí)：

當(dāng)1＜m＜M-2時(shí)：

當(dāng)m=M-1時(shí)：

當(dāng)m=M-2時(shí)：

同理，可得投影數(shù)據(jù)g(sk,um,vn)對探測器縱向單元求導(dǎo).

實(shí)際采集的投影數(shù)據(jù)中，由于CT探測系統(tǒng)噪聲包括暗電流噪聲、散粒噪聲、熱噪聲等[8]，而熱噪聲和散粒噪聲是高斯白噪聲，這些噪聲會弱化邊緣信息，如果通過直接求導(dǎo)可以放大噪聲信息，而利用五點(diǎn)公式求導(dǎo)可以減少求導(dǎo)誤差，增強(qiáng)邊緣信息，減少噪聲對重建圖像質(zhì)量的影響.

5 實(shí)驗(yàn)數(shù)據(jù)及重建結(jié)果

5.1 實(shí)驗(yàn)數(shù)據(jù)

為驗(yàn)證以上算法有效性，對修訂的三維Shepp-Logan頭部模型利用解析法生成螺旋錐束投影數(shù)據(jù)，并對螺旋錐束投影數(shù)據(jù)添加高斯白噪聲，利用改進(jìn)的求導(dǎo)方法和改進(jìn)Katsevich重建方法錐束對Shepp-Logan頭部模型進(jìn)行重建，射線源到探測器的距離D為150像素，螺旋掃描軌跡的螺距為25像素，探測器的高和寬都為100像素，每圈采樣360次.

5.2 重建結(jié)果

以經(jīng)典Shepp-Logan頭部模型為例，在實(shí)驗(yàn)中添加高斯白噪聲，因?yàn)樵贑T探測系統(tǒng)中，噪聲主要成分有低頻噪聲與白噪聲[5].不僅通過實(shí)驗(yàn)用原始Katsevich重建算法重建圖像，還用以上介紹的結(jié)合回避投影角度求導(dǎo)以及五點(diǎn)公式插值多項(xiàng)式對探測器坐標(biāo)求導(dǎo)的Katsevich重建算法重建圖像，并對比兩種方法的重建圖像質(zhì)量.圖2為Shepp-Logan頭部模型重建圖像z=0切片圖，圖3為Shepp-Logan頭部模型x=0切片圖.圖2和圖3中，（a）為投影數(shù)據(jù)沒有添加噪聲Katsevich重建圖；（b）為投影數(shù)據(jù)沒有添加噪聲本文方法Katsevich重建圖；（c）為投影數(shù)據(jù)添加高斯白噪聲后Katsevich重建圖；（d）為投影數(shù)據(jù)添加高斯白噪聲后本文Katsevich重建圖.

5.3 圖像評價(jià)

在圖像中，某一方向的灰度級變化率大，它的梯度也就大.因此可以用平均梯度值來衡量圖像的清晰度，反映圖像對細(xì)節(jié)對比的表達(dá)能力，還同時(shí)反映出圖像中微小細(xì)節(jié)反差和紋理變換特征，所以本文用平均梯度來評價(jià)重建圖像的清晰度.平均梯度越大，圖像層次越多，清晰度越高.平均梯度定義如（17）式：

圖2 Shepp-Logan頭部模型重建圖像z=0切面

圖3 Shepp-Logan頭部模型重建圖像x=0切面

式中：F(i,j)為圖像的第i行、第j列的灰度值；M、N分別為圖像的總行數(shù)和總列數(shù).表1列出了分別用Katsevich直接求導(dǎo)法與本文Katsevich方法，對投影數(shù)據(jù)求導(dǎo)重建的切片的平均梯度值.

表1 圖2、圖3中各圖像的平均梯度

由表1可以看出，用本文方法的重建切片比原始Katsevich的平均梯度值大，重建圖像的清晰度更高.

6 結(jié)論

在實(shí)際CT成像系統(tǒng)中，獲得投影數(shù)據(jù)混有各種噪聲，然而Katsevich這種螺旋錐束精確重建算法對噪聲敏感，所以有必要對這種算法進(jìn)行改進(jìn)，提高它的抗噪能力，所以應(yīng)盡可能減少噪聲對重建圖像的影響.對噪聲敏感的原因是投影數(shù)據(jù)求導(dǎo)會放大噪聲的影響，并且需要進(jìn)行高采樣來達(dá)到更為精確的重建結(jié)果.本文結(jié)合了避開旋轉(zhuǎn)角求導(dǎo)和對探測器橫縱坐標(biāo)五點(diǎn)方法求導(dǎo)方式.避開旋轉(zhuǎn)角求導(dǎo)可以減少旋轉(zhuǎn)角的采樣，減少對病人的射線輻射.對探測器橫縱坐標(biāo)五點(diǎn)公式插值多項(xiàng)式求導(dǎo)方式，能夠較好地處理邊緣數(shù)據(jù)，從而使重建結(jié)果更加清晰，還能減少噪聲對重建質(zhì)量的影響.

[1]KATSEVICHA.Animprovedexactfilteredbackprojectionalgo?rithmforspiralcomputedtomography[J].AdvancesinApplied Mathematics,2004,32(4):681-697.

[2]YUHY,WANGG.StudiesonimplementationoftheKatsevichal?gorithmforspiralcone-beamCT[J].JournalofX-rayScienceand Technology,2004(12):97-116.

[3]KATSEVICHA.Ontwoversionsof3PIalgorithmforspiralCT[J]. PhysicsinMedicineandBiology,2004,49(11):2129-2143.

[4]YANGJS,GUOXH,KONGQ,etal.Parallelimplementationof Katsevich'sFBPalgorithm[J].InternationalJournalofBiomedical Imaging,2006:17463.doi:10.1155/IJBI/2006/17463.

[5]馬建華,陳武凡.不含旋轉(zhuǎn)角度微分的螺旋錐束CT重建[J].中國圖象圖形學(xué)報(bào),2008,4(13):647-653.

[6]NOOF,PACKJ,HEUSCHERD.Exacthelicalreconstructionus?ingnativecone-beamgeometries[J].PhysicsinMedicineandBiolo?gy,2003,48(23):3787-3818.

[7]曾理,牛小明,鄒曉兵.螺旋錐束CT重建Katsevich算法在含噪投影中的應(yīng)用[J].計(jì)算機(jī)工程與應(yīng)用,2007,43(26):187-189.

[8]王玨,陳教澤,譚輝,等.工業(yè)CT探測系統(tǒng)噪聲特性研究[J].核電子學(xué)與探測技術(shù),2010,7(30):929-934.

【編校：許潔】

ASpiralConeBeamCTKatsevichReconstructionAlgorithmforNoisedProjection

HUIMiao
(InstituteofInformationEngineering,SanmingUniversity,Sanming,Fujian365004,China)

Inpracticalengineeringapplications,theprojectiondatacontainsalotofnoise.TheKatsevichalgorithmisan exactimagingalgorithmwithgoodreconstructionimage.Thisalgorithmneedsderivationoperationonprojectiondatarota?tionangleandtwo-dimensionalcoordinatesofthedetector.Thederivationwillamplifythenoiseanddistorttherecon?structionimage.Combinedwiththeavoidderivationoftherotationangleandtheuseoffive-pointdifferentialmethodin?steadofdirectderivationtothedetectorcoordinate,thisalgorithmcanreducethesamplingrateofrotationangleandthe radiationtothepatientandX-raycost.Itcanlessentheerrorofthederivativeandthenegativeimpactcausedbythe noise.Thus,itimprovestheimagequality.

spiralconebeamCT;imagereconstruction;Katsevichalgorithm;projectiondataderivation

TP391.4

A

1671-5365(2015)12-0028-04

惠苗.一種含噪投影的螺旋錐束CTKatsevich重建算法[J].宜賓學(xué)院學(xué)報(bào),2015,15(12):28-31. HUIM.ASpiralConeBeamCTKatsevichReconstructionAlgorithmforNoisedProjection[J].JournalofYibinUniversity, 2015,15(12):28-31.

2015-04-21修回：2015-06-29

資金項(xiàng)目：福建省教育廳科技項(xiàng)目(JA12302)

惠苗（1980-），女，講師，研究方向?yàn)樾畔⑻幚砼c重建、無損檢測等

時(shí)間：2015-07-0110:15

http://www.cnki.net/kcms/detail/51.1630.Z.20150701.1015.002.html