基于改進(jìn)卷積反投影模型的CT成像偽影消除方法*

韓寅奔， 董泉潤(rùn)

(西北工業(yè)大學(xué) a. 教務(wù)處， b. 動(dòng)力與能源學(xué)院， 西安 710072)

隨著現(xiàn)代工業(yè)的發(fā)展，計(jì)算機(jī)斷層成像(computed tomography，CT)技術(shù)被越來(lái)越多地運(yùn)用到無(wú)損檢測(cè)和逆向工程等行業(yè)中[1].CT重建是一種無(wú)損檢測(cè)技術(shù)，其利用精確準(zhǔn)直的X射線(xiàn)在不破壞被測(cè)物體的前提下，根據(jù)被測(cè)物體不同組織結(jié)構(gòu)吸收射線(xiàn)能力的差異性對(duì)被測(cè)物體進(jìn)行斷層成像，能省時(shí)、無(wú)損地探測(cè)物體內(nèi)部情況，獲取內(nèi)部結(jié)構(gòu)信息[2-4].

目前，常用的CT重建算法主要有迭代法和解析法兩種[5].迭代法能處理不同方式的采樣數(shù)據(jù)，且能獲得稀疏采樣數(shù)據(jù)的重建圖像，但該種方法需要較大的計(jì)算量，重建速度較慢[6-7]；解析法中最具代表性的是卷積反投影法，其通過(guò)處理完全等間距的投影數(shù)據(jù)來(lái)獲得更好的成像質(zhì)量和更快的成像速度，但這種方法得到的重建圖像容易產(chǎn)生偽影[8-9].

相比于迭代法，卷積反投影法的應(yīng)用范圍更廣.許多專(zhuān)家和學(xué)者提出了諸多方法來(lái)提高卷積反投影法的重建速度與成像質(zhì)量，消除偽影.在提升重建速度方面，文獻(xiàn)[10]基于極坐標(biāo)的特性，提出了一種PCBP重建算法，但其在坐標(biāo)轉(zhuǎn)換過(guò)程中需要較大的計(jì)算量；文獻(xiàn)[11]基于正余弦函數(shù)的特性，同時(shí)重建16幅采樣圖像，明顯提升重建速度；文獻(xiàn)[12]則基于三角函數(shù)的對(duì)稱(chēng)性，使用對(duì)稱(chēng)優(yōu)化算法同時(shí)重建64幅圖像.而在提升成像質(zhì)量方面，文獻(xiàn)[13]通過(guò)分析不同偽影的種類(lèi)和成因，分別消除各種醫(yī)學(xué)圖像的偽影；文獻(xiàn)[14]則基于周期性原理，使用正切函數(shù)來(lái)優(yōu)化重建圖像，但這種方法得到的重建圖像仍存在部分偽影.

針對(duì)上述問(wèn)題，本文提出了一種基于卷積反投影算法的CT成像偽影消除方法.首先分析并標(biāo)定給定模板的測(cè)量參數(shù)；然后使用Radon和傅里葉變換得到反投影方程，并計(jì)算每一點(diǎn)吸收率與投影值的積分關(guān)系，從而得到較清晰的重建圖像.

1 參數(shù)標(biāo)定

使用CT系統(tǒng)進(jìn)行無(wú)損檢測(cè)時(shí)，保持射線(xiàn)源和探測(cè)器不動(dòng).CT系統(tǒng)探測(cè)模板與著色圖及其坐標(biāo)系示意圖如圖1所示.為采集不同旋轉(zhuǎn)角度的投影信息，探測(cè)物包含橢圓和圓形兩個(gè)模板.然而，理論上旋轉(zhuǎn)中心與射線(xiàn)源焦點(diǎn)的連線(xiàn)是垂直于探測(cè)器的，但是在安裝CT系統(tǒng)時(shí)，難以實(shí)現(xiàn)這一要求，通常存在誤差，并導(dǎo)致重建圖像存在偽影.因此，需要對(duì)安裝好的CT系統(tǒng)進(jìn)行參數(shù)標(biāo)定，即借助模板標(biāo)定CT系統(tǒng)的參數(shù)，并使用這些參數(shù)對(duì)未知結(jié)構(gòu)的樣品進(jìn)行成像.

本文首先繪制模板投影信號(hào)著色圖，得到模板旋轉(zhuǎn)軌跡.圖1b中用黃色表示投影信號(hào)強(qiáng)度為零，并用深淺不同的藍(lán)色表示不同強(qiáng)度的投影信號(hào)，強(qiáng)度越高，顏色越深；然后以著色圖左下角為原點(diǎn)，橫軸表示旋轉(zhuǎn)角度，縱軸表示其探測(cè)器端點(diǎn)的距離，建立直角坐標(biāo)系；最后結(jié)合著色圖的幾何特征，建立其與實(shí)際照射情況的對(duì)應(yīng)關(guān)系，并在此基礎(chǔ)上標(biāo)定CT系統(tǒng)的旋轉(zhuǎn)中心、探測(cè)器單元間距離及探測(cè)方向等參數(shù).

圖1 探測(cè)模板與著色圖及其坐標(biāo)系Fig.1 Detection template，coloring diagram and coordinate system

1.1 探測(cè)器單元間距離

當(dāng)X射線(xiàn)平行于橢圓長(zhǎng)軸入射時(shí)，光在介質(zhì)中走過(guò)的路徑最長(zhǎng)，衰減幅度最大，接收到的信號(hào)最強(qiáng).此時(shí)，橢圓模板在接收器上的投影長(zhǎng)度最短且恰為短軸長(zhǎng)，即圖1b中A、B兩點(diǎn)的縱坐標(biāo)之差.測(cè)得AB之間共有L個(gè)單元格，即對(duì)應(yīng)L-1個(gè)探測(cè)單元間距，故可得探測(cè)器單元間距離d=m/(L-1)，其中，m為橢圓模板的短軸長(zhǎng).

1.2 系統(tǒng)照射方向

模板每旋轉(zhuǎn)360°后將恢復(fù)為初始狀態(tài)，著色圖上的曲線(xiàn)則完成一個(gè)首尾相接的周期.由模板形狀可知，其旋轉(zhuǎn)0°～180°所得數(shù)據(jù)與旋轉(zhuǎn)180°～360°所得數(shù)據(jù)相同，因此，只需要將探測(cè)物轉(zhuǎn)過(guò)180°即可獲得實(shí)驗(yàn)所需的所有數(shù)據(jù).

設(shè)X射線(xiàn)與旋轉(zhuǎn)中心和模板上小圓圓心連線(xiàn)之間所夾銳角為φ，則小圓圓心在探測(cè)器上投影點(diǎn)距起始點(diǎn)的距離l可表示為

l=l0+d0sinφ

(1)

式中：l0為旋轉(zhuǎn)中心投影點(diǎn)到探測(cè)器起始點(diǎn)的距離；d0sinφ為旋轉(zhuǎn)中心與圓心連線(xiàn)在探測(cè)器上的投影長(zhǎng)度.

由式(1)可知，模板與探測(cè)器始端的距離可以表示為三角函數(shù)形式.而式(1)則為著色圖中等寬曲線(xiàn)的數(shù)學(xué)描述，故可使用傅里葉級(jí)數(shù)擬合該曲線(xiàn)并求得其跨度q.由CT系統(tǒng)的工作原理可知，系統(tǒng)每次轉(zhuǎn)過(guò)的角度相同，故探測(cè)器每次轉(zhuǎn)過(guò)180°/q.

1.3 旋轉(zhuǎn)中心標(biāo)定

如圖1b所示，C點(diǎn)對(duì)應(yīng)等寬曲線(xiàn)的最高點(diǎn)，設(shè)θ1、θ2分別為A、C兩點(diǎn)在坐標(biāo)圖上對(duì)應(yīng)的橫坐標(biāo)，A、C兩點(diǎn)的橫軸間距θm為

θm=θ1-θ2

(2)

在實(shí)際情況下，θm與光線(xiàn)轉(zhuǎn)過(guò)夾角相等，對(duì)應(yīng)狀態(tài)幾何關(guān)系如圖2所示.

圖2 A、C兩點(diǎn)對(duì)應(yīng)狀態(tài)的幾何關(guān)系Fig.2 Geometric relation of corresponding state between A and C points

由于探測(cè)器每次轉(zhuǎn)過(guò)180°/q，故通過(guò)觀(guān)察A、C兩點(diǎn)之間對(duì)應(yīng)列的列數(shù)I，即可求得

θm=I180°/q

(3)

當(dāng)?shù)葘捛€(xiàn)位于C點(diǎn)時(shí)，模板圓心與旋轉(zhuǎn)中心的間距被完全投影到探測(cè)器上，表明此時(shí)圓心投影點(diǎn)距探測(cè)器始端最遠(yuǎn)；而當(dāng)位于其對(duì)稱(chēng)位置時(shí)，距探測(cè)器始端最近.故只需計(jì)算出圖1b中等寬曲線(xiàn)兩極值點(diǎn)豎直方向上相差的單元格數(shù)目T，即可求出圓心與旋轉(zhuǎn)中心之間的距離D，即

l(T-1)=2D

(4)

若以模板在圖2的幾何狀態(tài)為基準(zhǔn)，以左下角為原點(diǎn)建立直角坐標(biāo)系，則旋轉(zhuǎn)中心的坐標(biāo)值表示為

(5)

2 模型CT成像偽影消除

CT圖像斷層平面中某一點(diǎn)的密度值可以看作平面內(nèi)所有經(jīng)過(guò)該點(diǎn)的射線(xiàn)投影值之和的均值.然而觀(guān)測(cè)樣品各點(diǎn)的形貌與吸收率各不相同，有些點(diǎn)比較突出，經(jīng)卷積反投影處理后，這些點(diǎn)依舊突出，并可據(jù)此還原斷層；但有些原本不為零的地方現(xiàn)在成為了零，即產(chǎn)生偽跡，因此，需要在反投影之前對(duì)投影進(jìn)行濾波處理.

2.1 Radon變換及傅里葉變換

Pφ(ρ)=P(ρ，φ)

(6)

式中：w1、w2為傅里葉變換[15]的特征參量；F為傅里葉變換映射；P為投影函數(shù)p的傅里葉變換.可由傅里葉反變換求得待重建圖像，即

a(x，y)=F-1[E(w1，w2)]=

(7)

式中：α為實(shí)際某點(diǎn)投影角度；r為弧長(zhǎng).對(duì)于式(7)進(jìn)行積分，可以整理為

g[rcos(α-φ)，φ]

(8)

故待建圖像轉(zhuǎn)化為

(9)

由前述推導(dǎo)過(guò)程可得

(10)

反投影得到托盤(pán)上各點(diǎn)的吸收率為

(11)

考慮到實(shí)際分析處理問(wèn)題時(shí)需要將介質(zhì)劃分為256×256的微元格點(diǎn)研究，故需對(duì)上述連續(xù)模型進(jìn)行離散化處理，則有

(12)

濾波因子采樣得到

(13)

其取值為

(14)

可將離散條件下的反投影積分轉(zhuǎn)化為級(jí)數(shù)求和，即

(15)

式中，N為一個(gè)微元滿(mǎn)足的最小旋轉(zhuǎn)角的最大倍數(shù).將式(15)離散化則有

(16)

式中，Δx、Δy為橫縱坐標(biāo)最小變化量.

2.2 偽影消除

使用2.1節(jié)變換操作可確定絕大部分點(diǎn)的介質(zhì)吸收率，然而，實(shí)際操作中發(fā)現(xiàn)有少部分點(diǎn)并沒(méi)有X射線(xiàn)穿過(guò)而無(wú)法確定其吸收率[16].若使用近鄰插值方法求取這些點(diǎn)的吸收率，雖仍能保持樣品的形貌特征，但重建圖像部分邊緣較模糊，并形成偽跡.因此為抑制偽跡出現(xiàn)，一方面要對(duì)源數(shù)據(jù)進(jìn)行降噪處理；另一方面需要對(duì)變換后的數(shù)據(jù)進(jìn)行插值賦值，減小數(shù)據(jù)的偏離幅度.

本文使用3階樣條法對(duì)式(16)得到的重建圖像進(jìn)行插值處理.3階樣條法除了要求在兩端與給定數(shù)據(jù)值相等外，還要求樣條系統(tǒng)的光滑性，即內(nèi)部鄰接數(shù)據(jù)點(diǎn)的1階及2階導(dǎo)數(shù)必須匹配.在外側(cè)端點(diǎn)1階導(dǎo)數(shù)沒(méi)有變化，2階導(dǎo)數(shù)為零.樣條函數(shù)定義為

(17)

其導(dǎo)數(shù)為

(18)

為滿(mǎn)足數(shù)據(jù)匹配條件，則有

(19)

內(nèi)部數(shù)據(jù)點(diǎn)1階、2階導(dǎo)數(shù)匹配，故在兩函數(shù)銜接點(diǎn)gn(i2)處有

(20)

在外側(cè)端點(diǎn)1階導(dǎo)數(shù)沒(méi)有變化，故2階導(dǎo)數(shù)為0，即

(21)

聯(lián)立式(18)～(21)可解得具體的插值函數(shù)，并可將[x1，x3]區(qū)間內(nèi)任一點(diǎn)帶入相應(yīng)坐標(biāo)求得其值.

3 算例分析

本文使用圖1a所示的模板標(biāo)定CT系統(tǒng)的參數(shù)，使用文中提出的方法求得探測(cè)器單元間距離l為0.280 4 mm，旋轉(zhuǎn)中心的坐標(biāo)為(58.633 1 mm，42.395 3 mm)，著色圖中等寬曲線(xiàn)在整個(gè)數(shù)據(jù)點(diǎn)區(qū)間內(nèi)跨度大約為180°，探測(cè)器每次轉(zhuǎn)過(guò)1°，傅里葉變換特征參量w1、w2分別取30和60，濾波因子由式(14)求得.

使用標(biāo)定后的系統(tǒng)參數(shù)進(jìn)行測(cè)試實(shí)驗(yàn)，得到如圖3所示的反投影結(jié)果.從圖3可以看出，得到的重建圖較為模糊，且輪廓周?chē)河羞^(guò)渡白色.使用3階樣條法進(jìn)行插值處理后得到的結(jié)果如圖4所示，處理后的圖像明顯清晰，內(nèi)部分布有兩個(gè)彼此分離的橢圓空洞、兩個(gè)相交的更高吸收率的橢圓形小區(qū)域.

圖3 未經(jīng)去噪、插值處理的圖像Fig.3 Images without denoising and interpolating processes

圖4 降噪、插值后樣品最終圖像Fig.4 Final image of sample after noise reduction and interpolation

本文也比較了消除偽影前后10個(gè)給定位置處的介質(zhì)吸收率，如表1、2所示.根據(jù)表1結(jié)果可以看出，部分吸收率為負(fù)值，但其值均較小，表明實(shí)際吸收率應(yīng)為0，存在一定的測(cè)量誤差.比較表1、2的吸收率可知，消除偽影后能測(cè)得每個(gè)位置的吸收率，并能得到更清晰的重建圖像.

表1 偽影消除前給定10個(gè)位置處吸收率Tab.1 Absorption rates at given 10 positions before artifact removal

表2 偽影消除后給定10個(gè)位置處吸收率Tab.2 Absorption rates at given 10 positions after artifact removal

3.1 算法精度分析

實(shí)驗(yàn)中的誤差將降低測(cè)試精度，這些誤差主要包括截?cái)嗾`差、公式化誤差、舍入誤差和測(cè)量誤差.在參數(shù)標(biāo)定時(shí)，本文討論不涉及多項(xiàng)式近似，且忽略了次要條件的影響，并不會(huì)大幅度影響最終的結(jié)果.因此，本文算法的主要誤差來(lái)源為給定實(shí)驗(yàn)數(shù)據(jù)、計(jì)算軟件的舍入誤差以及數(shù)據(jù)收集過(guò)程中不精確性引入的測(cè)量誤差.下面以圓形遮光區(qū)域?yàn)槔?，推?dǎo)吸光區(qū)尺寸與誤差的定量關(guān)系.

設(shè)k為圓遮擋的X射線(xiàn)數(shù)，設(shè)z為真實(shí)值與測(cè)量值之間的偏差，則有

kd+z=2L

(22)

式中，2L為測(cè)量值，故相對(duì)誤差為

(23)

由于一般情況下誤差會(huì)控制在較小范圍內(nèi)，即kd≈2L，則有

(24)

由式(24)可知，隨著物體半徑L的增大，誤差η在減小，表明增大模板吸光區(qū)半徑，可有效減小誤差.如將半徑擴(kuò)大為原來(lái)的二倍，則誤差為

可見(jiàn)半徑擴(kuò)大為原來(lái)的二倍，誤差幾乎變?yōu)樵瓉?lái)的1/2，因此，對(duì)于小圓半徑相對(duì)于整個(gè)平板較小的情況，參數(shù)的測(cè)量精度還有較大的提升空間.但由于偏差z值較小，即使誤差減半，也不會(huì)對(duì)最終結(jié)果影響過(guò)大，故本文測(cè)得的參數(shù)精度是可以接受的.

3.2 算法穩(wěn)定性分析

本文假設(shè)X射線(xiàn)與橢圓中心、小圓圓心的連線(xiàn)恰好垂直，如圖2所示.但實(shí)際上取樣是非連續(xù)的，只能近似為垂直關(guān)系，由于探測(cè)器間距較近，這種近似導(dǎo)致的誤差雖會(huì)隨具體情況變化，但其值較小，在一定程度上是可接受的.

本文測(cè)試了以旋轉(zhuǎn)中心坐標(biāo)為例的參數(shù)值對(duì)系統(tǒng)參數(shù)的敏感性.將式(5)中的旋轉(zhuǎn)中心橫縱坐標(biāo)對(duì)夾角θm求導(dǎo)，得到

(25)

4 結(jié) 論

本文提出了一種基于卷積反投影算法的CT成像偽影消除方法.利用實(shí)物與投影的對(duì)應(yīng)關(guān)系，根據(jù)投影方便、快捷地求出探測(cè)器間距等待系統(tǒng)參數(shù)；使用卷積反投影技術(shù)和三次樣條插值算法，精確地求出每一點(diǎn)的吸收率并得到較清晰的重建圖像.算例仿真測(cè)試的精確度、穩(wěn)定性分析結(jié)果表明，提出的算法具有簡(jiǎn)明易懂、操作方便和精確度高等諸多優(yōu)點(diǎn)，且其相較于傳統(tǒng)方法具有更高的運(yùn)算效率.