基于正弦圖的計算機(jī)斷層圖像配準(zhǔn)

孫晶晶 楊 民 劉靜華 吳文晉

(北京航空航天大學(xué) 機(jī)械工程及自動化學(xué)院,北京 100191)

基于正弦圖的計算機(jī)斷層圖像配準(zhǔn)

孫晶晶 楊 民 劉靜華 吳文晉

(北京航空航天大學(xué) 機(jī)械工程及自動化學(xué)院,北京 100191)

為了解決工業(yè) X射線無損檢測中圖像配準(zhǔn)的問題,以計算機(jī)斷層(CT,Computerized Tomography)圖像中物體的位置變化與采集的投影數(shù)據(jù)之間的理論關(guān)系為基礎(chǔ),提出了基于正弦圖的 CT圖像配準(zhǔn)算法.該算法結(jié)合實際的投影采集系統(tǒng)對投影信號進(jìn)行預(yù)處理,并利用投影信號的相關(guān)性尋找物體的位置變化,可以解決二維平行束和扇束投影采集方式下物體二維剛性變換的配準(zhǔn)問題.由于提出的算法是在重建之前的投影域內(nèi)進(jìn)行,因此相比傳統(tǒng)的圖像域內(nèi)的配準(zhǔn)算法適用性更高,尤其當(dāng)投影數(shù)據(jù)不足、質(zhì)量不高、噪聲較大、重建圖像有嚴(yán)重的偽影時,該方法更加有效可靠.對某一封裝零件的配準(zhǔn)結(jié)果證明了算法的可行性.

計算機(jī)斷層成像;圖像配準(zhǔn);投影

計算機(jī)斷層(CT,Computerized Tomography)圖像的配準(zhǔn)首先是在 20世紀(jì) 90年代初由醫(yī)學(xué)圖像處理發(fā)展起來的一個重要分支,目的是將兩幅圖像中具有解剖意義的診斷點和手術(shù)感興趣的點進(jìn)行匹配,以達(dá)到對病灶精確手術(shù)的目的[1-2].工業(yè) CT應(yīng)用中,也需要將測試件和標(biāo)準(zhǔn)件進(jìn)行對比,以進(jìn)行自動缺陷判別和多余物判別.但是由于測試件和標(biāo)準(zhǔn)件在沒有標(biāo)準(zhǔn)夾具的情況下,掃描時的放置位置很難保證完全相同,從而造成待配準(zhǔn)圖像與模板圖像(即標(biāo)準(zhǔn)試件的斷層圖像)之間產(chǎn)生位置上的差異,因此準(zhǔn)確的圖像配準(zhǔn)是正確判別的基礎(chǔ).傳統(tǒng)的圖像配準(zhǔn)算法[3-8]大多是基于圖像域的配準(zhǔn),對圖像質(zhì)量的依賴性強(qiáng).然而工業(yè) CT檢測的對象為密度較大的金屬件,重建的圖像存在較為嚴(yán)重的環(huán)形偽影、條狀偽影以及噪聲,這些因素會極大影響基于圖像域的配準(zhǔn)算法中極值的尋找,從而得到錯誤的配準(zhǔn)結(jié)果.

為了解決上述問題,本文首先分析了 CT圖像中物體的位置變化與采集的投影數(shù)據(jù)之間的理論關(guān)系[9],給出了基于正弦圖的 CT圖像配準(zhǔn)算法;然后結(jié)合實際投影采集系統(tǒng),分析了投影采集中常見的中心偏移以及采集條件變化問題對該算法的影響,并提出了解決方案.

本文提出的基于正弦圖的 CT圖像配準(zhǔn)算法可以對二維平行束和扇束投影采集方式下兩次掃描物體的相對位置變化做出準(zhǔn)確估計,解決物體二維剛性變換的配準(zhǔn)問題[10-11].由于本文提出的算法是在重建之前的投影域內(nèi)進(jìn)行的,因此相比傳統(tǒng)的圖像域內(nèi)的配準(zhǔn)適用性更高,對于投影數(shù)據(jù)不完全、散射硬化現(xiàn)象嚴(yán)重、掃描系統(tǒng)存在偏差的情況時也同樣適用.

1 理論分析

平行束下的投影采集方式如圖 1所示,掃描過程為:探測器和射線源不動,掃描臺繞中心旋轉(zhuǎn),探測器采集掃描臺上的物體斷層在不同旋轉(zhuǎn)角(這里稱為投影角)下的投影數(shù)據(jù).探測器在每個投影角下采集的投影信息為一維數(shù)據(jù),將 360°范圍內(nèi)的所有投影角下的一維投影數(shù)據(jù)依次排列起來,便構(gòu)成二維數(shù)據(jù),稱為正弦圖.

正弦圖坐標(biāo)系為(s,θ),θ為當(dāng)前的投影角,s為投影到旋轉(zhuǎn)中心的距離,如圖 2所示.經(jīng)過空間某一點 (xr,yr)的投影 R[f](θr,sr)如式(1)所示[9]:

其中,δ()為 Dirac函數(shù);u(xr,yr)為線性衰減系數(shù) ;(θr,sr)為投影的坐標(biāo) ;θr為投影角度;sr為投影在探測器上的橫向距離.

由于 CT掃描中,物體位置的變化多是在垂直于旋轉(zhuǎn)軸的平面(xy平面)內(nèi),因此本文只討論物體在 xy平面內(nèi)的二維剛體運(yùn)動和投影變化之間的關(guān)系.物體的二維剛性運(yùn)動是通過基本的平移和旋轉(zhuǎn)變換組合而成的,每一個變換都可以表示成為矩陣的形式,剛性變換則可以表示為一系列矩陣的乘積.點(xr,yr)在 xy平面內(nèi)繞旋轉(zhuǎn)軸 z軸轉(zhuǎn)動 α角,并在 xy平面內(nèi)沿方向角 β移動距離

圖1 平行束下的投影采集方式

圖2 正弦圖及其坐標(biāo)系

旋轉(zhuǎn)矩陣 Rz(α)和平移矩陣 T(r,β)分別為r,坐標(biāo)變?yōu)?xf,yf),相應(yīng)的矩陣變換為

因此坐標(biāo)(xr,yr)與(xf,yf)間的關(guān)系為

此時在平行束投影采集方式下,經(jīng)過點(xf,yf)的投影 R[f](θf,sf)如式 (6)所示:

其中,u(xf,yf)表示線性衰減系數(shù);(θf,sf)表示投影的坐標(biāo);θf表示投影角度;sf表示投影在探測器上的橫向距離.

由式(1)、式 (5)、式 (6)推導(dǎo)可得

由式(7)可以總結(jié)出空間中某點(xr,yr)的投影R[f](θr,sr)和它在 xy平面內(nèi)進(jìn)行二維剛體變換后的點(xf,yf)的投影 R[f](θf,sf)存在著式(8)所示的關(guān)系:

式(8)說明兩幅投影的灰度曲線是一樣的,只不過在 θ坐標(biāo)上相差了剛性旋轉(zhuǎn)變換的角度α,并且 s軸坐標(biāo)之差是與平移變換的幅值、方向角有關(guān)的正弦函數(shù).因此通過對正弦圖的投影信號進(jìn)行信號的自相關(guān)以及互相關(guān)分析,找到兩幅正弦圖信號中相關(guān)的部分,就可以得到物體二維剛性變換的參數(shù),這也是本文基于正弦圖的 CT圖像配準(zhǔn)算法的基本思路.

2 基于正弦圖的配準(zhǔn)算法

由于實際計算機(jī)斷層掃描系統(tǒng)中存在著旋轉(zhuǎn)中心偏移,采集電流、電壓不穩(wěn)定的情況,因此在實際采集的投影中應(yīng)用信號相關(guān)運(yùn)算還需要進(jìn)行以下的修正:

1)CT圖像的重建坐標(biāo)系是以旋轉(zhuǎn)中心為坐標(biāo)原點的,因此如果兩次采集時,系統(tǒng)的旋轉(zhuǎn)中心存在偏移,就意味著兩個重建坐標(biāo)系是不同的,因此不能直接進(jìn)行正弦圖的配準(zhǔn),而是應(yīng)該首先對旋轉(zhuǎn)中心偏移進(jìn)行校正,將兩個重建坐標(biāo)系放在同一個坐標(biāo)原點上,才能進(jìn)行后續(xù)的基于正弦圖的配準(zhǔn);

2)由于在兩次采集的過程中,射線源的電流、電壓和探測器的吸收系數(shù)不能保證完全相同,因此采集到的投影雖然形狀不變,但是在幅值上會有變化,為了更好地分析正弦圖中信號的相似性,需要先對正弦圖進(jìn)行歸一化處理,以背景為 0值,在[0,255]之間對正弦圖進(jìn)行歸一化,可以大大提高相關(guān)性分析的準(zhǔn)確性.

設(shè)工件在位置 1時采集的投影為 Pr,旋轉(zhuǎn)中心偏移量為 mr,在位置 2時采集的投影為 Pf,旋轉(zhuǎn)中心偏移量為 mf,兩次采集的步進(jìn)角均為 Δθ,每次采集投影的幅數(shù) M=360/Δθ.本文提出的基于正弦圖的 CT圖像配準(zhǔn)算法具體描述如下.

1)預(yù)處理.首先利用投影 Pr和 Pf求取各自的旋轉(zhuǎn)中心偏移量 mr和,然后按照 Pr(i,j)=Pr(i,j-mr)和 Pf(i,j)=Pf(i,j-mf)進(jìn)行旋轉(zhuǎn)中心偏移的校正,然后將扇束掃描結(jié)構(gòu)下的投影重排為平行束掃描結(jié)構(gòu)下的投影,最后對正弦圖進(jìn)行歸一化處理.

2)正弦圖的相關(guān)運(yùn)算分析.利用第 1節(jié)所述的待檢測物在平行束投影中做剛體運(yùn)動時,投影數(shù)據(jù)與位置變化之間的關(guān)系,可以首先采用自相關(guān)運(yùn)算處理平行束的正弦圖得到平移不變函數(shù)序列,然后利用物體旋轉(zhuǎn)相當(dāng)于正弦圖起始位置變化的原理,利用互相關(guān)信息找到待檢測正弦圖和模板正弦圖的對應(yīng)角度,然后再次利用自相關(guān)計算找到各個角度下的投影平移距離,這樣就得到了待檢測物體的位置變化參數(shù).其實現(xiàn)流程圖如圖 3所示.

圖3 正弦圖的相關(guān)運(yùn)算流程圖

3 實驗結(jié)果

本文對圖 4所示的某一封裝零件進(jìn)行計算機(jī)斷層掃描以判斷該封裝零件中是否存在多余物.該零件是密封的,內(nèi)部結(jié)構(gòu)不可見,因此實驗中同時提供了一個確定不含任何多余物的標(biāo)準(zhǔn)件作為參考模板,將測試件的斷層圖像與標(biāo)準(zhǔn)件的同層斷層圖像相減,判斷測試件中是否存在多余物.

圖4 測試件與標(biāo)準(zhǔn)件的實物圖

采集標(biāo)準(zhǔn)件和測試件的相同斷層處(圖 4中黑色直線所示)的 720幅投影數(shù)據(jù),射線源到探測器的距離為 1100mm,采集到的正弦圖如圖 5、圖 6所示.采用 R L濾波重建,重建圖像如圖 7、圖 8所示.事實上測試件這一層并不存在多余物,圖 7、圖 8的特征理論上應(yīng)該是一致的.但是由于測試件與標(biāo)準(zhǔn)件都是密封的,并不清楚其內(nèi)部結(jié)構(gòu),再加上沒有專用夾具,兩次掃描兩個工件在轉(zhuǎn)臺上的放置方位沒有辦法保持完全一致,因此圖7和圖 8中圖像的特征出現(xiàn)了位置變化.

采用傳統(tǒng)的 Fourier-Mellin配準(zhǔn)算法對測試件的斷層圖像(圖 8)和模版的斷層圖像(圖 7)進(jìn)行配準(zhǔn),配準(zhǔn)后的圖像如圖 9所示.配準(zhǔn)后圖像的殘差高達(dá) 76.6%,從圖 9中也可以看出相同的特征并沒有得到匹配.這是由于封閉件內(nèi)存在著彈簧、電線等密度較大的物質(zhì),散射、硬化現(xiàn)象嚴(yán)重,在斷層圖像中存在嚴(yán)重的偽影,Fourier-Mellin算法估計參數(shù)時,難以準(zhǔn)確地收斂到全局極值,所以配準(zhǔn)的結(jié)果并不準(zhǔn)確.

圖5 標(biāo)準(zhǔn)件在斷層處采集的正弦圖

圖6 測試件在斷層處采集的正弦圖

圖7 標(biāo)準(zhǔn)件的斷層圖像

圖8 測試件的斷層圖像

采用本文提出的基于正弦圖的配準(zhǔn)方法進(jìn)行配準(zhǔn)后的圖像如圖 10所示,配準(zhǔn)后的殘差只有7.87%,從圖 10中也可以看出相同的特征很好地匹配在了一起.采用兩種方法計算所得的變換參數(shù)如表 1所示.

圖9 采用Fourier-Mellin算法配準(zhǔn)后的圖像

圖10 采用本文提出的算法配準(zhǔn)后的圖像

表 1 兩種配準(zhǔn)算法計算的變換參數(shù)

表 1可見,本文提出的基于正弦圖的 CT圖像配準(zhǔn)方法,由于是在投影域內(nèi)進(jìn)行的,因此即使重建圖像中存在嚴(yán)重的偽影,基于圖像域的傳統(tǒng)配準(zhǔn)算法已經(jīng)失效的情況下,也能估計出較為準(zhǔn)確的變換參數(shù),使特征得到較好的匹配.

4 結(jié) 論

通過理論分析與實驗驗證,本文提出的基于正弦圖的 CT圖像配準(zhǔn)方法具有以下特點:

1)當(dāng)重建圖像中存在嚴(yán)重的偽影時,基于圖像的配準(zhǔn)算法已經(jīng)失效,該算法依然取得了較好的配準(zhǔn)效果,配準(zhǔn)后圖像殘差僅為 7.87%;

2)本文提出的算法是利用平移后的自相關(guān)序列的殘差最小來尋找旋轉(zhuǎn)角度的,精度以投影采集的步進(jìn)角為單位,如果在精度要求較高的情況下,可以通過對投影進(jìn)行內(nèi)插減小這一誤差.

References)

[1]Maintz JB A,Viergeveramax A.A survey of medical image registration[J].Medical Image Analysis,1998,2(1):1-36

[2]Mark J,Peter B,Michael B,et al.Improved optimization for the robust and accurate linear registration and motion correction of brain images[J].Neuro Image,2002,17(2):825-841

[3]Guo Xiaoxin,Xu Zhiwen.An application of Fourier-Mellin transform in image registration[C]//Proceedings of the Fifth International Conference on Computer and Information Technology.Washington DC:IEEE Computer Society,2005:619-623

[4]Yang Z,Cohen F S.Image registration and object recognition using affine invariants and convex hulls[J].IEEE Transition on Image Processing,1999(8):934-946

[5]Foroosh H,Zerubia JB,Berthod M.Extension of phase correlation to subpixel registration[J].IEEE Transition on Image Processing,2002(11):188-200

[6]Alhichri H S,Kamel M.Virtual circ les:A new set of feature for fast image registration[J].Pattern Recognition Letters,2003(4):1181-1190

[7]陳靈娜,盛利元,陳俊熹.CT圖像配準(zhǔn)算法的研究與實現(xiàn)[J].醫(yī)療設(shè)備信息,2004,19(8):32-35 Chen Lingna,Sheng Liyuan,Chen Junxi.Research on CT image registration algorithm and its implementation[J].Information on Medical Equipment,2004,19(8):32-35(in Chinese)

[8]Reddy B S,Chatterji B N.An FFT-based technique for translation,rotation,and scale-invariant image registration[J].IEEE Transition on Image Processing,1996,5(8):1266-1271

[9]MaoWeihua,Li Tianfang,Wink Nicole.CT image registration in sinogram space[J].Medical Physics,2007,34(9):3596-3602

[10]Cain S C,Hayat M M,Armstrong E E.Projection-based image registration in the presence of fixed-pattern noise[J].IEEE Transitions on Image Processing,2001(10):1860-1872

[11]Lanzavecchia S,Tosoni L,Bellon P L.Fast sinogram computation and sonogram-based alignment of images[J].Oxford Journals,1996(12):531-537

[12]李保磊,傅健,黃巧珍,等.基于正弦圖的工業(yè) CT系統(tǒng)轉(zhuǎn)臺旋轉(zhuǎn)中心的自動確定方法[J].航空學(xué)報,2009,30(7):1341-1345 Li Baolei,Fu Jian,Huang Qiaozhen,et al.Method for automatic determination of center of rotation in industrial computed tomography systemsbased on sinogram[J].Acta Aeronautica EtAstronautica Sinica,2009,30(7):1341-1345(in Chinese)

(編 輯:文麗芳)

Computerized tomography image registration based on sinogram

Sun Jingjing Yang Min Liu Jinghua Wu Wenjin

(School of Mechanical Engineering and Automation,Beijing University of Aeronautics and Astronautics,Beijing 100191,China)

In order to resolve the image registration problem in X-ray nondestructive testing,a new computerized tomography image registration algorithm based the mathematical relationship between the parts displacement and the change of the projection data was proposed.The algorithm preprocesses the projection data considering the practical computed tomography system,and utilizes the correlation between the projection data to search for the parts displacement.The algorithm can deal with registration problem of the two dimensional rigid transformation in the parallel and fan beam geometry.Because the proposed algorithm is done in the projection field before image reconstruction,it is more adaptive comparing with the conventional registration algorithms in the image field.Especially when the projection quality is low,there ism is sing projection and high level noise and there are severe artifacts in the reconstructed images,the proposed algorithm is more effective and reliable.The registration results of a sealed apparatus verify the performance of the algorithm.

computerized tomography;image registration;projection

TP 391.41

A

1001-5965(2011)02-0223-04

2009-12-30

國家自然科學(xué)基金資助項目(60872080);航天科技創(chuàng)新基金資助項目(CASC 0410);北京市科技新星基金資助項目(2005A 14)

孫晶晶(1982-),女,新疆哈密人,博士生,jinger0925@126.com.