趙芳薇，邱 鈞，劉 暢

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

0 引 言

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

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

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

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

1 重建點離散化模型

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

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

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

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

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

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

2 反投影模型

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

設待重建區(qū)域被離散化為N個理想小區(qū)域，有M個有效投影數(shù)據(jù).反投影模型刻畫了依據(jù)不同掃描射線對某重建點的貢獻因子大小，進行圖像重建的過程.記tji是第i根射線對第j個重建點的貢獻因子，則第j個重建點像素值可記為式中：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ù)對重建點的貢獻因子已知的情形下，圖像重建問題歸結(jié)為由全體投影向量b和反投影系數(shù)矩陣T直接計算圖像X的問題.

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

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

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

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

式中：si是第i根射線被重建區(qū)域截得的長度；dij是第j個重建點到第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ù)的估計.算法的分量迭代格式為

式中：M為投影射線總數(shù)；N為重建點的總數(shù)；bi為投影數(shù)據(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），計算投影值與投影數(shù)據(jù)bi，按照式（5）進行迭代.k=0，1，2，….

5 算法驗證的數(shù)據(jù)實驗

本節(jié)采用兩組模擬數(shù)據(jù)驗證算法，數(shù)據(jù)實驗中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

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

5.1 模擬數(shù)據(jù)實驗1

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

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

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

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

5.2 模擬數(shù)據(jù)實驗2

模擬數(shù)據(jù)實驗2選取鸚鵡模型如圖8和圖9所示.迭代5，10，20，30次的重建結(jié)果及第200行對應的像素曲線如圖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 實驗2反投影模型下迭代算法迭代不同輪次的重建結(jié)果 Fig.10 Reconstruction results of the new iterative algorithm under different numbers of iterations for experiment two

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

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

6 結(jié) 論

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

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

［1］莊天戈.CT原理與算法［M］.上海：上海交通大學出版社，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］.電子與信息學報，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ù)重建算法的快速實現(xiàn)［J］.核電子學與探測技術(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理論與應用研究，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］.北京：首都師范大學，2009.

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

［15］孫守思，邱鈞，桂志國，等.重建點模型的EM迭代成像［J］.中北大學學報（自然科學版），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］王超，王化祥.電阻抗斷層圖像重建算法研究——預迭代算法提出［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ù)斷層成像理論與應用研究［D］.北京：中國科學院，2005.

［18］柯麗，曹馮秋，杜強.MIT中反投影矩陣的計算與數(shù)據(jù)處理方法［J］.儀器儀表學報，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"）