基于Radon變換與濾波反投影的CT模型及應(yīng)用

何世宇

（河海大學(xué)計(jì)算機(jī)與信息學(xué)院，南京 211100）

0 引言

CT（Computed Tomography）技術(shù)，也稱計(jì)算機(jī)斷層成像技術(shù)，是利用X射線對(duì)不同厚度物品穿透性與衰減性的不同，在不破壞待測(cè)品內(nèi)部結(jié)構(gòu)的情況下,根據(jù)物體周邊所獲取的某種物理量（如波速、X線光強(qiáng)、電子束強(qiáng)等）的投影數(shù)據(jù)，運(yùn)用一定的數(shù)學(xué)方法，通過(guò)計(jì)算機(jī)處理，重建物體特定層面上的二維圖像以及依據(jù)一系列上述二維圖像構(gòu)成三維圖像的技術(shù)。

從1971年Hounsfielld發(fā)明頭顱CT到20世紀(jì)80年代，CT技術(shù)的發(fā)展主要在于掃描部位的延伸，即從單一的頭部檢查拓展到體部檢查。CT技術(shù)的發(fā)展大致分為以下四個(gè)階段：第一代CT裝置通過(guò)平行束旋轉(zhuǎn)掃描獲得投影數(shù)據(jù)；第二代使用小角度扇形射線束來(lái)代替平行束；第三代設(shè)備僅包含扇形束的旋轉(zhuǎn)掃描動(dòng)，不再采用平移動(dòng)作；第四代設(shè)備將探測(cè)器固定在360°的圓周上，僅旋轉(zhuǎn)X射線源以解決幻影偽像問(wèn)題。由于被測(cè)物與掃描環(huán)境的復(fù)雜性，CT掃描方式日趨靈活。

由于CT技術(shù)的特點(diǎn)，使其廣泛應(yīng)用于醫(yī)療檢查。除了醫(yī)學(xué)領(lǐng)域，CT在物質(zhì)探測(cè)方面同樣具有的巨大的優(yōu)勢(shì)，尤其在工業(yè)、地球物理、工程、農(nóng)業(yè)、安全檢測(cè)等行業(yè)得到了廣泛的應(yīng)用。

1 基礎(chǔ)理論

1.1 Radon變換

Radon變換是CT技術(shù)的主要理論基礎(chǔ)，1917年由數(shù)學(xué)家Radon證明。Radon變換是一種積分變換，它的本質(zhì)是對(duì)原來(lái)的函數(shù)做空間轉(zhuǎn)換，將XY平面上的點(diǎn)映射到AB平面，原來(lái)XY平面上的某條直線的所有的點(diǎn)在AB平面上都位于同一點(diǎn)，即對(duì)物體進(jìn)行線積分。根據(jù)AB平面上點(diǎn)的量化指標(biāo)，可得XY平面在該條線上的部分特性。這便是Radon變換的實(shí)質(zhì)。

假設(shè)有一條穿過(guò)空間中某物體的X射線，其方程為：

其中θ影響射線的斜率，ρ影響射線的位置。同時(shí)我們假設(shè)空間中該物體對(duì)射線能量的吸收強(qiáng)度可以表示成函數(shù)：f(x,y)，那么物體在該條X射線上所吸收的能量可以表示為積分：

其中ds是該射線的微分。上述積分常借助Delta函數(shù)的抽樣性進(jìn)行化簡(jiǎn)求解，表示為：

其離散近似可寫為：

圖1 物體在L的線積分

由此可見(jiàn)只要給定一組 ρ和θ，即可以得出一個(gè)沿L的積分值。因此，Radon變換就是函數(shù) f(x,y)的線積分。

從一條射線推廣到多條，假設(shè)有多條平行與L的線，它們有相同的θ，不同的 ρ,我們對(duì)每一條這樣的平行線都做 f(x,y)的線積分，會(huì)產(chǎn)生很多投影線，如圖2所示。也就是說(shuō)對(duì)一幅圖像在某一特定角度下的Ra?don變換會(huì)產(chǎn)生N個(gè)線積分值，而每一個(gè)線積分值會(huì)對(duì)應(yīng)一個(gè)自己的位置ρ。各個(gè)角度的Radon變換值匯總在一起就構(gòu)成一幅Radon變化圖。

圖2 物體在射線束下的投影

1.2 濾波反投影

濾波反投影（FBP）算法是在傅立葉中心切片理論基礎(chǔ)上的一種算法。它在反投影前將每一個(gè)采集角度下的投影進(jìn)行卷積處理，從而改善復(fù)原后函數(shù)的模糊現(xiàn)象，重建的圖像質(zhì)量較好。

根據(jù)傅立葉中心切片定理的定義可知，對(duì)投影的一維傅立葉變換等效于對(duì)原圖像進(jìn)行二維的傅立葉變換。有了該理論的支持，就能通過(guò)在投影上執(zhí)行傅立葉變換，得到物體的的二維傅立葉變換。

由此，我們可以列出投影重建的步驟，在投影重建的過(guò)程中，先把線陣探測(cè)器上獲得的投影（衰減后的射線能量）數(shù)據(jù)g(ρ ,θ)進(jìn)行一次一維傅立葉變換得到G(w ,θ)，再將 g(ρ ,θ)與濾波器函數(shù)H(w)進(jìn)行卷積運(yùn)算，得到各個(gè)方向卷積濾波后的投影數(shù)據(jù)g'(ρ ,θ)；然后把它們沿各個(gè)方向進(jìn)行反投影，即按原路徑分配到每一矩陣單元上，進(jìn)行重疊后得到每一矩陣單元的CT值；再經(jīng)過(guò)適當(dāng)處理后得到被掃描物體的斷層圖像。具體算法流程見(jiàn)圖3。

圖3 FBP算法流程圖

由此我們可以得到物體對(duì)射線能量吸收強(qiáng)度的表達(dá)式：

2 MATLAB函數(shù)

在 MATLAB中，有專門的radon（）函數(shù)和 iradon（）函數(shù)來(lái)實(shí)現(xiàn)Radon變換與濾波反變換的算法，下面將對(duì)這兩個(gè)函數(shù)的用法以及應(yīng)用效果進(jìn)行簡(jiǎn)單的介紹。

2.1 radon函數(shù)

radon函數(shù)的常見(jiàn)用法有以下兩種:

（1）R=radon（I,theta）：若 theta為一個(gè)數(shù)值，R 值則是計(jì)算圖像I在theta方向上的Radon變換；若theta為一個(gè)角度范圍，那么R則為一個(gè)矩陣，表示在該范圍內(nèi)的Radon變，其一個(gè)維度代表角度（一般為0到180°），一個(gè)維度代表積分值。

（2）[R,xp]=radon（I,theta）：R 的返回值含義同上，矩陣“xp”中包含了射線的位置信息，即ρ信息。

Shepp-Logan是測(cè)試CT系統(tǒng)的常用圖片，我們使用它來(lái)演示radon函數(shù)的使用效果，圖4為Shepp-Lo?gan測(cè)試模型。

圖4 Shepp-Logan測(cè)試模型

output_size=2*floor（size（R,1）/（2*sqrt（2）））

我們對(duì)上一小節(jié)Radon變換產(chǎn)生的四個(gè)圖像進(jìn)行濾波反投影變換，得到如圖6的四幅圖片，可以看出隨著旋轉(zhuǎn)角度變小——也就是增大投影個(gè)數(shù)，圖像重建越好，越接近實(shí)際效果。

圖6 不同投影個(gè)數(shù)的圖形重建情況

下面我們通過(guò)在radon函數(shù)中設(shè)置不同范圍的theta值，對(duì)圖4進(jìn)行Radon變換，我們將180°分別分為了5、10、18、180個(gè)角度，變換效果見(jiàn)圖5。可以看投射的角度越多，產(chǎn)生的Radon變換越清晰。

圖5 不同投影次數(shù)的Radon變換

2.2 iradon函數(shù)

iradon函數(shù)是MATLAB中提供的CT復(fù)原函數(shù)，其原理是使用濾波反變換對(duì)切片圖像進(jìn)行復(fù)原，用法如下：

I=iradon（R,theta,interp,filter,frequency_scaling,output_size）

其中R為反投影數(shù)據(jù)；theta為描述角度的變量，既可以是包含角度的矢量（給出每次旋轉(zhuǎn)的度數(shù)）也可以提供一組行向量（如 0、1、2…178、179）；interp 是字符串，定義了重建函數(shù)時(shí)選用的內(nèi)插函數(shù)，默認(rèn)使用“l(fā)inear”線性插值；filter指定了濾波反投影時(shí)選擇的濾波器，默認(rèn)為“Ram-lak”濾波器，它可以恢復(fù)圖像的頻率成分，頻率越高，恢復(fù)越多。frequency_scaling是處于（0，1）范圍內(nèi)的標(biāo)量，通過(guò)改變頻率軸的比例來(lái)修改濾波器，默認(rèn)1，若小于1濾波器將被壓縮。output_size是標(biāo)量，規(guī)定了重建圖像的行數(shù)和列數(shù)。默認(rèn)值為：

3 實(shí)例分析

2017年全國(guó)大學(xué)生數(shù)學(xué)建模競(jìng)賽的A題即為CT圖像復(fù)原，在這一節(jié)我們將對(duì)該題目的難點(diǎn)，即圖像轉(zhuǎn)動(dòng)與缺失問(wèn)題進(jìn)行分析，并介紹如何使用MATLAB處理上述兩類問(wèn)題。

3.1 圖像轉(zhuǎn)動(dòng)問(wèn)題的解決

題目中給出了模型的形狀及其吸收率，其中黑色部分吸收率為0，白色部分吸收率為1，我們對(duì)模板進(jìn)行標(biāo)準(zhǔn)的Radon變換，模板及變換的結(jié)果見(jiàn)圖7。此外題目還給出了該模板經(jīng)過(guò)某一般條件下的Radon變換數(shù)據(jù)，如圖8。

圖7 模板及其Radon變換

圖8 某一般情況下的Radon變換

由圖8可以得知，由于紅線位置的圖形寬度最窄，吸收的射線強(qiáng)度最弱。假設(shè)從模板正下方投影為0°，紅線位置即為模板180°方向的射線束的線積分。因此我們得知該對(duì)模板Radon變換的初始角度非0度。經(jīng)過(guò)等比例換算得知，Radon變換的初始角度約為30度。下面我們對(duì)原數(shù)據(jù)進(jìn)行濾波反投影變換（見(jiàn)圖9），左圖為角度范圍0°到180°，右圖將角度進(jìn)行修正，角度范圍為 30°到 210°。

圖9 不同角度范圍的濾波反投影圖像

可以發(fā)現(xiàn)右圖相比作圖位置得到更正，位置由未修正前的傾斜恢復(fù)到了模板形狀，但是圖像信息發(fā)生了損失，及丟失掉了右方圓的信息，下一小節(jié)將會(huì)討論如何恢復(fù)未顯示的信息。

3.2 圖像缺失問(wèn)題的解決

為了確定是數(shù)據(jù)丟失還是顯示范圍過(guò)小，我們對(duì)濾波反變換后的圖像進(jìn)行三維顯示（圖10），我們發(fā)現(xiàn)濾波反變換后的數(shù)據(jù)信息中已經(jīng)包含了小圓的信息。

經(jīng)過(guò)對(duì)模板與復(fù)原后圖像的形狀等參數(shù)的分析，我們發(fā)現(xiàn)由于模板的幾何中心與旋轉(zhuǎn)中心存在一定的位置偏移，使復(fù)原后的幾何中心偏離范圍。因此復(fù)原后的圖片信息不能在模板范圍內(nèi)顯示完全，需要比模板更大的范圍才能顯示完備的信息。經(jīng)過(guò)對(duì)iradon函數(shù)的研究發(fā)現(xiàn)，其最后一個(gè)輸入可選項(xiàng)output_size即為可控制圖像大小的參量，我們通過(guò)對(duì)圖形大小的擴(kuò)展，即可復(fù)原出原圖片顯示丟失的內(nèi)容。經(jīng)過(guò)擴(kuò)展后的圖像見(jiàn)圖11。

4 結(jié)語(yǔ)

計(jì)算機(jī)斷層成像作為醫(yī)學(xué)領(lǐng)域重要的技術(shù)之一具有很高的研究?jī)r(jià)值。若僅會(huì)使用MATLAB中的工具箱，對(duì)旋轉(zhuǎn)角度與尺度變換等問(wèn)題需要配置一些額外的參數(shù)才看可以很好的解決，因此對(duì)Radon變換與濾波反投影算法深入的學(xué)習(xí)是必不可少的。此外，旋轉(zhuǎn)中心與幾何中心位置的確定等問(wèn)題，目前還需要針對(duì)具體問(wèn)題進(jìn)行逐一建模才可以解決，如何得到普遍性的模型還有待進(jìn)一步的研究。

圖10 濾波反變換結(jié)果的三維化

圖11 擴(kuò)展范圍后的濾波反投影圖像