扇束CT幾何偽影的校正方法

1、第 5 期2011 年 5 月電子學(xué)報(bào)vol . 39 no . 5may 2011acta el ectronica sinica扇束 ct 幾何偽影的校正方法王亮 ,張朋( 首都師范大學(xué)數(shù)學(xué)科學(xué)學(xué)院 ,北京 100048)摘要 : 實(shí)際工業(yè) ct 系統(tǒng)很難滿足經(jīng)典 ct 圖像重建公式所要求的射線源焦點(diǎn) 、轉(zhuǎn)臺(tái)中心及線陣探測器之間嚴(yán)格的對(duì)準(zhǔn)關(guān)系. 不滿足對(duì)準(zhǔn)關(guān)系將導(dǎo)致經(jīng)典 ct 圖像重建公式所重建的圖像中出現(xiàn)幾何偽影. 本文引入可以描述一般情況下 ct 掃描幾何關(guān)系的一組幾何參數(shù) ,給出了由金屬絲的投影數(shù)據(jù)間接估計(jì)這些幾何參數(shù)的方法 ,還給出了基 于這些幾何參數(shù)的圖像重建公式. 本文只假設(shè)轉(zhuǎn)

2、臺(tái)的旋轉(zhuǎn)軸垂直于由射線源焦點(diǎn)與線陣探測器所形成的扇束 ,而不要 求掃描系統(tǒng)滿足其它幾何對(duì)準(zhǔn)關(guān)系. 對(duì)于掃描系統(tǒng)不對(duì)準(zhǔn)的情況 ,先按照本文方法估計(jì)幾何參數(shù) ,再將其代入本文提 出的圖像重建公式 ,便可以正確重建出被測物的斷層 ct 圖像.關(guān)鍵詞 :扇束 ct ; 圖像重建算法 ; 掃描幾何參數(shù) ; 幾何偽影校正中圖分類號(hào) : tp391文獻(xiàn)標(biāo)識(shí)碼 :文章編號(hào) :037222112 (2011) 0521143207aco rre ctio n of geo metric artif a ct s fo r fa n2b e a m ctwang liang ,zhang peng( school

3、 of mathematical sciences , capital normal university , beijing 100048 , china)ab stract : for a real industrial ct system , it is very difficult to meet the rigorous aligned relations among the x2ray f ocus , the rotation center of turntable and the linear detector , which is required by the conven

4、tional image reconstruction f ormula . misalign2 ment will result in geometric artif acts in the image reconstructed by the conventional image reconstruction f ormula . in this paper , we first introduce a set of geometric parameters which f ully describe general ct scanning configuration , and then

5、 p rop ose an indirect method to estimate the set of parameters f rom the p rojection data of a metal silk . we also p rop ose a generalized filtered back2p rojec2 tion ( fb p) reconstruction f ormula based on such set of parameters . in this paper , we assume that the rotation axis of the ct turnta

6、ble is vertical to the f an2beam f ormed by the x2ray f ocus and the linear detector array , but do not require that the ct scanning system satisfies other alignment relations . for a misaligned scanning system , we can first estimate the parameters by our method , and then reconstruct the image cor

7、rectly using our reconstruction f ormula with the estimated parameters .key word s : f an2beam ct ; image reconstruction algorithm ; scanning geometric parameters ; geometric artif act correction參數(shù) :射線源焦點(diǎn)到轉(zhuǎn)臺(tái)中心的距離、射線源焦點(diǎn)到探測器的距離. 此時(shí) ,用理想的 fbp 重建公式所重建的 ct圖像將包含幾何偽影 ,導(dǎo)致圖像分辨率降低 , 影響對(duì)圖 像的正確判讀. 當(dāng)掃描視野與射線源焦點(diǎn)到探

8、測器的距 離之比較大時(shí) ,圖像的幾何偽影將加重 , 特別是對(duì)圖像中遠(yuǎn)離中心的部分. 因此 ,如何減輕或消除圖像的幾何偽影是 ct 設(shè)備研發(fā)的關(guān)鍵技術(shù)之一.有關(guān)減輕或消除 ct 圖像幾何偽影的研究結(jié)果公 開報(bào)道的相對(duì)較少. 1986 年 gullberg 等2 提出一種帶有轉(zhuǎn)臺(tái)中心偏移量、射線源焦點(diǎn)到轉(zhuǎn)臺(tái)中心和到探測器的距離等多個(gè)參數(shù)的扇束重建公式 ,并給出了參數(shù)的估計(jì) 方法. 該估計(jì)方法需要求解多個(gè)變量的非線性最優(yōu)化問 題 ,魯棒性較差. 陳明3 將轉(zhuǎn)臺(tái)中心的偏移歸結(jié)為探測1 引言基于線陣探測器的斷層工業(yè) ct 目前仍然是工業(yè) 無損檢測的主流設(shè)備. 三代掃描模式是斷層 ct 最為常 用的掃描模

9、式. 三代 ct 掃描模式常用的圖像重建公式 為經(jīng)典的扇束濾波反投影( fbp) 重建公式. 該公式要求 ct 系統(tǒng)滿足嚴(yán)格的幾何成像對(duì)準(zhǔn)關(guān)系1 : (1) 轉(zhuǎn)臺(tái)的旋 轉(zhuǎn)軸垂直于由射線源焦點(diǎn)與線陣探測器所形成的扇束; (2) 射線源焦點(diǎn)與探測器中心的連線垂直于探測器 , 且 轉(zhuǎn)臺(tái)中心位于該垂線上. 此外 , 還要求精確測量 : ( 1) 射 線源焦點(diǎn)到轉(zhuǎn)臺(tái)中心的距離 ; ( 2) 射線源焦點(diǎn)到探測器 的距離. 然而實(shí)際 ct 系統(tǒng)機(jī)械安裝時(shí)很難同時(shí)滿足上 述對(duì)準(zhǔn)關(guān)系 ,也難以直接精確測量圖像重建所需的兩個(gè)收稿日期 :2010211212 ;修回日期 :2011202218基金項(xiàng)目 :國家自然科

10、學(xué)基金(no . 60971131) ;北京市高層次人才項(xiàng)目(no . rhr20110509) ;北京自然科學(xué)基金(no . 3112006)器的偏轉(zhuǎn) ,給出了帶有探測器偏轉(zhuǎn)參量的重建公式 , 但沒有考慮探測器中心的偏移量. 孫靈霞等4 設(shè)計(jì)了一 種射束中心位置測試裝置 , 未涉及其它參數(shù)的測量. 對(duì)于錐束 ct 圓軌跡掃描系統(tǒng) ,noo 等人5 和陳煉等人6分別采用求質(zhì)心的方法 , 求解掃描系統(tǒng)的部分幾何參 數(shù). 兩篇文獻(xiàn)都做了探測器平面平行于轉(zhuǎn)軸的假定. 孫軸 , 其方向與 os 相同; 過 o 且垂直于軸的直線為軸 , 使 o為右手直角坐標(biāo)系. 轉(zhuǎn)臺(tái)坐標(biāo)系 : 以 o 為原 點(diǎn) , o

11、軸順時(shí)針旋轉(zhuǎn)一個(gè)角度得到 ox 軸 , 建立右手直角坐標(biāo)系 oxy . 探測器坐標(biāo)系 : 以 o 為原點(diǎn) , 以探測器d所在直線為坐標(biāo)軸建立坐標(biāo)系 od u .記 r = | os | , o到 ox 的順時(shí)針夾角為, o到怡等人7 ,8o u 的順時(shí)針夾角為. 設(shè) o為 o與 o u 的交點(diǎn) , 記在“4 小球模體所在平面垂直于轉(zhuǎn)臺(tái)平面”的dd假設(shè)下 ,利用一個(gè)角度下的投影數(shù)據(jù)給出了錐束 ct 掃描系統(tǒng)參數(shù)的一種測量方法. 針對(duì)錐束 ct 螺旋掃描系 統(tǒng) ,牛小明等9 給出了帶有探測器中心偏移量的重建方法. 該方法不需要對(duì)投影數(shù)據(jù)進(jìn)行重排就能夠減小 投影數(shù)據(jù)的偏移對(duì)重建圖像的影響. 針對(duì)微焦

12、點(diǎn)錐束ct ,patel 等10 提出一種不用測量模體的自校準(zhǔn)方法 ,該方法對(duì)部分參數(shù)進(jìn)行了測量與校正. 楊民等11 提出了 一種精確測量射線源焦點(diǎn)投影坐標(biāo)的方法.針對(duì)三代 ct 圖像的幾何偽影問題 , 本文引入了可 以完全描述一般情況下的 ct 掃描幾何關(guān)系的一組幾何參數(shù) ,給出了由金屬絲的投影數(shù)據(jù)間接估計(jì)這組幾何參數(shù)的方法 , 還給出了基于這組幾何參數(shù)的圖像重 建公式. 我們給出的幾何參數(shù)間接估計(jì)方法的特點(diǎn)是:(1) 全部掃描幾何參數(shù)均可以從金屬絲模體的投影數(shù) 據(jù)中間接估計(jì)出來 , 克服了幾何參數(shù)無法直接精確測量的困難 ; ( 2) 所有幾何參數(shù)通過解析公式求出 , 不需 要求解非線性最優(yōu)

13、化問題 , 因而求解結(jié)果具有很好的魯棒性. 我們給出的基于這些幾何參數(shù)的圖像重建公式 ,只要求轉(zhuǎn)臺(tái)的旋轉(zhuǎn)軸垂直于由射線源焦點(diǎn)與線陣 探測器所形成的扇束( 否則無法得到同一斷層的掃描數(shù)據(jù)) ,而不要求掃描系統(tǒng)滿足其它幾何對(duì)準(zhǔn)關(guān)系. 仿 真和實(shí)測數(shù)據(jù)的驗(yàn)證結(jié)果表明 : 在射線源焦點(diǎn)、轉(zhuǎn)臺(tái)中心和線陣探測器之間的幾何位置不滿足傳統(tǒng)對(duì)準(zhǔn)關(guān)系的情況下 ,利用我們的方法可以從金屬絲的投影數(shù)據(jù) 中精確估計(jì)出幾何參數(shù) , 將這些幾何參數(shù)代入我們給 出的圖像重建公式就可以正確重建 ct 圖像 ,有效消除 幾何偽影. 當(dāng)掃描視野與射線源焦點(diǎn)到探測器的距離 之比較大時(shí) ,重建的 ct 圖像同樣具有較高的精度.2一般情

14、況下 ct 掃描幾何參數(shù)的間接估計(jì) 方法本節(jié)首先給出一般情況下 , 如何選取幾何參數(shù)和 坐標(biāo)系來刻畫 ct 掃描系統(tǒng); 其次給出所選取的全部幾 何參數(shù)的間接估計(jì)方法.211 一般情況下的掃描幾何參數(shù)和坐標(biāo)系如圖 1 所示 , 以 s , o 和 od 分別表示射線源焦點(diǎn)、 轉(zhuǎn)臺(tái)中心和探測器中心的位置. 我們將按如下方式建 立扇束坐標(biāo)系 ( 固定坐標(biāo)系) 、轉(zhuǎn)臺(tái)坐標(biāo)系和探測器坐標(biāo)系. 扇束坐標(biāo)系 :選取 o 為原點(diǎn); 過 o 與 s 的直線為h = | ood| , d = | so| ( 注意 : 當(dāng)射線源焦點(diǎn)與轉(zhuǎn)臺(tái)中心的連線不垂直于線陣探測器時(shí), d 表示轉(zhuǎn)臺(tái)中心在 探測器上投影點(diǎn) o到射線

15、源焦點(diǎn) s 的距離 , 而不是射 線源焦點(diǎn)到探測器的距離) . 固定參數(shù) r , d , h ,與角 度參數(shù)一起 , 唯一刻畫了 ct 掃描狀態(tài).212 任意點(diǎn) p0 的投影地址設(shè) p0 為物體上遠(yuǎn)離轉(zhuǎn)臺(tái)中心的任意一點(diǎn), 其在oxy 下的坐標(biāo)記為( x0 , y0) ; 當(dāng) ox 與 o夾角為時(shí) , p0在 o下的坐標(biāo)記為 (0 ,0 ,) ; p0 在 od u 上的投影 地址記為 u0 ,. 下面用( x0 , y0) 及參數(shù) r , d , h ,和來 表示 u0 ,和(0 ,0 ,) .由 s 發(fā)出且過點(diǎn) p0 的射線交探測器于點(diǎn) u0 , 在o坐標(biāo)系下可表示為ou0 , = r (

16、0 ,1) + t0 ,( (0 ,0 ,) - r (0 , 1) )( 1)另一方面 , 點(diǎn) u0 ,也可表示為ou0 , = - ( d - r) ( 0 ,1) + ( h + u0 ,) ( n1 , n2) ( 2) 其中( n1 , n2) = ( cos, sin) . 由式( 1) 與( 2) 求解可得 p0 在 od u 上的投影地址 d0 ,u0 , = - h + n( 3)0 ,) 20 ,1 ( r -+ n0 , = x0cos+ y0 sin其中( 4)0 , = y0cos - x0 sin2138 個(gè)對(duì)稱點(diǎn)的投影地址滿足的方程組為了簡單起見 , 我們做變量替

17、換. 令n1n2x0y0n1 = d , n2 = d , x0 = r , y0 =r( 5)0 , = x 0cos+ y0 sin0 , = y 0cos- x0 sin第 5 期王 亮 :扇束 ct 幾何偽影的校正方法1145記( x0 , y0) 所在的坐標(biāo)系為 oxy , (0 ,0 ,) 所在的坐標(biāo)系為 o. ( x 0 , y0 ) 與 ( x0 , y0) , (0 ,0 ,) 與 (0 ,0 ,) 的關(guān)系是固定比例的縮放關(guān)系 , 圖像經(jīng)過該縮放后結(jié)構(gòu)不變. 原來的幾何參數(shù) r , d , h , n1 , n2 中 n1 與 n2 是不獨(dú)立的 , 而新的幾何參數(shù) r , h

18、 , n1 , n2 是彼此獨(dú) 立的. 由 r , h , n1 , n2 與角度參數(shù)一起 , 可以唯一刻畫 ct 掃描狀態(tài). 做變量替換后 , 由式 ( 3) 可以得到 (0 , 0 ,) 在探測器坐標(biāo)下的投影地址0 ,u3 + h2 - m1 + m2u7 + h = -(13)2 + m - m1 2這 4 個(gè)方程中包含 3 個(gè)未知數(shù) h , m , m , 從中任選 3 個(gè)1 2方程聯(lián)立 , 便可消去 m1 和 m2 求得 h . 于是 h 可有以下4 個(gè)表達(dá)式h =- ( u0 + u4) u1 u3 + 2 u0 ( u3 u5 -u1 u7)- 12 u u + 2 u u -

19、u0 (3)(14)1 31h =- ( u2 + u6) u1 u3 - 2 u2 ( u5 u7 -u1 u3)- 12 u1 u3 - 2 u2 ( u1 + u3)u0 , = - h +( 6)(15)n1 ( 1 - 0 ,) + n20 ,h =- ( u1 + u5) u0 u2 + 2 u1 ( u0 u2 -u4 u6)當(dāng) ox 與 o 夾角為時(shí) , ( x 0 , y0 ) 在 o下的坐標(biāo)為(0 ,0 ,) . 令 ( x0 , y 0 ) 順時(shí)針每次旋轉(zhuǎn) / 4 , 旋轉(zhuǎn) 7- 12 u0 u2 - 2 u1 ( u0 + u2)(16)h =次后 得 到 8 個(gè) 對(duì)

20、稱 點(diǎn). 該 8 個(gè) 點(diǎn) 相 應(yīng) 的 (0 ,+ k ,- ( u3 + u7) u0 u2 - 2 u3 ( u2 u4 -u0 u6)4- 12 u0 u2 - 2 u3 ( u0 - u2)0 ,+ k) 之間有如下關(guān)系(17)4其中 u0 = u0 - u4 , u1 = u1 - u5 , u2 = u2 - u6 , u3 = u3 -u7 . 由于投影地址之間存在著內(nèi)在的相容關(guān)系, 因此在cos k- sin k44(0 ,+ k ,0 ,+ k) = (0 ,0 ,),sin kcos k44( x , y 的投影地址 u ( k = 0 ,1 ,2 ,7) 不含誤差時(shí) , 由

21、0)0k44, 7式( 1417) 計(jì)算出的 h 是相同的. 在投影地址 u 存在kk = 0 ,1 ,2 ,( 7)誤差時(shí) , 將式( 1417) 計(jì)算出的 h 進(jìn)行平均 , 可以降低投影地址 uk 誤差引起的 h 的計(jì)算誤差.求解出參數(shù) h 之后 , 現(xiàn)在考慮求解參數(shù) n1 和 n2 . 為此 , 聯(lián)立 u0 與 u4 , u2 與 u6 分別解出兩組0 ,和0 , 的表達(dá)式利用式(7) , 上述 8 個(gè)對(duì)稱點(diǎn)的投影地址可表示為0 ,+ k4 ,u = - h + 0 ,+ kn ( 1 - ) + n 0 ,+ k 2 0 ,+ k144k = 0 ,1 ,2 ,4, 7( 8)214由

22、投影地址估計(jì)幾何參數(shù)的公式本小節(jié)將導(dǎo)出由 8 個(gè)對(duì)稱點(diǎn)的投影地址估計(jì)幾何 參數(shù) h , n1 , n2 和 r 的公式.- 10 , = - 2 h0 h4 n1 u0(18)0 , =0 , =- 1- 2 h0 h4 n2 + h0 + h4u0- 1為簡單起見 , 令 uk = u , k = 0 ,1 ,2 , 7 . 由 u00 ,+ k2 h2 h6 n2 - h2 - h6u24(19)與 u4 , u1 與 u5 , u2 與 u6 , u3 與 u7 可得如下方程組- 10 , = - 2 h2 h6 n1 u2其中 h0 = h + u0 , h1 = h + u1 , h

23、2 = h + u2 , h3 = h + u3 ,h4 = h + u4 , h5 = h + u5 , h6 = h + u6 , h7 = h + u7 . 從中消去0 ,和0 , 即可求解出 n1 和 n2u + h0 ,+ k 4 + hu 0 ,+ ( k + 4) 4n1 ( 1 - 0 ,+) + n20 ,+n1 = u0 u2 ( h h ( h + h ) - h h ( h + h )0 ,+ k( k + 4) ( k + 4) 0 4 2 6 26 0 4444(2 h02 h42 u22 + 2 h22 h62 u02) - 1=(20),n1 ( 1 - k)

24、+ n2 kn2 = ( u02 h h h + h+ u 2 h h h + h0 ,+ 40 ,+ 46 (6)4 (4) )0 ,+ ( k + 4) 2 22 0 04(2 h02 h42 u22 + 2 h22 h62 u02) - 1(21)( 9)k = 0 ,1 ,2 , 3由另外 4 個(gè)角度下的投影 u 與 u , u 與 u 同樣可以1 5 37n2n2解出兩組0 ,和0 ,的表達(dá)式 , 進(jìn)而可解得 n1 和 n2令 m1 = n 0 , +0 ,m2 = n 0 , - 0 ,11n1 = u1 u3h1 h5 ( h3 + h7) - h3 h7 ( h1 + h5)則

25、方程組( 9) 可化簡為2 h12 h52 u32 + 2 h32 h72 u12- 1(22) u0 + h1 - m2u4 + h = - 1 + m2( 10)u12 h3 h7 ( h3 + h7) + u32 h1 h5 ( h1 + h5)n =22 2 22 2 2 - 1 u2 + h1 - m12 h1 h5 u3 + 2 h3 h7 u1(23)u6 + h = - 1 + m1( 11)同樣 , 由于投影地址間存在著內(nèi)在相容關(guān)系 , 在投影地址 uk ( k = 0 ,1 , 2 ,7) 不含誤差時(shí) , 由式( 20) ( 22) , u1 + h 2 - m1 - m2

26、( 12)u5 + h = -( 21) ( 23) 分別計(jì)算出的 n , n 是相同的. 在投影地址 u2 + m + m1 21 2k存在誤差時(shí) , 將式 ( 20) ( 22) , ( 21) ( 23) 計(jì)算出的 n1 , n2進(jìn)行平均 , 可以降低投影地址 uk 誤差引起的 n1 , n2 的 計(jì)算誤差.在實(shí)際應(yīng)用中 , 我們選取一根細(xì)的鉛或鎢的金屬 絲( 或水銀溫度計(jì)) 盡量垂直地立于轉(zhuǎn)臺(tái)上遠(yuǎn)離中心的位置 , 并保持金屬絲 ( 或水銀溫度計(jì)) 的投影位于探測器視野內(nèi). 我們采用三代 ct“轉(zhuǎn)?！辈蓸幽Ｊ?, 即轉(zhuǎn)臺(tái)旋 轉(zhuǎn)一個(gè)角度后停止 , 探測器開始采集穿過被測物的射 線 , 當(dāng)采

27、集結(jié)束后再旋轉(zhuǎn)到下一個(gè)角度位置采集數(shù)據(jù),直至完成一周采樣. 旋轉(zhuǎn)一周的角度采樣數(shù)為8 的倍數(shù) ( 如 720 或 1800) . 從金屬絲( 或水銀溫度計(jì)) 掃描數(shù)據(jù)圖 像( sinogram) 中提取出其中心的投影地址. 利用式( 1423) , 從每個(gè) 0 ,/ 4) 相應(yīng)的 8 個(gè)投影地址 , 分別計(jì) 算出參數(shù) h , n1 和 n2 ; 為提高數(shù)值計(jì)算的穩(wěn)定性 , 將相 應(yīng)不同求得的參數(shù) h , n1 和 n2 進(jìn)行平均.接下來考慮計(jì)算參數(shù) r . 注意到由 h , n1 , n2 和式(18 ,19) 可以求出0 ,和0 , 因此 , 利用兩根標(biāo)定了距 離的金屬絲就可以求出參數(shù) r

28、, 具體方法如下 :設(shè)兩根平行的金屬絲之間的距離為 d , ( x1 , y1 ) 和( x2 , y2) 分別為兩根金屬絲中心在 ct 斷層中的坐標(biāo) , 它 們?cè)?o下的坐標(biāo)分別記為(1 ,1 ,) 和 (2 ,2 ,) ,利用 213 小節(jié)中變換式( 5) 可得 , 在 o下的坐標(biāo)分別記為(1 ,1 ,) 和(2 ,2 ,) . 易見有下列關(guān)系下面通過建立平行束下參數(shù)( r , ) 與扇束下參數(shù)( u ,) 之間的對(duì)應(yīng)關(guān)系 , 并利用坐標(biāo)變換 , 推導(dǎo)基于參 數(shù) h , n1 , n2 的扇束 ct 濾波反投影圖像重建公式.第一步 , 推導(dǎo) ( r , 0的方法可得. 用求極值r ( h

29、+ u) n1(28)r =( h + u) 2 n12 +2d - ( h + u) n2設(shè) q 為原點(diǎn) o 到射線 su 的垂足 , 可求得 oq 在 o下的表達(dá)式為d - ( h + u) n2r ( h + u) n1(,) =,2( h + u) 2 n12 +d - ( h + u) n22d - ( h + u) n2rr -(29)2( h + u) 2 n12 +d - ( h + u) n2另一方面cossin- sin cosr ( cos , sin ) = (,)由式( 2830) 解得(30)( h + u) n1cos+ ( d - ( h + u) n2) si

30、n = arctan(31)- ( h + u) n1 sin+ ( d - ( h + u) n2) cos式( 28) ( 31) 即是參數(shù)( r , ) 與參數(shù)( u ,) 之間的關(guān)系.第二步 , 確定 d rd 與 d ud的關(guān)系. 由式( 28) 和式( 31) 可求得雅克比行列式( x1 , y1) - ( x2 , y2)d =(1 ,1 ,) - (2 ,2 ,)(1 ,1 ,) - (2 ,2 ,)= r( 24)d - ( h + u) n25 r , rdn1 ( )利用本節(jié)方法可以求出幾何參數(shù) h , n1 , n2 , 而由 h , n1 ,n2 和式( 18 、19

31、) 可以求出(1 ,1 ,) 和(2 ,2 ,) , 因此可得到射線源焦點(diǎn)到轉(zhuǎn)臺(tái)中心的距離為d= (32)5 ( u ,)322 n 2 +d - ( h + u) n2( h + u)21第三步 , 將 x0cos + y0 sin - r 用 u , r , d , h , n1 ,n2 表示. 點(diǎn) p0 在 oxy 下的坐標(biāo)為 ( x0 , y0 ) , 由式 ( 4) 可得(0 ,0 ,) 所表示的( x0 , y0) 為x0 = 0 ,cos- 0 ,sinr =( 25)(1 ,1 ,) - (2 ,2 ,)(33)3 基于參數(shù) h , n1 , n2 的圖像重建公式下面我們從平行

32、束的濾波反投影 ct 圖像重建公 式出發(fā) , 推導(dǎo)一般情況下扇束掃描基于參數(shù) h , n1 , n2 的濾波反投影圖像重建公式.平行束的濾波反投影 ct 圖像重建公式可以寫成 如下形式y(tǒng)0 = 0 ,sin+ 0 ,cos由式( 2830) 可得( cos , sin ) = - ( h + u) n1 sin+ ( d - ( h + u) n2) cos,( h + u) 2 n12 +2d - ( h + u) n2( h + u) n1cos+ ( d - ( h + u) n2) sin(34)2 + ( h + u) 2 n12 + ( d - ( h + u) n2) 2 1f

33、( x0 , y0) = 2 f ( r , ) h( x0cos + y0 sin - r) d rd 0 - 由式( 28 , 33 ,34 得)( 26) 其中 , ( x0 , y0) 為待重建點(diǎn) , r 為坐標(biāo)原點(diǎn) o 到射線 su 的距離 , 為射線法向與 x 軸方向的夾角 ( 如圖 1 所- hn1 ( r - 0 ,) + ( d -hn2) 0 ,x0cos + y0 sin - r =n1 ( r - 0 ,) + n20 , u- 12 22 ( h + u)n1 + ( d - ( h + u) n2)(35)2+ 示) , f ( r , ) 是以( r , ) 為變

34、量的投影數(shù)據(jù) , h( r) = 由投影地址公式( 3) 得- hn1 ( r - 0 ,) + ( d - | | ei2rd.hn2) 0 ,op第 5 期王亮 :扇束 ct 幾何偽影的校正方法1147n1 ( r - 0 ,) + n20 , u0 ,所代表的長度進(jìn)行標(biāo)定時(shí)才需要 r 的值. 如果模體的長度已知 , 通過比較該模體 ct 圖像的像素單位長度和 它的實(shí)際長度 , 也可推出 r 的值.4實(shí)驗(yàn)驗(yàn)證本節(jié)我們分別采用仿真數(shù)據(jù)和實(shí)測工業(yè) ct 數(shù)據(jù) ,驗(yàn)證本文提出的幾何偽影校正方法的有效性.411仿真數(shù)據(jù)驗(yàn)證用直徑為 0175mm 的小圓模擬金屬絲的橫截面. 系 統(tǒng)的 掃 描 幾 何

35、 參 數(shù) 分 別 設(shè) 定 為 r = 1000mm , d =1200mm. 線陣探測器由 1400 個(gè)長度為 0125mm 的單元組成 ,總長度 350mm. 此時(shí)放大比為 112 , 像素單元尺寸為012mm ,掃描視野直徑為 2881614mm. 角度采樣數(shù)為 1800 .為驗(yàn)證本文參數(shù)間接估計(jì)方法的正確性和穩(wěn)定 性 , 選取 3 組參數(shù)進(jìn)行數(shù)值實(shí)驗(yàn). no . 1 組參數(shù)值及金屬絲的 位 置 分 別 設(shè) 定 為 : h = 2mm , n1 = 0199996 , n2 =0100873 (= 015) , x0 = 130mm , y0 = 40mm , 模擬金屬絲 的斷層掃描數(shù)據(jù)

36、, 如圖 2 所示. 圖中橫向?yàn)樘綔y器坐標(biāo) ,縱向?yàn)椴蓸咏嵌? no . 2 組和 no . 3 組參數(shù)值及金屬絲的 位置分別設(shè)定為 : h = 4mm , n1 = 0199985 , n2 = 0101745 (= 1) , x0 = - 120mm , y0 = - 60mm 和 h = 6mm , n1 =( 36)將式(36) 代入(35) 右端 , 則n1 ( r - 0 ,) + n20 ,( u0 , - u)x0cos + y0 sin - r =2( h + u) 2 n12 +d - ( h + u) n2( 37)又濾波函數(shù) h( 3 ) 為負(fù)二次齊次 , 因此有h(

37、x0cos + y0 sin - r) =2( h + u) 2 n12 +d - ( h + u) n2h( u0 , - u)( 38)n1 ( r - 0 ,) + n20 ,2將式( 32 , 38) 代入平行束重建公式( 26) , 得到帶有掃描幾何參數(shù) r , d , h , n1 , n2 的扇束濾波反投影重建 公式2 + 1f ( x0 , y0) = 2 p ( u ,)0 - d - ( h + u) n22( h + u) 2 n12 +d -( h + u) n2rdn1 h( u0 , - u)d ud( 39)n1 ( r - 0 ,) + n20 ,2其中 p (

38、 u ,) 是以( u ,) 為變量的扇束投影數(shù)據(jù).將重建公式( 39) 右端第一個(gè)分式的分子分母同時(shí) 除以 d , 第二個(gè)分式的分子分母同時(shí)除以 r2 d2 , 公式兩 端再同時(shí)乘 r , 由變量替換式( 5) 可得0199939 , n = 0103489 (= 2) , x = 95mm , y= - 95mm.200其斷層掃描數(shù)據(jù)圖與圖 2 相似 , 只是起始位置不同.2 1rf ( x0 , y0) = 2 p ( u ,)0 - 1 - ( h + u) n2( h + u) 2 n12 + ( 1 - ( h + u) n2) 2 n1 h( u0 , - u) 2 d ud(

39、40)( n1 ( 1 - 0 ,) + n20 ,)將式( 40) 右端記作新的函數(shù) f ( x0 , y 0 ) , 得到基于參數(shù) h , n1 和 n2 的扇束 ct 濾波反投影圖像重建公式2 f ( x0 , y 0) = 12 0 - p ( u ,)1 - ( h + u) n2采用圖像處理的方法提取質(zhì)心 , 可以得到投影圖中金屬絲中心的投影地址 , 利用本文第 2 節(jié)提出的間接( h + u) 2 n12 + ( 1 -n1 h( u0 , - u)( h + u) n2) 2估計(jì)方法計(jì)算出參數(shù) h , n 和 n . 為保證計(jì)算結(jié)果的穩(wěn)1 2定性 , 對(duì)于每組實(shí)驗(yàn) , 將 36

40、0 度內(nèi)不同角度對(duì)應(yīng)的估計(jì)值進(jìn)行平均后作為參數(shù)的最終計(jì)算結(jié)果, 結(jié)果見表 1 .由表 1 中 no . 1 組參數(shù)值的實(shí)驗(yàn)結(jié)果可知 : (1) 計(jì)算 出的 h 平均值為 211165mm , 與設(shè)定值 2mm 相比 , 誤差為011165mm. 實(shí)驗(yàn)所設(shè)定的線陣探測器單元長度為0125mm , 該誤差相當(dāng)于探測器單元長度的 01466 倍. ( 2)根據(jù)估計(jì)值 n1 = 813328 10 - 4 , n2 = 718246 10 - 6 可以2 d ud( 41)( n1 ( 1 - 0 ,) + n20 ,)于是有如下關(guān)系f ( x0 , y0) = rf ( x0 , y0)( 42)從

41、關(guān)系 x0 = x0/ r , y0 = y0/ r 和式 ( 42) 可知 , r 的值不影響f ( x 0 , y 0) 和 f ( x0 , y0) 量化后的圖像值 , 而只影響它 們相應(yīng)圖像的整體縮放比例 , 即 f ( x0 , y0 ) 和 f ( x0 , y0 )相應(yīng)圖像是相似的 , 因此可將重建 f ( x0 , y0 ) 改為重建f ( x0 , y 0) .與式(39) 相比 , 式 ( 41) 中只含有 h , n1 和 n23 個(gè)參 數(shù) , 而且這 3 個(gè)參數(shù)是彼此獨(dú)立的. 此外 , 從式 ( 41) 可 知 , 在圖像重建中并不需要知道 r 的值. 對(duì)圖像上像素1計(jì)

42、算 出 d = ( n 2 + n 2 ) - = 12001024mm , 與 設(shè) 定 值21 21200mm 相比 , 誤差為 01024mm. ( 3) 根據(jù) n , n 可以計(jì)1 2算出探測器的偏轉(zhuǎn)角 = arctan ( n2/ n1 ) = 010094 弧度(01538 度) , 與設(shè)定值 015 度相比 , 誤差為 01038 度. 類 似地 , 可以計(jì)算出表1 中no . 2 ,no . 3組參數(shù)值所對(duì)應(yīng)的h , 根據(jù) n1 , n2 計(jì)算出 d 以及探測器的偏轉(zhuǎn)角.表 1 參數(shù)間接估計(jì)結(jié)果通過本文第 2 節(jié)方法計(jì)算出的 h , n1 和 n2 , 如表 21所示. 由參數(shù)

43、n 和 n 可以計(jì)算出 d = ( n 2 + n 2) - =21 21 2no . 參數(shù)設(shè)定值估計(jì)值( 平均值)絕對(duì)誤差10191165mm , 探測器偏轉(zhuǎn)角= arctan ( n / n ) = 01024852 1h ( mm)n 1n 228 . 333 10 - 47 . 272 10 - 62 . 11658. 3328 10 - 47. 8246 10 - 60 . 11652 . 16719 10 - 85 . 52508 10 - 7弧度( 11424 度) .為驗(yàn)證參數(shù)估計(jì)結(jié) 果的正確性 , 在 同 樣 的掃描條件下 , 掃 描 兩 個(gè)1表 2 ct 系統(tǒng)參數(shù)估計(jì)結(jié)果h

44、 ( mm)n 1n 248 . 332 10 - 41 . 454 10 - 54. 127238 . 33178 10 - 41. 4356 10 - 50 . 127238 . 21602 10 - 81 . 88642 10 - 72易拉罐 一個(gè)靠近轉(zhuǎn)臺(tái)(h ( mm)n 1n 268 . 328 10 - 42 . 908 10 - 56. 128488 . 32762 10 - 42. 9267 10 - 50 . 12848- 6. 34234 10 - 81 . 84054 10 - 7中心 , 一個(gè)遠(yuǎn)離 轉(zhuǎn) 臺(tái) 中3心) , 獲得其掃描數(shù)據(jù). 圖 3 給出了下列三種方法的重建

45、結(jié)果 :方法 1 為理想的 fbp 重建算法; 方法 2 為調(diào)節(jié) 探測器中心偏移量(offset) 的重建算法; 方法 3 為本文提412實(shí)測數(shù)據(jù)驗(yàn)證實(shí)際 ct 系統(tǒng)所用的探測器單元寬度為 01083mm ,探測器單元個(gè)數(shù)為 3710 ; x 射線源焦點(diǎn)的 1/ 10 高寬為3mm ;射線源焦點(diǎn)到轉(zhuǎn)臺(tái)中心的距離約為 800mm ; 射線 源焦 點(diǎn) 到 探 測 器 的 距 離 約 為 1000mm. 將 直 徑 約 為015mm 的細(xì)鉛絲立于轉(zhuǎn)臺(tái)上遠(yuǎn)離旋轉(zhuǎn)中心的位置處進(jìn)行掃描 ,角度采樣數(shù)為 1800 . 從鉛絲的掃描數(shù)據(jù)中可以 提取出不同角度下的中心位置.出的基于參數(shù) h , n 和 n 的扇

46、束濾波反投影重建算1 2法. 圖 3 ( a) 為方法 1 的重建結(jié)果 , 從其局部放大圖 3( a1) 、3 ( a2) 中可以看出易拉罐壁存在明顯的幾何錯(cuò)位 偽影 , 這是由掃描幾何誤差所引起的. 圖 3 ( b) 為方法 2的重建結(jié)果 , 從其局部放大圖 3 ( b1) 、3 ( b2) 中可以看出靠近轉(zhuǎn)臺(tái)中心處的易拉罐壁基本正確, 但遠(yuǎn)離轉(zhuǎn)臺(tái)中 心處的易拉罐壁仍有明顯的幾何錯(cuò)位偽影. 圖 3( c ) 為方法 3 的重建結(jié) 果 , 從 其 局 部 放 大 圖 3( c1) 、3 ( c2) 中可以看出 重建的易拉罐壁在靠近轉(zhuǎn)臺(tái)中心處和遠(yuǎn)離轉(zhuǎn)臺(tái) 中心處均無明顯偽影 , 說明本文方法可以對(duì)

47、掃描 系統(tǒng)幾何位置偏差所引起的偽影進(jìn)行有效的校正.5總結(jié)本文提出了扇束 ct 幾何偽影的校正方法 ,該 方法只要求轉(zhuǎn)臺(tái)的旋轉(zhuǎn) 軸垂直于由射線源焦點(diǎn) 與線陣探測器所形成的 扇束 ,而不要求掃描系統(tǒng) 滿足其它幾何對(duì)準(zhǔn)關(guān)系. 本文重建公式所需 的幾何參數(shù)可以從金屬 絲的投影數(shù)據(jù)中間接估 計(jì) ,從而克服了經(jīng)典重建 公式中所需參數(shù)無法直 接精確測量的困難. 由于參數(shù)估計(jì)值( 平均值)h (mm)0. 328265n19 . 809 10 - 4n22 . 438 10 - 5第 5 期王 亮 :扇束 ct 幾何偽影的校正方法1149所有參數(shù)可通過解析公式求出 , 不需要求解非線性最優(yōu)化問題 ,因而求解結(jié)

48、果具有很好的魯棒性. 而且該間 接估計(jì)方法還降低了對(duì)模體制作精度和復(fù)雜度的要求.實(shí)驗(yàn)表明在射線源焦點(diǎn)、轉(zhuǎn)臺(tái)中心和線陣探測器不 滿足傳統(tǒng)對(duì)準(zhǔn)關(guān)系的情況下 ,利用本文提出的幾何參數(shù)的間接估計(jì)方法及基于這些幾何參數(shù)的 ct 圖像重建公 式 ,仍可正確重建 ct 圖像. 當(dāng)掃描視野與射線源焦點(diǎn)到探測器的距離之比較大時(shí) ,重建的 ct 圖像也具有較高 精度.本文間接估計(jì)方法假定了角度采樣數(shù)是8 的倍數(shù) ,即投影具有 8 對(duì)稱性. 實(shí)際系統(tǒng)如采用閉環(huán)控制的“轉(zhuǎn) ?！蹦Ｊ?,則可以精確做到角度采樣數(shù)是8 的倍數(shù). 如果采用“連續(xù)”模式 , 轉(zhuǎn)臺(tái)的運(yùn)動(dòng)和探測器采樣可能不耦 合 ,造成一周內(nèi)的角度采樣數(shù)可能不是

49、8 的倍數(shù). 此時(shí) ,可以先采用閉環(huán)控制的“轉(zhuǎn)?！辈蓸幽Ｊ讲杉饘俳z的 掃描數(shù)據(jù) ,進(jìn)而間接估計(jì)出全部幾何參數(shù). 在“連續(xù)”掃描模式下 ,將這些幾何參數(shù)代入重建公式中即可重建高精度的 ct 圖像. 也可以將 ( 8) 改為具有準(zhǔn)對(duì)稱關(guān)系 的 8 個(gè)投影 ,導(dǎo)出一般的幾何參數(shù)間接估計(jì)公式.致謝 :感謝首都師范大學(xué)檢測成像實(shí)驗(yàn)室張慧滔、 陳德峰老師以及朱溢佞、趙云松同學(xué)的有益討論及在 數(shù)值實(shí)驗(yàn)方面給予的幫助 , 感謝趙星老師對(duì)本文修改 的建議. 6 陳煉 , 吳志芳 , 劉錫明 , 等. 錐束 ct 系統(tǒng)幾何參數(shù)校正的解析計(jì)算j . 清華大學(xué)學(xué)報(bào) ( 自然科學(xué)版) , 2010 , 50 ( 3)

50、 :418 - 421.chen lian , wu zhif ang ,liu ximing , et al . analytical geometric parameter calibration algorithm f or cone2beam ct j . j ournalof tsinghua university ( science and technology ) , 2010 , 50(3) :418 - 421. (in chinese) 7 sun y , hou y , zhao f y , et al . a calibration method f or mis2 a

51、ligned scanner geometry in cone2beam comp uted tomograp hy j . nd t & e international , 2006 , 39 (6) :499 - 513. 8 sun y , hou y , hu j s . reduction of artif acts induced by mis2 aligned geometry in cone2beam ct j . i eee transactions on biomedical engineering , 2007 , 54 (8) :1461 - 1471. 9 牛小明 ,