基于雙域?yàn)V波與距離變換的螺旋CT圖像風(fēng)車偽影校正

黃銅港， 李保磊， 湯少杰， 徐圓飛，沈建冬， 屈軍鎖， 范九倫

(1．西安郵電大學(xué) 自動化學(xué)院， 陜西 710121; 2．北京航星機(jī)器制造有限公司， 北京 100013;3．自動分揀技術(shù)研究中心(西安郵電大學(xué))，國家郵政管理局， 陜西 710121; 4．西安市先進(jìn)控制與智能處理重點(diǎn)實(shí)驗(yàn)室， 陜西 710121;5．西安郵電大學(xué) 通信與信息工程學(xué)院， 陜西 710121)

0 引言

當(dāng)今計算機(jī)斷層掃描(CT)技術(shù)在眾多領(lǐng)域有著不可取代的作用。隨著CT技術(shù)的發(fā)展，多層螺旋CT機(jī)在醫(yī)療及安檢領(lǐng)域有著廣泛的應(yīng)用，但是該掃描方式會產(chǎn)生偽影，這是由于掃描的三維性質(zhì)和使用近似算法如FDK重建算法[1]及其衍生的三維算法[2-3]造成的。雖然Grangeat算法[4]是一種精確重建算法,但是實(shí)現(xiàn)比較繁瑣，目前在實(shí)際CT成像系統(tǒng)中很少有應(yīng)用，所以本文仍采用廣泛使用的FDK重建算法。風(fēng)車偽影(Windmill Artifact)[5]是螺旋CT掃描所產(chǎn)生的一種偽影，形成原因是在螺旋CT數(shù)據(jù)采集時z軸方向采樣率不足引起的。風(fēng)車偽影其表現(xiàn)為螺旋CT(Helical CT)重建圖像中在沿著高對比邊緣方向呈現(xiàn)風(fēng)車狀明暗相間的條紋，最常出現(xiàn)在z軸解剖結(jié)構(gòu)和密度變化很快的位置。這些偽影的產(chǎn)生會降低圖像質(zhì)量，甚至?xí)?dǎo)致錯檢、誤檢。

針對風(fēng)車偽影校正問題，許多科研人員進(jìn)行了大量研究來降低風(fēng)車偽影影響。例如，合理選取螺距[6-7]。螺距大小對風(fēng)車偽影的影響十分明顯，當(dāng)螺距值越大時螺旋一周所涉及的用于重建的探測器排數(shù)就越多，這樣就容易導(dǎo)致采樣不足而引起風(fēng)車偽影。但在實(shí)際的應(yīng)用中，小螺距的選取是要付出代價的，這種方法會大大增加投影數(shù)據(jù)的采集時間，同時圖像重建濾波和反投影的計算量成數(shù)倍增加。并且小螺距換來的也僅僅是等中心處的雙倍采樣，所以效果受到一定的限制，因此，簡單的降低螺距并不能完全克服風(fēng)車偽影。風(fēng)車偽影校正研究還包括層厚的合理選取[8]。風(fēng)車偽影的嚴(yán)重程度與層厚密切相關(guān)，層厚包括準(zhǔn)直器層厚和重建圖像層厚。

有效的提高z軸采樣率是風(fēng)車偽影校正最有效的方法。如果能將平面內(nèi)的飛焦(Flying Focal Spot，F(xiàn)FS)技術(shù)應(yīng)用到z軸上，即通過球管焦點(diǎn)在z軸方向上周期性運(yùn)動來提高采樣率，就可以有效地克服風(fēng)車偽影。西門子公司的z軸動態(tài)飛焦[9](z-Flying Focal Spot，z-FFS)技術(shù)可在任意位置獲得雙倍采樣，打破了小螺距僅在等中心處獲得雙倍采樣的缺陷，擺脫了螺距限制，可在大視場(FOV)條件下有效地克服風(fēng)車偽影。然而，焦點(diǎn)沿z軸擺動的實(shí)現(xiàn)比平面內(nèi)的擺動更具挑戰(zhàn)。Hsieh J提出了一種自適應(yīng)插值方法[10]，使用加權(quán)插值組合，在采樣不足的情況下進(jìn)行自適應(yīng)。該方法能有效減少偽影，同時,將對層面敏感度輪廓(Slice Sensitivity Profile, SSP)的影響保持在最低水平。除上述風(fēng)車偽影校正方法外，Xie等[11]提出了一種基于點(diǎn)陣采樣框架的軸向CT數(shù)據(jù)采集優(yōu)化，利用算法實(shí)現(xiàn)探測器交錯來抑制風(fēng)車偽影。Fang等提出了雙字典學(xué)習(xí)(ASDDL)抑制偽影的方法[12]來減少多層螺旋CT中的風(fēng)車偽影。Shi等提出了一種結(jié)合低秩和TV正則化(LRTV)的超分辨率重建方法[13]，以提高M(jìn)DCT的z軸分辨率。張韡提出了一種安檢CT重建圖像的處理方法和裝置[14]，該方法是將安檢CT寬排探測器采集到的投影數(shù)據(jù)進(jìn)行排間擴(kuò)充，用擴(kuò)充后的投影數(shù)據(jù)進(jìn)行重建，增加了z軸采樣率。上述幾種方法從校正效果來看，還存在一些不足且不能完全消除風(fēng)車偽影，仍有偽影殘留。

在現(xiàn)有技術(shù)中，可在獲得CT重建圖像前，預(yù)先設(shè)置掃描條件和重建條件，從而在CT掃描過程中，增加CT重建圖像的厚度，以此來避免風(fēng)車偽影的影響。但是，即使改變掃描條件和重建條件，也可能無法完全消除風(fēng)車偽影，從而得到包含風(fēng)車偽影的CT重建圖像?？紤]到上述方法的局限性，本文利用雙側(cè)濾波(Bilateral Filtering)與短時傅里葉變換(STFT)的雙域?yàn)V波[15-16]進(jìn)行螺旋CT圖像風(fēng)車偽影校正。由于風(fēng)車偽影特征是隨空間位置變化的，因此，本文嘗試同時基于雙域?yàn)V波和距離變換對螺旋CT重建圖像風(fēng)車偽影進(jìn)行校正。

1 方法

1.1 雙域?yàn)V波

雙域?yàn)V波執(zhí)行過程具有非迭代和局部性特點(diǎn)，包括雙側(cè)濾波與STFT兩部分[15]。雙側(cè)濾波既以中心像素附屬鄰域內(nèi)像素間空間位置差異設(shè)計加權(quán)，又結(jié)合了像素值差異設(shè)計加權(quán)，可有效恢復(fù)高對比度細(xì)節(jié)；STFT可有效恢復(fù)低對比度細(xì)節(jié)。與其他濾波方法相比，雙域?yàn)V波主要優(yōu)點(diǎn)是能夠以合理計算量進(jìn)行圖像去噪，同時能夠較好地保持圖像細(xì)節(jié)[16-17]。

雙域?yàn)V波給定一個未知CT圖像x，假定圖像采集過程引入方差為σ2的干擾η，得到含干擾圖像y：

y=x+η

(1)

本文進(jìn)一步將y看作含有風(fēng)車偽影的CT圖像，x為原始理想CT圖像，η為風(fēng)車偽影。風(fēng)車偽影校正的目的就是盡量去除風(fēng)車偽影η，同時還原出原始理想CT圖像x。

1.1.1 空域雙側(cè)濾波

雙域?yàn)V波對導(dǎo)引圖像g與含干擾圖像y進(jìn)行雙側(cè)濾波[13]。本文中的引導(dǎo)圖像為輸入圖像本身，通過引導(dǎo)圖像對未知的輸出結(jié)果進(jìn)行約束[17]。將雙側(cè)濾波定義在一個以任意像素p為中心2r為邊長的正方形鄰域Np內(nèi)。加權(quán)系數(shù)是以該鄰域內(nèi)圖像內(nèi)容為基礎(chǔ)構(gòu)造的。對于中心像素p，加權(quán)系數(shù)kp,q由兩個加權(quán)因子相乘得到,如公式(4)所示。

(2)

(3)

其中雙側(cè)濾波核函數(shù)為：

(4)

參數(shù)σs和γr用于調(diào)整空間核與像素值核形狀。

1.1.2 STFT

(5)

(6)

(7)

1.1.3 頻域小波收縮

雙域?yàn)V波進(jìn)一步縮小傅里葉噪聲系數(shù)Sp,f，使用類似于公式(4)中的雙側(cè)濾波器像素值核的收縮因子：

(8)

(9)

其中收縮因子Kp,f的構(gòu)造利用了導(dǎo)引圖像譜Gp,f，而小波收縮因子調(diào)整參數(shù)γf的作用類似于雙側(cè)濾波像素值核調(diào)整參數(shù)γr。

所以校正后圖像可近似為兩個去噪層之和。文獻(xiàn)[16]又對上述雙域?yàn)V波進(jìn)一步進(jìn)行改進(jìn)，用含噪圖像與雙域?yàn)V波去噪圖像間差異圖像估計噪聲圖像，并迭代地提高估計精度；基于雙側(cè)濾波權(quán)系數(shù)估計非平穩(wěn)噪聲方差；改進(jìn)雙域?yàn)V波算法，使之更適應(yīng)實(shí)際圖像去噪應(yīng)用。

1.2 距離變換

在上述改進(jìn)雙域?yàn)V波的基礎(chǔ)上，本文提出了基于雙域?yàn)V波與距離變換的CT圖像風(fēng)車偽影校正方法。雙側(cè)濾波核定義在一個以任意像素p為中心2r為邊長的正方形鄰域Np內(nèi)，而實(shí)際螺旋CT圖像風(fēng)車偽影是隨空間位置而變化的。對高對比度物體周圍進(jìn)行距離變換，并隨著距離增大而不斷增大r，改進(jìn)雙域?yàn)V波算法，使其更適應(yīng)實(shí)際螺旋CT圖像風(fēng)車偽影校正：

(10)

式中，r0為初始雙側(cè)濾波核半徑的大??；bwdist為根據(jù)密度變化較大的高對比度物體而進(jìn)行的距離變換；k、s和指數(shù)n用來調(diào)節(jié)r。由于風(fēng)車偽影主要分布在高對比度物體附近，且其強(qiáng)度與頻率隨空間距離變化，通過對高低密度植入物進(jìn)行簡單閾值分割，再利用Matlab中bwdist函數(shù)進(jìn)行空間距離變換。

預(yù)先求出幾組不同r的風(fēng)車偽影校正圖像，然后利用線性插值，使得雙域?yàn)V波更適用于螺旋CT圖像風(fēng)車偽影校正。

2 仿真模擬

分別對高密度物體中各植入物進(jìn)行編號，如圖1所示。植入物參數(shù)設(shè)置如表1所示。其中材料1～5由線性插值得到，植入物7～13為同一種材質(zhì)。將低密度物體植入物均設(shè)置為真空，如圖5所示，植入物參數(shù)設(shè)置不再贅述。

表1 高密度物體植入物參數(shù)設(shè)置

圖1 高密度物體植入物編號(z=0)

通過開發(fā)高低密度物體的數(shù)學(xué)體模，利用計算機(jī)仿真生成圓軌道掃描CT投影數(shù)據(jù)重建得到無風(fēng)車偽影的圖像(如圖2與圖4所示)作為理想重建圖像，生成螺旋掃描CT投影數(shù)據(jù)重建得到含有風(fēng)車偽影的圖像(如圖3與圖5所示)，用來進(jìn)行風(fēng)車偽影校正性能研究。

圖2 高密度物體圓軌道掃描理想重建圖像

圖3 高密度物體周圍含有風(fēng)車偽影的重建圖像

圖4 低密度物體圓軌道掃描理想重建圖像

CT圖像重建采用了廣泛應(yīng)用的螺旋FDK重建算法。CT投影模擬與FDK圖像重建各參數(shù)設(shè)置完全依據(jù)了合作單位北京航星制造有限公司所生產(chǎn)的真實(shí)X射線安檢CT機(jī)參數(shù)設(shè)計的，具體參數(shù)如表2所示。

表2 CT投影仿真與FDK圖像重建參數(shù)

實(shí)驗(yàn)中同時采用原始雙域?yàn)V波[16](即未經(jīng)距離變換插值)對風(fēng)車偽影進(jìn)行校正，用于與本文所提方法(即基于雙域?yàn)V波與距離變換)進(jìn)行對比。k=1.3，s=3，n=1.2，r0=2，r=2、3、5、7、11、13、15，z軸切片位置=體模半徑×2z/z軸總切片數(shù) mm=90×2z/45 mm=4zmm，其中z表示切片序號[-22、-21、…、-1、0、1、…、21、22]。圖1～圖5灰度顯示范圍均為[0.15,0.3]。

圖5 低密度物體周圍含有風(fēng)車偽影的重建圖像

3 結(jié)果

本節(jié)對實(shí)驗(yàn)結(jié)果進(jìn)行比較，其中CT圖像灰度顯示范圍均為[0.185,0.215]，差值圖為校正前減去校正后的差異值，其顯示范圍均為[-0.01,0.01]。采用原始雙域?yàn)V波[15]對風(fēng)車偽影校正，結(jié)果如圖6所示，其中r=3。在圖6實(shí)驗(yàn)基礎(chǔ)上，對雙側(cè)濾波核半徑r進(jìn)行調(diào)整，當(dāng)r=13時實(shí)驗(yàn)結(jié)果如圖7所示。

圖6 當(dāng)r=3時，基于原始雙域?yàn)V波對高密度物體周圍風(fēng)車偽影校正

圖7 當(dāng)r=13時，基于原始雙域?yàn)V波對高密度物體周圍風(fēng)車偽影校正

由圖6明顯可見，原始雙域?yàn)V波校正[16]當(dāng)雙側(cè)濾波核的半徑r較小時，校正結(jié)果會存在風(fēng)車偽影殘留，校正效果較差；當(dāng)半徑r較大時，高密度物體間會產(chǎn)生次生影響，如圖7箭頭所指區(qū)域，從而導(dǎo)致校正性能下降，說明直接使用雙域?yàn)V波進(jìn)行風(fēng)車偽影校正具有一定的局限性?；谏鲜鰧?shí)驗(yàn)結(jié)果，本文提出了基于雙域?yàn)V波和距離變換的螺旋CT圖像風(fēng)車偽影校正，結(jié)果如圖8所示。

由圖8可見，基于雙域?yàn)V波與距離變換進(jìn)行風(fēng)車偽影校正效果要優(yōu)于原始雙域?yàn)V波，由圖9(z=-4)可見，校正后螺旋CT圖像曲線與理想圓軌道CT圖像曲線更為相似；但中心切片處(圖9(z=0))基本不存在風(fēng)車偽影，而未校正螺旋CT圖像與理想圓軌道CT圖像中同時存在一些其他原因?qū)е碌臈l紋偽影，本文方法對其也有一定的去除效果。但需注意，如圖8紅色圓圈內(nèi)所示，較小植入物經(jīng)校正后消失了，表明本文方法對CT圖像空間分辨率有一定負(fù)面影響；同時,較大植入物經(jīng)校正后灰度發(fā)生變化，表明本文方法對CT圖像灰度對比度也有一定負(fù)面影響。

圖8 基于雙域?yàn)V波和距離變換對高密度周圍風(fēng)車偽影校正

圖9 基于雙域?yàn)V波和距離變換對高密度物體周圍風(fēng)車偽影校正前后結(jié)果曲線

上述實(shí)驗(yàn)是對高密度物體周圍進(jìn)行風(fēng)車偽影校正，接下來嘗試對低密度物體周圍進(jìn)行風(fēng)車偽影校正。首先利用原始雙域?yàn)V波[16]進(jìn)行風(fēng)車偽影校正，當(dāng)r=3時校正結(jié)果如圖10所示，當(dāng)r=13時校正結(jié)果如圖11所示。由圖10與圖11可見，當(dāng)r較小時校正結(jié)果仍會存在風(fēng)車偽影殘留，當(dāng)r較大時校正結(jié)果有一定的負(fù)面影響，如圖11中箭頭所指區(qū)域。接下來利用本文基于雙域?yàn)V波和距離變換進(jìn)行風(fēng)車偽影校正，實(shí)驗(yàn)結(jié)果如圖12所示。

圖10 r=3時，基于原始雙域?yàn)V波對低密度周圍風(fēng)車偽影校正

圖11 r=13時，利用原始雙域?yàn)V波對低密度周圍風(fēng)車偽影校正

在上述原始雙域?yàn)V波對低密度物體周圍風(fēng)車偽影校正的基礎(chǔ)上增加距離變換插值，根據(jù)距離變換結(jié)果來改變雙側(cè)濾波半徑r大小，由圖12可見，基于雙域?yàn)V波與距離變換能有效地減少風(fēng)車偽影，由圖13(z=-4)可見，校正后的螺旋CT圖像曲線與理想圓軌道CT圖像曲線更為相似，而中心切處[圖13(z=0)]圖像曲線的特殊性前文已清楚說明不再贅述。同樣需注意，如圖12紅色圓圈內(nèi)所示，本文方法對CT圖像空間分辨率與灰度對比度均有一定負(fù)面影響。

圖12 基于雙域?yàn)V波和距離變換對低密度周圍風(fēng)車偽影校正

圖13 基于雙域?yàn)V波和距離變換對低密度物體的風(fēng)車偽影校正前后結(jié)果曲線

4 結(jié)論

本文利用雙域?yàn)V波進(jìn)行螺旋CT圖像風(fēng)車偽影校正，將雙側(cè)濾波定義在一個以任意像素為中心特定半徑的正方形鄰域內(nèi)。當(dāng)鄰域半徑較小時對風(fēng)車偽影的校正能力較弱; 當(dāng)鄰域半徑不斷增大時，風(fēng)車偽影校正能力不斷增強(qiáng)。但是，隨著鄰域半徑的不斷增大，校正后的CT圖像會丟失細(xì)節(jié)特征。為此，將雙域?yàn)V波與距離變換相結(jié)合，根據(jù)離高低密度物體的距離來調(diào)節(jié)雙側(cè)濾波的鄰域半徑，從而提升了雙域?yàn)V波對螺旋CT圖像風(fēng)車偽影的校正能力，但需注意，本文方法對CT圖像空間分辨率與灰度對比度均有一定負(fù)面影響。為了更好地保持原始CT圖像的細(xì)節(jié)特征，未來將改進(jìn)本文方法，進(jìn)一步提高風(fēng)車偽影校正性能，并對實(shí)際螺旋CT圖像風(fēng)車偽影校正效果進(jìn)行驗(yàn)證。