一種反投影模型下的圖像重建迭代算法

趙芳薇，邱 鈞，劉 暢

（北京信息科技大學(xué) 應(yīng)用數(shù)學(xué)研究所，北京100101）

0 引 言

圖像重建技術(shù)在工業(yè)、醫(yī)學(xué)等很多領(lǐng)域有著廣泛的應(yīng)用［1］，其本質(zhì)是利用獲取的投影數(shù)據(jù)重構(gòu)被檢測目標(biāo)內(nèi)部的高分辨率斷層圖像.常見的圖像重建方法有解析算法和迭代算法［2-3］，其中迭代算法將投影數(shù)據(jù)對應(yīng)的積分方程離散化，把圖像重建問題直接轉(zhuǎn)化為線性方程組求解問題，在處理不完全投影數(shù)據(jù)及噪聲污染較為嚴(yán)重的投影數(shù)據(jù)時(shí)較解析算法更有優(yōu)勢.目前迭代算法有諸多理論上的研究成果，是一種被廣泛應(yīng)用的圖像重建算法［4-8］.

這里求解線性方程組是不適定問題，且由于噪聲等其他因素的影響，傳統(tǒng)求解線性方程組的方法并不適用.基于超平面投影和目標(biāo)優(yōu)化理論建立的迭代算法，在由投影重建圖像中應(yīng)用廣泛.ART（Algebraic Reconstruction Technique）［9］、SIRT（Simultaneous Iteration Reconstruction Technique）［10］、CAV（Component Averaging）［11］等都是常見的圖像重建迭代算法，不同的迭代算法對應(yīng)著不同的迭代校正格式，也對應(yīng)著對誤差不同的分配校正方式.

現(xiàn)有的迭代校正格式通常是基于對投影模型及其幾何物理意義的分析，在迭代過程中需計(jì)算投影系數(shù)矩陣的投影系數(shù)［12］.其中，對于投影系數(shù)的計(jì)算，常用的一類方法是計(jì)算射線穿過剖分格的截格長度［1，13］，另一類是引入重建點(diǎn)模型，計(jì)算重建點(diǎn)對射線的貢獻(xiàn)因子［14-15］.通常的迭代重建算法可視為用迭代方法給出廣義敏感逆矩陣估計(jì)形式［16］，也有研究者通過重新刻畫圖像重建算子，直接尋找廣義敏感逆矩陣的近似表達(dá)形式，避免了迭代中出現(xiàn)的半收斂現(xiàn)象，并且有利于節(jié)省圖像重建時(shí)間［17］.

本文從與投影模型相對應(yīng)的反投影模型出發(fā)，研究全體掃描射線對重建點(diǎn)的貢獻(xiàn)因子，用迭代方法給出廣義逆矩陣估計(jì)形式，建立了廣義敏感逆矩陣在某種意義下的近似，形成了迭代成像算法.本文關(guān)于反投影模型下快速圖像重建算法的研究，進(jìn)一步可以推廣到廣義敏感逆矩陣的較為精確的迭代估計(jì)，以期形成基于掃描系統(tǒng)特征的廣義逆敏感矩陣形式的計(jì)算模版，從而形成有效的直接計(jì)算成像方法.

1 重建點(diǎn)離散化模型

依據(jù)圖像重建離散化方式的不同，線性方程組也會存在差別.常見的圖像重建離散化模型有基于網(wǎng)格剖分的離散化模型［1，13］和重建點(diǎn)離散化模型［14-15］，本文引入重建點(diǎn)離散化模型，如圖1所示.

圖1 重建點(diǎn)離散化模型示意圖 Fig.1 Diagram of discrete model based on reconstruction points

重建點(diǎn)離散化模型淡化剖分網(wǎng)格的概念，重建點(diǎn)對應(yīng)著像素點(diǎn)，重建點(diǎn)周圍的一個(gè)理想小區(qū)域內(nèi)像素值視為定值，用重建點(diǎn)到射線的距離的函數(shù)刻畫像素對掃描射線的貢獻(xiàn)因子.較之基于網(wǎng)格剖分的離散化模型，該模型更好地刻畫了顯示模式、掃描模式以及坐標(biāo)架之間的關(guān)系，減少了離散誤差，具有很好的適用性.

記待重建圖像被離散化為N維有限向量X=（x1，x2，…，xN）T，xj是第j個(gè)重建點(diǎn)處的像素值，投影數(shù)據(jù)向量為b=（b1，b2，…，bM）T.重建點(diǎn)離散化模型下，圖像重建問題可歸結(jié)為線性方程組求解問題

式中：C=（cij）M×N是投影系數(shù)矩陣，投影系數(shù)cij=f（dij，φ）表示第j個(gè)重建點(diǎn)對第i條射線Li衰減的貢獻(xiàn)因子，dij是第j個(gè)重建點(diǎn)到射線Li的垂直距離，φ標(biāo)定射線入射方向，i∈｛1，2，…，N｝，j∈｛1，2，…，M｝，如圖2所示.

圖2 投影系數(shù)矩陣示意圖 Fig.2 Diagram of the projection coefficient matrix

2 反投影模型

反投影算法的本質(zhì)是斷層平面內(nèi)某一點(diǎn)的值可看作全體是投影數(shù)據(jù)的加權(quán)求和.

設(shè)待重建區(qū)域被離散化為N個(gè)理想小區(qū)域，有M個(gè)有效投影數(shù)據(jù).反投影模型刻畫了依據(jù)不同掃描射線對某重建點(diǎn)的貢獻(xiàn)因子大小，進(jìn)行圖像重建的過程.記tji是第i根射線對第j個(gè)重建點(diǎn)的貢獻(xiàn)因子，則第j個(gè)重建點(diǎn)像素值可記為式中：bi是第i根射線的投影數(shù)據(jù)；tji為反投影系數(shù)，i∈｛1，2，…，N｝，j?｛1，2，…，M｝.如圖3所示.

圖3 反投影系數(shù)矩陣示意圖 Fig.3 Diagram of the back-projection coefficient matrix

式（2）的矩陣表示形式為

是由全體的反投影系數(shù)組成的反投影系數(shù)矩陣［18］.

在全體投影數(shù)據(jù)對重建點(diǎn)的貢獻(xiàn)因子已知的情形下，圖像重建問題歸結(jié)為由全體投影向量b和反投影系數(shù)矩陣T直接計(jì)算圖像X的問題.

3 反投影系數(shù)矩陣T的估計(jì)

利用反投影模型對圖像進(jìn)行重建的關(guān)鍵是對反投影系數(shù)矩陣給出精確的估計(jì)，本文依據(jù)反投影模型中掃描射線和重建點(diǎn)間的幾何位置關(guān)系，及其反投影系數(shù)的幾何物理意義給出一種合理的估計(jì)方法.

反投影系數(shù)tji是第i根射線對第j個(gè)重建點(diǎn)的貢獻(xiàn)因子，在計(jì)算全體射線對某固定重建點(diǎn)的貢獻(xiàn)因子時(shí)，由于距離重建區(qū)域中心不同的掃描射線穿過待重建區(qū)域的長度不同，需先對全體投影數(shù)據(jù)進(jìn)行均衡.記第i條射線Li被待重建區(qū)域所截長度是si，探測器檢測到的投影數(shù)據(jù)為bi，則投影數(shù)據(jù)可被均衡化為bi／si.均衡化后的M維投影向量可記b′=（b1／s1，…，bM／sM）T∈RM.

對于均衡化后的投影數(shù)據(jù)，不同射線對重建點(diǎn)的貢獻(xiàn)因子與射線到重建點(diǎn)的距離有關(guān)，且隨著射線到重建點(diǎn)的距離逐漸變小而逐漸變大.記f（dij）是第i根射線到第j個(gè)重建點(diǎn)的距離dij的函數(shù)，在考慮均衡化后的全體射線對第j個(gè)重建點(diǎn)的貢獻(xiàn)系數(shù)時(shí)，第i根射線對第j個(gè)重建點(diǎn)的貢獻(xiàn)因子，即投影系數(shù)tji的表達(dá)形式為

式中：si是第i根射線被重建區(qū)域截得的長度；dij是第j個(gè)重建點(diǎn)到第i根射線的距離；f（dij）是dij的函數(shù).si與dij如圖4所示.

圖4 反投影系數(shù)矩陣影響參數(shù)示意圖 Fig.4 Diagram of the influence parameters in back-projection coefficient matrix

4 反投影模型下的圖像重建迭代算法

反投影模型下的圖像重建迭代算法依賴于對反投影系數(shù)的估計(jì).算法的分量迭代格式為

式中：M為投影射線總數(shù)；N為重建點(diǎn)的總數(shù)；bi為投影數(shù)據(jù)為計(jì)算投影值；λk是松弛

因子，i=1，2，…，M，j=0，1，2，…，N，k=0，1，2，….

迭代算法步驟為：

1）賦予待重建圖像向量初始值X（0）=（0，0，…，0）T；

2）依據(jù)初始值X（0），計(jì)算投影值與投影數(shù)據(jù)bi，按照式（5）進(jìn)行迭代.k=0，1，2，….

5 算法驗(yàn)證的數(shù)據(jù)實(shí)驗(yàn)

本節(jié)采用兩組模擬數(shù)據(jù)驗(yàn)證算法，數(shù)據(jù)實(shí)驗(yàn)中f（dij）選取線性插值函數(shù).

模擬數(shù)據(jù)角度采樣數(shù)為576，平行線采樣數(shù)為721，模型和重建結(jié)果的圖像分辨率均為600×600.

圖5 Shepp-Logan原圖及第200行位置 Fig.5 Diagram for the Shepp-Logan and the position of line 200

實(shí)驗(yàn)引入歸一化均方距離d，歸一化平均絕對值距離r和最壞情況距離e對重建結(jié)果進(jìn)行評價(jià).

5.1 模擬數(shù)據(jù)實(shí)驗(yàn)1

模擬實(shí)驗(yàn)1采用經(jīng)典的Shepp-Logan模型.迭代5，10，20，30次的重建結(jié)果及第200行對應(yīng)的像素曲線如圖7所示，與原圖的距離如表1所示.

圖6 Shepp-Logan的投影數(shù)據(jù)位圖 Fig.6 The bitmap of the Shepp-Logan's projection data

圖7 實(shí)驗(yàn)1反投影模型下迭代算法迭代不同輪次的重建結(jié)果 Fig.7 Reconstruction results of the new iterative algorithm under different numbers of iterations for experiment one

表1 實(shí)驗(yàn)1反投影模型下迭代算法誤差分析 Tab.1 Errors analysis of the new iterative algorithm for experiment one

5.2 模擬數(shù)據(jù)實(shí)驗(yàn)2

模擬數(shù)據(jù)實(shí)驗(yàn)2選取鸚鵡模型如圖8和圖9所示.迭代5，10，20，30次的重建結(jié)果及第200行對應(yīng)的像素曲線如圖10所示，與原圖的距離如表2所示.

圖8 鸚鵡原圖及第200行位置 Fig.8 Diagram for the parrot and the position of line 200

圖9 鸚鵡的投影數(shù)據(jù)位圖 Fig.9 The bitmap of the parrot's projection data

圖10 實(shí)驗(yàn)2反投影模型下迭代算法迭代不同輪次的重建結(jié)果 Fig.10 Reconstruction results of the new iterative algorithm under different numbers of iterations for experiment two

兩組模擬實(shí)驗(yàn)的數(shù)據(jù)結(jié)果表明：隨著迭代輪次的增加，成像精度逐步提高，細(xì)節(jié)也愈加明顯.從像素曲線可以看出重建結(jié)果有較好的密度分辨率，空間分辨率也隨著迭代次數(shù)的增加有所改善.

表2 實(shí)驗(yàn)2反投影模型下迭代算法誤差分析 Tab.2 Errors analysis of the new iterative algorithm for experiment two

6 結(jié) 論

本文從與投影模型相對應(yīng)的反投影模型出發(fā)，給出了一種反投影模型下的圖像重建迭代算法.算法基于反投影系數(shù)矩陣的估計(jì)，是廣義敏感逆矩陣在某種意義下的近似，反投影系數(shù)的估算考慮了全體掃描射線對重建點(diǎn)的貢獻(xiàn)因子，實(shí)現(xiàn)了誤差的全局分配，重建結(jié)果具有較好的密度分辨率.全體掃描射線對重建點(diǎn)的貢獻(xiàn)因子，是投影與反投影模型中的基本刻畫因子，反映了投影過程中重建點(diǎn)的物理與幾何的基本性質(zhì)和特征，建立對其的刻畫模型是用迭代方法給出廣義敏感逆矩陣估計(jì)形式的關(guān)鍵.

反投影系數(shù)的估計(jì)方法直接影響算法的收斂速度與成像精度，如何給出更合理的反投影系數(shù)估計(jì)方法，進(jìn)一步有效地提高重建結(jié)果的空間分辨率，值得進(jìn)一步研究.本文提出的反投影模型下的圖像重建迭代方法，可以推廣到廣義敏感逆矩陣的較為精確的迭代估計(jì)，以期形成基于掃描系統(tǒng)特征的廣義逆敏感矩陣形式的計(jì)算模板，從而形成有效的快速實(shí)用的直接計(jì)算成像方法，進(jìn)一步豐富研究的內(nèi)容.

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

［2］Herman G T.Fundamentals of computerized tomography：image reconstruction from projections［M］.2ed.Germany：Springer，2009.

［3］Jiang Ming，Wang Ge.Development of iterative algorithms for image reconstruction［J］.Journal of X-ray Science and Technology，2002，10（1／2）：77-86.

［4］Herman G T.Algebraic reconstruction techniques can be made computationally efficient［J］.Medical Imaging（0278-0062），1993，12（3）：600-609.

［5］Jiang Ming，Wang Ge.Convergence of iterative algo-rithms for image reconstruction［C］.Bio-medical Imaging，2002：685-688.

［6］邱鈞，徐茂林.由投影重建圖像的對稱塊迭代算法［J］.電子與信息學(xué)報(bào)，2007，29（10）：2296-2300.Qiu Jun，Xu Maolin.A method of symmetric block-iterative for image reconstruction［J］.Journal of Electronics ffInformation Technology（1009-5896），2007，29（10）：2296-2300.（in Chinese）

［7］畢丹丹，李筠.有關(guān)ART算法中松馳因子的研究［J］.信息技術(shù)，2014（9）：121-124，128.Bi Dandan，Li Jun.Research on the relaxation parameter of algebraic reconstruction technique［J］.Information Technology，2014（9）：121-124，128.（in Chinese）

［8］黨曉軍，韋孟伏.CT代數(shù)重建算法的快速實(shí)現(xiàn)［J］.核電子學(xué)與探測技術(shù)，2011，31（10）：1104-1106.Dang Xiaojun，Wei Mengfu.A fast implementation of CT statistical reconstruction methods［J］.Nuclear Electronics ffDetection Technology，2011，31（10）：1104-1106.（in Chinese）

［9］Gorden R，Bender R，Herman G T.Algebraic reconstruction techniques（ART）for three-dimensional electron microscopy and X-ray photograph［J］.Journal of Theoretical Biology，1970，29（3）：471-476.

［10］Gilbert P.Iterative methods for the three-dimensional reconstruction of an object from projections［J］.Journal of Theoretical Biology，1972，36（2）：105-117.

［11］Censor Y，Gordon D，Gordon R.Component averaging：an efficient iterative parallel algorithm for large and sparse unstructured problems［J］.Parallel Computing，2001，27：777-808.

［12］陳建林，閆鑌，李雷，等.CT重建中投影矩陣模型研究綜述［J］.CT理論與應(yīng)用研究，2014，23（2）：317-328.Chen Jianlin，Yan Bin，Li Lei，et al.Reviews of the model of projection matrix of CT reconstruction algorithm［J］.CT Theory and Applications，2014，23（2）：317-328.（in Chinese）

［13］查國震.基于正六邊形像素的扇束等距濾波反投影及平行束插值代數(shù)重建算法研究［D］.北京：首都師范大學(xué)，2009.

［14］徐鍵，邱鈞.一種基于計(jì)算點(diǎn)的圖像重建離散化模型［C］.應(yīng)用數(shù)學(xué)暑期研討會論文集，電子工業(yè)出版社，2011.

［15］孫守思，邱鈞，桂志國，等.重建點(diǎn)模型的EM迭代成像［J］.中北大學(xué)學(xué)報(bào)（自然科學(xué)版），2014，35（2）：209-217.Sun Shousi，Qiu Jun，Gui Zhiguo，et al.The EM iterative imaging on the discrete reconstruction points［J］.Journal of North University of China（Natural Science Edition）.2014，35（2）：209-217.（in Chinese）

［16］王超，王化祥.電阻抗斷層圖像重建算法研究——預(yù)迭代算法提出［J］.信號處理，2002，18（6）：547-550.Wang Chao，Wang Huaxiang.Research on the reconstruction algorithm of electrical impedance tomography——the pre-iteartion reconstruction algorithm［J］.Singal Processing，2002，18（6）：547-550.（in Chinese）

［17］張兆田.不完全數(shù)據(jù)斷層成像理論與應(yīng)用研究［D］.北京：中國科學(xué)院，2005.

［18］柯麗，曹馮秋，杜強(qiáng).MIT中反投影矩陣的計(jì)算與數(shù)據(jù)處理方法［J］.儀器儀表學(xué)報(bào)，2014，35（10）：2256-2262.Ke Li，Cao Fengqiu，Du Qiang.Back-projection matrix calculation and data processing methods used in magnetic induction tomography［J］.Chinese Journal of Scientific Instrument，2014，35（10）：2256-2262.（in Chinese"）