何世宇

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

0 引言

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

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

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

1 基礎(chǔ)理論

1.1 Radon變換

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

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

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

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

其離散近似可寫為：

圖1 物體在L的線積分

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

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

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

1.2 濾波反投影

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

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

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

圖3 FBP算法流程圖

由此我們可以得到物體對射線能量吸收強度的表達式：

2 MATLAB函數(shù)

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

2.1 radon函數(shù)

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

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

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

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

圖4 Shepp-Logan測試模型

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

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

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

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

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

2.2 iradon函數(shù)

iradon函數(shù)是MATLAB中提供的CT復(fù)原函數(shù)，其原理是使用濾波反變換對切片圖像進行復(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ù)時選用的內(nèi)插函數(shù)，默認使用“l(fā)inear”線性插值；filter指定了濾波反投影時選擇的濾波器，默認為“Ram-lak”濾波器，它可以恢復(fù)圖像的頻率成分，頻率越高，恢復(fù)越多。frequency_scaling是處于（0，1）范圍內(nèi)的標量，通過改變頻率軸的比例來修改濾波器，默認1，若小于1濾波器將被壓縮。output_size是標量，規(guī)定了重建圖像的行數(shù)和列數(shù)。默認值為：

3 實例分析

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

3.1 圖像轉(zhuǎn)動問題的解決

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

圖7 模板及其Radon變換

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

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

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

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

3.2 圖像缺失問題的解決

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

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

4 結(jié)語

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

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

圖11 擴展范圍后的濾波反投影圖像