孫守思，邱 鈞，桂志國(guó)，劉 暢

(1.中北大學(xué)電子測(cè)試技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，山西太原030051;2.中北大學(xué)儀器科學(xué)與動(dòng)態(tài)測(cè)試教育部重點(diǎn)實(shí)驗(yàn)室，山西太原030051;3.北京信息科技大學(xué)應(yīng)用數(shù)學(xué)研究所，北京100101)

發(fā)射式斷層成像［1］中，核素在被測(cè)物內(nèi)發(fā)生衰變時(shí)釋放出正電子，正電子極易與周?chē)娮影l(fā)生湮滅，發(fā)射出共線(xiàn)光子對(duì).探測(cè)器通過(guò)捕獲這樣一對(duì)光子來(lái)偵測(cè)出湮滅事件的發(fā)生，并且捕獲的光子對(duì)數(shù)就是湮滅事件累計(jì)發(fā)生的次數(shù).斷層內(nèi)任意一點(diǎn)處發(fā)生湮滅事件是隨機(jī)的，由于湮滅過(guò)程伴隨著光子對(duì)的出射，因此任意一點(diǎn)發(fā)射出的光子對(duì)也是隨機(jī)的.設(shè)斷層內(nèi)某點(diǎn)處在各個(gè)方向下的光子發(fā)射密度均為λ，該點(diǎn)向任意方向發(fā)射光子這一事件均以固定密度隨機(jī)獨(dú)立地出現(xiàn)，這一隨機(jī)事件在單位時(shí)間內(nèi)出現(xiàn)的次數(shù)近似滿(mǎn)足泊松分布.

通過(guò)泊松概率模型來(lái)描述放射性物質(zhì)發(fā)射粒子這一物理現(xiàn)象，能很好地刻畫(huà)觀測(cè)數(shù)據(jù)的統(tǒng)計(jì)特性，使之更為貼近實(shí)際情況.本文針對(duì)傳統(tǒng)模型［2-3］在離散過(guò)程中存在誤差這一缺點(diǎn)，引入重建點(diǎn)模型［4］，討論了新模型下加窗基函數(shù)［5］的應(yīng)用.最后，給出了基于重建點(diǎn)模型的EM［6-8］迭代格式，并進(jìn)行了實(shí)驗(yàn)驗(yàn)證.

1 傳統(tǒng)離散化模型

在傳統(tǒng)離散化模型中，待重建區(qū)域被剖分成若干矩形像素格，每個(gè)像素格內(nèi)各點(diǎn)的發(fā)射密度均相同.對(duì)于像素格j，設(shè)光子在各個(gè)方向下的發(fā)射密度均為 λj.

在像素格j內(nèi)，任意一點(diǎn) (x，y)向各個(gè)方向發(fā)射的光子均以固定的平均密度 λj隨機(jī)獨(dú)立地出現(xiàn)，那么，這一隨機(jī)事件在單位時(shí)間內(nèi)發(fā)生的次數(shù)(該點(diǎn)在任意方向下發(fā)射出的光子數(shù))Bj(x，y)則服從泊松分布

同樣，在像素格j內(nèi)每一點(diǎn)沿探測(cè)線(xiàn)i方向發(fā)射的光子都是一個(gè)隨機(jī)事件，故在被探測(cè)線(xiàn)i穿過(guò)的這些點(diǎn)處，則發(fā)生了一系列的隨機(jī)事件.這些隨機(jī)事件在單位時(shí)間內(nèi)發(fā)生的次數(shù)(像素格j發(fā)射出的光子被探測(cè)器接收到的數(shù)目)Bij同樣滿(mǎn)足泊松分布

對(duì)式(2)整理可得

式中:aij為探測(cè)線(xiàn) i被像素格 j截取的長(zhǎng)度.

由于傳統(tǒng)的離散模型和中心旋轉(zhuǎn)的掃描模式不匹配，導(dǎo)致了不同角度下到像素格中心點(diǎn)距離相等的射線(xiàn)與像素格的截格長(zhǎng)度差別較大(如圖1所示)，這會(huì)對(duì)線(xiàn)積分的離散化近似帶來(lái)誤差.

圖1 到中心點(diǎn)距離相等的射線(xiàn)Fig.1 Rays with equal distance from center

2 重建點(diǎn)模型

如圖2所示，從待重建圖像中抽出N個(gè)具有代表性的點(diǎn)，這些點(diǎn)均勻地排成橫豎相等的陣列.每個(gè)點(diǎn)都是其周?chē)粋€(gè)理想小區(qū)域的中心點(diǎn)，即每個(gè)小區(qū)域可以由該點(diǎn)描述.在新的模型中，沒(méi)有像素格的概念，每個(gè)重建點(diǎn)直接對(duì)應(yīng)圖像重建中要顯示的像素.

圖2 重建點(diǎn)模型Fig.2 Reconstruction points model

在離散重建點(diǎn)模型中，假設(shè)理想?yún)^(qū)域內(nèi)中心重建點(diǎn)發(fā)射密度已知，那么該區(qū)域內(nèi)各點(diǎn)的發(fā)射密度均可由中心重建點(diǎn)插值求得.本文選用基函數(shù)法進(jìn)行插值，進(jìn)而得到連續(xù)的重建圖像.

2.1 基函數(shù)法

圖像的基函數(shù)可以理解為圖像的基，任何圖像均可由基函數(shù)線(xiàn)性組合而成.對(duì)任一圖像f(x，y)，設(shè)它對(duì)應(yīng)的基函數(shù)為 b(x，y)，則必存在實(shí)數(shù)λ1，λ2，…，λj，使得

成立.這里的λ1，λ2，…，λj即是要求的光子發(fā)射密度.

在重建點(diǎn)離散化模型中，所選取的基函數(shù)b(x，y)要盡量滿(mǎn)足局部完整，即在重建點(diǎn)理想?yún)^(qū)域內(nèi)，基函數(shù)是完整或近似完整的，并且在理想?yún)^(qū)域邊緣處保持連續(xù)平滑.基函數(shù)定義如下:

基函數(shù)在空域內(nèi)的圖形如圖3所示.

圖3 基函數(shù)示意圖Fig.3 Diagram of basis function

在圖3中，基函數(shù)在理想?yún)^(qū)域內(nèi)旋轉(zhuǎn)對(duì)稱(chēng)，理想?yún)^(qū)域半徑為R，中心重建點(diǎn)j到探測(cè)線(xiàn)i的垂直距離為dij，陰影部分為該理想?yún)^(qū)域的基函數(shù)被射線(xiàn)所在平面截得的截面.

在使用基函數(shù)來(lái)描述重建點(diǎn)理想?yún)^(qū)域時(shí)，通常需要對(duì)其進(jìn)行截?cái)?最簡(jiǎn)單的方法是加矩形窗，但會(huì)導(dǎo)致窗外的時(shí)域信息全部消失，從而引起頻域頻譜泄漏.在實(shí)際應(yīng)用中，要選用合理的加窗函數(shù)，從而使頻譜的擴(kuò)散減到最少.

在選擇窗函數(shù)時(shí)，應(yīng)盡量滿(mǎn)足主瓣較寬，旁瓣高度較低，同時(shí)還要結(jié)合基函數(shù)的特點(diǎn)來(lái)加窗，不僅能取得較好的效果，還便于計(jì)算.設(shè)加窗函數(shù)為 w(x，y)，窗的長(zhǎng)度為 2R，加窗后的基函數(shù)為

2.2 基函數(shù)實(shí)例

R.M.Lewitt［9］選取徑向?qū)ΨQ(chēng)函數(shù)作為基函數(shù)，定義如下:

然后選取與基函數(shù)形似的Kaiser-Bessel窗函數(shù)

最后對(duì)基函數(shù)加窗，得到參數(shù)可調(diào)的廣義基函數(shù)，定義如下:

在m=2，α=4和a=1的情形下，對(duì)基函數(shù)、Kaiser-Bessel窗函數(shù)和加窗基函數(shù)的切片做頻譜分析，仿真結(jié)果如圖4所示.

圖4表明，加窗基函數(shù)旁瓣高度明顯降低，頻譜擴(kuò)散也得到了抑制.通過(guò)對(duì)截?cái)嗟幕瘮?shù)加窗，不僅使基函數(shù)在截?cái)嗵幾兊闷交?，在頻域上也會(huì)盡量壓低基函數(shù)旁瓣的高度，使能量集中在主瓣，高頻分量減少，有效地抑制了基函數(shù)之間因旁瓣相互疊加而造成的頻譜混疊現(xiàn)象.因此，加窗基函數(shù)方法在實(shí)際應(yīng)用中切實(shí)可行.

圖4 加窗基函數(shù)頻譜圖Fig.4 Windowed basis function in frequency domain

3 重建點(diǎn)模型的EM迭代重構(gòu)

在重建點(diǎn)離散化模型的理想?yún)^(qū)域內(nèi)，任意一點(diǎn)(x，y)向各個(gè)方向發(fā)射光子均以固定的平均密度λjbw(x，y)隨機(jī)獨(dú)立出現(xiàn)，那么該事件在單位時(shí)間內(nèi)發(fā)生的次數(shù)(該理想?yún)^(qū)域內(nèi)某點(diǎn)在任意方向下發(fā)射出的光子數(shù))Bj(x，y)則服從泊松分布，即

同樣，在重建點(diǎn)j理想?yún)^(qū)域內(nèi)，每一點(diǎn)沿探測(cè)線(xiàn)i方向發(fā)射光子都是一個(gè)隨機(jī)事件，故在被探測(cè)線(xiàn) i穿過(guò)的這些點(diǎn)處，則發(fā)生了一系列的隨機(jī)事件.這些隨機(jī)事件在單位時(shí)間內(nèi)發(fā)生的次數(shù)(該理想?yún)^(qū)域發(fā)射出的光子被探測(cè)器 i接收到的數(shù)目)Bij同樣滿(mǎn)足泊松分布，即

對(duì)式(11)整理得

設(shè)隨機(jī)變量bi為探測(cè)單元i上記錄的光子數(shù)，則有

由于Bij～Poisson(cijλj)相互獨(dú)立，因此有

式中:bi是可觀測(cè)到的光子數(shù)，它是由不可觀測(cè)但卻完整的 Bij組成的.

由于bi是獨(dú)立的泊松變量，因此其聯(lián)合概率密度函數(shù)為

對(duì)聯(lián)合概率密度函數(shù)取對(duì)數(shù)可得

對(duì)ln g(b，λ)求關(guān)于參數(shù)λ的二階偏導(dǎo)，可得

由式 (17)可知，lng(b，λ)呈嚴(yán)格凹性.

在經(jīng)典EM算法中，需要考慮完備的投影數(shù)據(jù)，故對(duì)不可觀測(cè)變量Bij求其聯(lián)合概率密度

對(duì)式(18)取對(duì)數(shù)可得

EM算法:

1)E-step.對(duì)式(19)取在觀斥數(shù)據(jù) bi和當(dāng)前估計(jì)參數(shù)λn下的條件期望.

其中:

式中:R與參數(shù)λ無(wú)關(guān)，視為常量.

2)M-step.使條件期望最大化，對(duì)期望步求偏導(dǎo)得

令式(22)為零并進(jìn)行整理，可得到EM迭代格式

4 數(shù)據(jù)實(shí)驗(yàn)

4.1 模擬數(shù)據(jù)的數(shù)值實(shí)驗(yàn)

本文采用一組模擬數(shù)據(jù)，選擇傳統(tǒng)模型下的經(jīng)典EM算法和重建點(diǎn)模型下的EM算法重建圖像.該模擬數(shù)據(jù)的角度采樣數(shù)為576，平行線(xiàn)采樣數(shù)為721.該實(shí)例中重建圖像采用的分辨率為600×600，如圖5～圖7所示.

圖5 鸚鵡原圖Fig.5 Original image of parrot

圖6 鸚鵡的投影數(shù)據(jù)位圖Fig.6 Projection data bitmap of parrot

圖7 經(jīng)典EM迭代不同輪次的重建結(jié)果Fig.7 Results of classical EM algorithm under different numbers of iterations

本文采用下面兩個(gè)圖像距離的測(cè)量值來(lái)評(píng)價(jià)不同算法的重建圖像質(zhì)量.

1)歸一化均方距離判斷d

式中:tu，v和 γu，v分別表示測(cè)試模型和重建后圖像中第 u行、第v列的像素密度;ˉt為測(cè)試模型密度的平均值;圖像的像素格個(gè)數(shù)為N×N個(gè).d=0表示重建后圖像真實(shí)地再現(xiàn)測(cè)試模型圖像，d值愈大表示兩者的偏差愈大.

2)歸一化平均絕對(duì)值距離判斷r

r=0說(shuō)明沒(méi)有誤差，r增大說(shuō)明誤差增大.經(jīng)典EM迭代一定次數(shù)后的誤差如表1所示.

表1 經(jīng)典EM迭代一定次數(shù)后的誤差分析表Tab.1 Errors of classical EM under some number of i terations

這里選取m=1時(shí)對(duì)應(yīng)的基函數(shù)，對(duì)不加窗基函數(shù)和加窗基函數(shù)分別進(jìn)行了兩組數(shù)據(jù)實(shí)驗(yàn).①不加窗的基函數(shù)，選取α=0的情形.②加窗基函數(shù)，選取α=4的情形.重建點(diǎn)模型不加窗基函數(shù)EM的誤差和重建結(jié)果如表2和圖8所示.重建點(diǎn)模型加窗基函數(shù)EM的誤差和重建結(jié)果如表3和圖9所示.

表2 重建點(diǎn)模型不加窗基函數(shù)EM的誤差分析表Tab.2 Errors of EM with no window basis function under new discrete model

圖8 重建點(diǎn)模型不加窗基函數(shù)EM的重建結(jié)果Fig.8 Results of EM with no window basis function under new discrete model

圖9 重建點(diǎn)模型加窗基函數(shù)EM的重建結(jié)果Fig.9 Results of EM with windowed basis function under new discrete model

表3 重建點(diǎn)模型加窗基函數(shù)EM的誤差分析表Tab.3 Errors of EM with windowed basis function under new discrete model

實(shí)驗(yàn)表明，在相同迭代次數(shù)時(shí)，傳統(tǒng)模型下EM重建圖像與原圖誤差較大，而重建點(diǎn)模型下EM重建圖像的誤差相對(duì)更小，重建后圖像的質(zhì)量也要好于前者.隨著迭代次數(shù)的增加，重建圖像與原圖像的距離逐漸減小，圖像越來(lái)越趨近于原圖，重建效果越來(lái)越好.對(duì)比兩組實(shí)驗(yàn)發(fā)現(xiàn)，加窗基函數(shù)EM重建出的圖像誤差要小于不加窗基函數(shù)EM的，這表明窗函數(shù)有減小誤差的作用.

4.2 實(shí)測(cè)數(shù)據(jù)的數(shù)值實(shí)驗(yàn)

本節(jié)采用實(shí)測(cè)陶瓷葉片的數(shù)據(jù)，重建的圖像需要強(qiáng)調(diào)細(xì)小裂紋特征.選取重建點(diǎn)EM算法重建圖像，比較在不同迭代次數(shù)下圖像的裂紋變化情況.該陶瓷葉片的實(shí)測(cè)數(shù)據(jù)量是576×241，該實(shí)例中采用的圖像分辨率為256×256.實(shí)測(cè)投影數(shù)據(jù)位圖如圖10所示，經(jīng)典EM迭代不同輪次的重建結(jié)果如圖11所示.

同樣，這里選取m=1時(shí)對(duì)應(yīng)的基函數(shù)，并且分別對(duì)不加窗的基函數(shù)和加窗基函數(shù)進(jìn)行兩組數(shù)據(jù)實(shí)驗(yàn).①不加窗的基函數(shù)，選取α=0的情形，如圖12所示.②加窗基函數(shù)，選取α=4的情形，如圖13所示.

由圖12和圖13的像素曲線(xiàn)圖可以看出，隨著迭代次數(shù)的增加，圖像的灰階變化越來(lái)越明顯，細(xì)節(jié)越來(lái)越突出，與經(jīng)典EM算法相比，曲線(xiàn)也較為平滑，這表明重建點(diǎn)模型下的EM算法能突出陶瓷葉片的細(xì)小裂紋特征，并且有抑制噪聲的作用.

圖10 實(shí)測(cè)投影數(shù)據(jù)位圖Fig.10 Measured projection data bitmap

模擬數(shù)據(jù)和實(shí)測(cè)數(shù)據(jù)的數(shù)值實(shí)驗(yàn)表明，重建點(diǎn)模型較之傳統(tǒng)模型減小了誤差，EM算法對(duì)該模型的重構(gòu)是有效的，可以在一定程度上提高成像的精度和質(zhì)量.

圖11 經(jīng)典EM迭代不同輪次的重建結(jié)果Fig.11 Results of classical EM algorithm under different numbers of iterations

圖12 重建點(diǎn)模型不加窗基函數(shù)EM的重建結(jié)果Fig.12 Results of EM with no window basis function under new discrete model

圖13 重建點(diǎn)模型加窗基函數(shù)EM的重建結(jié)果Fig.13 Results of EM with windowed basis function under new discrete model

5 結(jié)論

本文用泊松概率模型對(duì)發(fā)射式斷層成像進(jìn)行了描述，針對(duì)傳統(tǒng)模型的缺點(diǎn)，引入離散重建點(diǎn)模型.在離散重建點(diǎn)模型中，投影系數(shù)矩陣中元素對(duì)應(yīng)著新的物理意義，用重建點(diǎn)到射線(xiàn)的距離的權(quán)重來(lái)代替原來(lái)的截格長(zhǎng)度，能更準(zhǔn)確地描述線(xiàn)積分的離散表達(dá).

在離散重建點(diǎn)模型中引入了加窗基函數(shù)，對(duì)該模型進(jìn)行重構(gòu)，推導(dǎo)并建立了基于新模型的EM算法，這一迭代重構(gòu)方法反映了重建模型和掃描模式的本質(zhì)刻畫(huà).實(shí)驗(yàn)結(jié)果表明算法有效可行.在新模型下，如何選取更為精確刻畫(huà)投影分布特性的基函數(shù)，如何優(yōu)化和調(diào)整迭代構(gòu)造，深入分析算法的收斂性［10-11］，以及建立相應(yīng)的 OSEM［12-15］算法構(gòu)造等值得進(jìn)一步研究，有助于形成應(yīng)用于實(shí)際檢測(cè)領(lǐng)域的快速實(shí)用迭代算法.

［1］Shepp L A，Vardi Y.Maximum likelihood restoration for emission tomography［J］.IEEE Trans.Med.Imaging，1982，1:113-122.

［2］莊天戈.CT原理與算法［M］.上海:上海交通大學(xué)出版社，1992.

［3］ Herman G T.Fundamentals of Computerized Tomography:Image Reconstruction from Projections［M］.Springer:Second Edition，2009.

［4］劉暢.基于計(jì)算點(diǎn)的圖像重建離散化模型及其相關(guān)算法研究［D］.北京:北京信息科技大學(xué)，2012.

［5］Lewitt R M.Alternatives to voxels for image representation in iterative reconstruction algorithms［J］.Physics in Medicine and Biology，1992，37(3):705-716.

［6］Lange K，Carson R.EM reconstruction algorithms for emission and transmission tomography［J］.Journal of Computer Assisted Tomography，1984，8:302-316.

［7］ Yan Ming.EM-type algorithms for image reconstruction with background emission and poisson noise［J］.Lecture Notes in Computer Science，2011，69(38):33-42.

［8］ Teng Yueyang，Zhang Tie.Generalized EM-type reconstruction algorithms for emission tomography［J］.IEEE Transactions on Medical Imaging，2012，31(9):1724-1733.

［9］Lewitt R M.Multidimensional digital image representations using generalized Kaiser-Bessel window functions［J］.Journal of the Optical Society of America，1990，7(10):1834-1846.

［10］ Jiang Ming，Wang Ge.Development of iterative algorithms for image reconstruction［J］.Journal of X-ray Science and Technology(Invited Review)，2002，10:77-86.

［11］Jiang Ming，Wang Ge.Convergence studies on iterative algorithms for image reconstruction［J］.IEEE Transactions on Medical Imaging，2003，22(5):569-579.

［12］邱鈞，徐茂林.由投影重建圖像的對(duì)稱(chēng)塊迭代算法［J］.電子與信息學(xué)報(bào)，2007，29(10):2293-2300.Qiu Jun，Xu Maolin.A method of the symmetric block iterative for image reconstruction［J］. Journal electronics＆ Information Technology，2007，29(10):2293-2300.(in Chinese)

［13］Hudson H M，Larkin R S.Accelerated image reconstruction using ordered subsets of projection Data［J］.IEEE Transaction on Medical Imaging，1994，13(4):601-609.

［14］劉暢，邱鈞.一種基于對(duì)稱(chēng)結(jié)構(gòu)優(yōu)化的 OSEM快速重建算法［J］.CT理論與應(yīng)用研究，2009，18(4):1-8.Liu Chang，Qiu Jun.An symmetric ordered subset expectation maximization accelerated algorithm［J］.CT Theory and Applications，2009，18(4):1-8.(in Chinese)

［15］Dikaios N，F(xiàn)ryer T D.Acceleration of motion-compensated PET reconstruction:ordered subsets-gates em algorithms and a priori reference gate information［J］.Physics in Medicine and Biology，2011，56(6):1695-1715.