基于CT技術的投影重建模型研究*

張 霞 陳佳音 高英喬 董家旋 劉 明

(天津工業(yè)大學數(shù)學科學學院，天津 300387)

0 引 言

計算機斷層成像技術，也就是CT技術作為一種非接觸式檢測技術，能夠在不影響樣品的情況下，獲取其內部結構組織的相關信息，在醫(yī)療檢測、地質勘探領域有廣泛應用．然而，CT系統(tǒng)在安裝的過程中往往存在誤差，這種誤差會影響成像質量，進而影響測量分析結果．

本文借助已知結構的樣品(模板)對安裝好的CT系統(tǒng)進行參數(shù)標定，能有效提高成像結果的精度和穩(wěn)定性．圖1為一種典型的二維CT系統(tǒng)，發(fā)射端n束平行射線垂直于探測器平面入射對應接收點，圖2為模板示意圖，圖3為射線強度衰減規(guī)律示意圖．CT標定問題就是通過探測器采集到的信息，確定CT系統(tǒng)的旋轉中心、探測器間距以及系統(tǒng)每次轉動的角度，進而確定介質各點的吸收率，分析介質的組織結構．

圖1 CT樣機工作原理

圖2 模板示意圖(單位/mm)

圖3 射線強度衰減規(guī)律示意圖

1 參數(shù)的標定

1.1 原理和定義

當射線穿透介質時，攜帶的能量被介質吸收而衰減，故其衰減量與介質的空間分布有關．對于介質均勻且光程極短的厚度范圍內，衰減量與入射強度、穿透厚度成正比，見圖3，具體關系可描述為：

ΔI=-μI0Δl.

(1)

假定Il為射線穿過介質后的強度，I0為原始強度，l為吸收介質的厚度，μ為本文所述的吸收率，其大小取決于介質種類和射線能量且遵循Lambert-Beer定律．對于非均勻復合介質，μ還與介質分布有關，在計算時常分段處理．比如，對于含有n個同等厚度的Δl,有：

I(l)=I0e-Δlμ(l0)e-Δlμ(l1)…e-Δlμ(ln).

(2)

化簡后可得到：

(3)

(4)

對所有吸收量進行如上處理轉換為射線在介質中的光程長，以便于模型后續(xù)的建立和參數(shù)求解．為了確定CT系統(tǒng)旋轉中心的位置、探測器單元間距，根據(jù)模板的幾何特征選取了射線可能經過的3個特殊方向(圖4)，并會出扎3個方向下n條射線的吸收曲線(圖5)．

圖4 X射線發(fā)射端的3個特殊方向(a)方向1; (b)方向2; (c)方向3

圖5 3個特殊方向射線的吸收曲線(a)方向1; (b)方向2; (c)方向3

由于方向1和方向2中的最大吸收分別出現(xiàn)橢圓長軸和短軸上，定義產生最大吸收量的發(fā)射器Si1和Si2，即可確定這兩個發(fā)射器的位置．方向3吸收曲線中的兩個峰相交于橫軸上一點(Si3, 0)，由此可確定發(fā)射器Si3的位置．為了提高射線利用率、降低安裝難度，一般選擇整個系統(tǒng)的幾何中心作為旋轉中心，使得有一條射線(多為第n/2條射線)始終經過旋轉中心．

1.2 探測器單元間距和旋轉中心的確定

建立平面直角坐標系如圖6，用于求解和描述探測器間的距離．當射線方向為方向1時，選取過橢圓長軸和其附近與橢圓相交的若干組射線，其中的一些相關參數(shù)，滿足下列關系：

li=2xi,

(5)

(6)

圖6 用于求解探測器間距的平面直角坐標系

則對應第i個探測器和第i+1個探測器之間的距離為：

di=|yi+1-yi|.

(7)

建立以正方形托盤的左下角為原點、以兩條臨邊為橫軸和縱軸的如下平面直角坐標系(圖7)，用于求解和描述旋轉中心：

圖7 用于求解旋轉中心的平面直角坐標系

通過查閱文獻可知, 第n個探測器到第m個探測器之間的距離為：

(8)

所以，在方向1下，第223條射線在正方形托盤中的位置為橢圓的長軸長，與第256.5條等效射線的間距為：

(9)

通過計算可得，方向1下的第256.5條射線在所建立的平面直角坐標系中的位置為：

x=40.910 3.

(10)

同理，在方向2下，第235條射線在正方形托盤中的位置為橢圓的長軸長，與第256.5條等效射線的間距為：

(11)

通過計算可得，方向1下的第256.5條射線在所建立的平面直角坐標系中的位置為：

y=55.833 6.

(12)

綜上求得CT系統(tǒng)旋轉中心在正方形托盤中的位置為(40.910 3，55.833 6)．

利用方向3對應的位置和數(shù)據(jù)，發(fā)現(xiàn)方向3下的第256.5條射線經過點(40.914 6，55.835 7)，與上文中求得的中心位置十分接近，誤差與探測器間距等其他參數(shù)的線度相比可以忽略，故此模型可靠．

在上述結果的基礎上，求解該CT系統(tǒng)使用的X射線的180個方向的難度大大降低．選取第256.5條等效射線，滿足幾何關系如下：

yi-y0=tanθ(xi-x0),i=1,2,

(13)

(14)

l2=(x1-x2)2+(y1-y2)2.

(15)

其中，l為經過橢圓形介質的光程長．但由于所給數(shù)據(jù)是總光程長，且無法直接去除小圓部分對數(shù)據(jù)的影響，所以使用了基于密度的異常值檢驗法，鎖定了受小圓影響的數(shù)據(jù)區(qū)域，經3次樣條插值修正后，使數(shù)據(jù)符合上述模型的使用條件．

將修正后的數(shù)據(jù)帶入上述模型(圖8)，即可得到180個θ的取值．其中，起始角度為29.604 9°，每個角度的變化量為1(±5%)．

圖8 射線方向滿足的幾何關系

2 投影重建模型的建立與實例分析

2.1 線性衰減重建模型

參考了Radon變換的原理，在其現(xiàn)有的基礎上進行改造，使穿過每個單元的射線方向連續(xù)可調，最后建立了適用于本題的反投影重建模型．

R2為一二維平面，設函數(shù)f(x,y)為該平面上定義的重建圖像密度函數(shù)(正比于吸收率)．直線L為平面內某一條直線，Radon變換算子R即為密度函數(shù)沿L的線積分：

(16)

也可以表示為：

p(s,θ)=Rf(x,y)

(17)

其中，p(s,θ)為f(x,y)的投影函數(shù)，s為原點到直線L的距離，δ為狄拉克函數(shù)，其中s滿足：

s=xcosθ+ysinθ.

(18)

由上述原理可知，CT成像的基本原理為：已知投影數(shù)據(jù)p(s,θ)，計算重構圖像的密度函數(shù)f(x,y)，所以Radon逆變換就是CT成像的過程．由以上Radon變換可得到逆變換為：

f(x,y)

(19)

由投影數(shù)據(jù)計算吸收率μ(x,y)分布最直接的方法就是求解方程組．

參考了這種原理，在其現(xiàn)有的基礎上進行改造，使穿過每個單元的射線方向連續(xù)可調，最后建立了線性衰減重建模型．由此可以用推廣的到n×n像素的圖像，求解此方程即可得到吸收率的二維分布．設Xjk為在(j，k)所在區(qū)域內單位長度介質的吸收率，Pijk為該射線通過該區(qū)域吸收能量的值，該區(qū)域是否被射線穿過用rjk來表示．

(20)

像素Xjk吸收的能量為：

Pijk=rjkXjkdjk.

(21)

該射線的總能量吸收為：

(22)

矩陣表示為P=XRD,式中：

P=[p1,p2,…pn]T.

(23)

吸收率矩陣、判斷該射線是否穿過某區(qū)域的判斷矩陣、射線穿過該區(qū)域的光程矩陣的表達方式與第二問完全相同．選取了一組典型吸收數(shù)據(jù)，利用MATLAB編程，得到該未知介質在正方形托盤中的位置和幾何形狀如圖9：

圖9 線性衰減重建模型重構圖像

2.2 二維濾波反投影模型的建立及求解

參考二維傅里葉重建算法，發(fā)現(xiàn)二維傅里葉重建算法存在不適用于局部目標重建、現(xiàn)有頻域空間插值算法等缺陷，所以針對這些缺陷對二維傅里葉重建算法進行了改進，形成了二維濾波反投影模型．二維傅里葉重建算法的基本原理是：二維圖像一維投影的傅立葉變換精準地等于通過該圖像傅立葉變換中心的直線．當投影角度變化時，過傅立葉變換中心的那條直線隨之變化．具體證明過程與問題二中完全相同．針對這些缺陷，對其進行了改進，采用濾波反投影來重建圖像：先把檢測器上獲得的原始數(shù)據(jù)與一個濾波函數(shù)進行卷積運算，得到個方向上的投影函數(shù)，反投影后經適當處理就能得帶斷層圖像，并保證重構圖像邊緣清晰且內部均勻，在二維傅里葉逆變換的基礎上，利用對稱關系：

P(ω,θ+π)=P(-ω,θ)

(24)

將f(x,y)化為：

(25)

最后得出：

(26)

其中，|ω|表示濾波函數(shù)Haming窗．沿用線性衰減重建模型中實例的吸收數(shù)據(jù)帶入此模型，得到該未知介質在正方形托盤中的位置和幾何形狀如圖10：

圖10 求解方程組重構圖像

3 標定模板的改進

由于標定過程中精度和穩(wěn)定性會直接影響標定結果，在經典模板的基礎上自行設計新模板，建立對應的標定模型，以提高標定精度和穩(wěn)定性．

(1)精確確定探測單元之間的距離.在多次仿真中發(fā)現(xiàn)有某些探測器發(fā)出的射線始終不通過已知結構的樣品(模板)，所以難以精確確定這些探測單元之間的距離，并且在測定新模板時，由于不知道這些探測器單元的精確距離以至于難以保證還原圖像的質量．為了更精準地確定探測單元間距，在原有模板的基礎上增加了兩個圓柱體，位置如圖11(a)所示，通過增加這兩個圓柱體，使所有探測器發(fā)出的射線在不同角度均有機會打到模板上，從而可以精確確定探測單元之間的距離．

圖11 第1～5次改進模板

(2)確定CT系統(tǒng)的旋轉中心在正方形托盤中位置的確定.在上一小問的基礎上改進了模板，在保證所有探測器單元都有機會打到模板的基礎上更改模板圖形的位置，增加更多的特殊位置.如相切，通過更多的特殊位置，可以更加精確的確定CT系統(tǒng)的旋轉中心在正方形托盤中的位置，新建立的模型如圖11(b)所示.

(3)精準確定X射線旋轉180次的方向.

(4)改變了模板上物體的密度.對于與原模板構成物質不同的物體也可以精確測定吸收率，所以改變了模板上物體的密度,得到的模板如圖11(c)所示.

(5)加噪及去噪.對穩(wěn)定性的另一層理解為當外界環(huán)境變化時或因系統(tǒng)某些因素在某范圍或某些數(shù)據(jù)產生異常值或造成數(shù)據(jù)缺失時對系統(tǒng)參數(shù)標定產生的影響大?。畬Ψ€(wěn)定性的另一層理解為當外界環(huán)境變化時或因系統(tǒng)某些因素在某范圍或某些數(shù)據(jù)產生異常值或造成數(shù)據(jù)缺失時對系統(tǒng)參數(shù)標定產生的影響大?。?/p>

所以在上一模板的基礎上對新建模板添加了白噪聲，通過測定CT系統(tǒng)對加噪聲之后得到的接受信息對噪聲的反映是否靈敏來判斷系統(tǒng)的穩(wěn)定度．添加白噪聲后得到的新模板如圖11(d、e)所示：

4 模型的推廣

CT技術可以在不破壞樣品的情況下，利用樣品對射線能量的吸收特性對生物組織和工程材料的樣品進行斷層成像，由此獲取樣品內部的結構信息，是近20年來國際上最先進的非接觸式檢測技術之一.本文就由CT系統(tǒng)的接收信息，通過二維濾波反投影和Radon算法重建了二維斷面圖像，此算法可以推廣到鐵路、地質勘探、石油勘探、航空航天等各個領域．

目前使用的主流標定方法是利用結構己知的標定模體，通過探測器上標定模體的投影對CT成像系統(tǒng)的幾何參數(shù)進行計算，再調節(jié)CT成像系統(tǒng)的關系，提高斷層圖像的重建質量．本文就如何設計一個簡單的標定模型，并盡可能提高標定參數(shù)精度和穩(wěn)定性給出了一個方案，可以在工業(yè)和醫(yī)學檢測中調節(jié)CT成像系統(tǒng)推廣與應用．