改進投影數據訪問順序的CT圖像重建研究

張蕾+徐伯慶

摘 要：代數重建算法是CT圖像重建中較為常用的圖像重建算法，重建圖像的質量和速度是評價重建算法優劣的兩個重要標準。傳統的代數重建算法有重建質量低、速度慢的缺點。針對這些缺點，提出了一種改進投影數據訪問順序的代數重建算法。該算法通過對投影訪問順序的選擇，減小了投影數據間的相關性，同時保證了投影角度的均勻分布，避免了密集訪問情況的出現。實驗結果證明，采用改進的算法進行CT圖像重建不僅能減少圖像偽影，而且能很好地改善重建圖像質量，明顯加快重建速度。

關鍵詞：CT圖像重建；代數重建；投影數據；訪問順序；最小相關性

DOIDOI：10.11907/rjdk.171337

中圖分類號：TP317.4

文獻標識碼：A 文章編號文章編號：1672-7800（2017）008-0189-04

0 引言

計算機斷層成像（Computed Tomography， CT）技術從外部測量物體得到數據重建內部信息，普遍應用于醫學、無損檢測等領域[1]。代數重建法（ART）和濾波反投影法（FBP）是CT圖像重建的兩類經典算法[2]。FBP算法重建的圖像質量好且速度快，但其條件是有充足的投影數據，而現實操作中卻很難做到。ART算法不需要這一條件，但其重建速度會減慢[3]。將現實問題轉化為運用代數方式解線性方程組的數學問題是ART算法的核心思想[4]。初值選取、投影數據訪問順序、松弛參數選擇等因素對重建過程都有不同程度的影響[5-7]。其中，投影數據訪問順序對重建速度和質量影響較大。為此，學者進行了大量研究，如固定角度訪問順序、Van Dijke[8]提出的隨機訪問順序等。通過深入分析與研究，本文提出一種基于改進投影數據訪問順序的ART算法，以求改善CT圖像重建的精度和速率。

1 ART算法原理

ART算法首先將被重建的物體離散化為一個個像素或體素，這些圖像數據就相當于一個未知的數學矩陣，通過不同投影角度獲取的投影數據聯合這一未知的數學矩陣建立線性方程組，用迭代的方式求解方程組便可得到原始圖像矩陣[9]。本文以二維圖像為例，將待重建圖像離散化為N=m×n個圖像像素，待重建區域劃分為一個個網格，并用xi（1≤i≤N）表示第i個網格內像素的平均值，如圖1所示。

在此，采用平行束投影方式進行數據采集。平行射線的總數為M，圖像像素的總數為N，待重建區域由N個邊長為1的正方形網格組成。按照CT圖像重建原理以及構造的數學模型，投影數據與像素值之間的關系可用下式表示：

W11x1+W12x2+…+W1NxN=p1W21x1+W22x2+…+W2NxN=p2…WM1x1+WM2x2+…+WMNxN=pM （1）

式（1）中，Wij為權重因子，即第i條射線對第j個像素的奉獻值，pi（1≤i≤M）為待重建圖像在第i條射線上的投影值，即待重建圖像沿此射線方向的線積分。由于M、N的值通常都很大，且每條射線與像素的交點個數較少，權重因子矩陣成為大型稀疏矩陣，所以只能采用迭代方式求解。

ART 算法通過對未知圖像矩陣賦初值，再將射線掃描圖像得到的實際投影數據與估計的投影數據進行比較，使用迭代公式對初值進行校正與更改，便可完成第一次迭代。若一次迭代后未滿足收斂條件，則將更新后的圖像矩陣作為下一次的初值，繼續以上步驟直至達到收斂條件[10]。ART算法的迭代公式如下：

x（n + 1）j = x（n）j + λWij ∑Nj = 1W2ij pi -∑Nj = 1Wij x（n）j（2）

式（2）中，n為迭代次數，λ為松弛參數（0<λ<2），且i=1，2，…，M，j=1，2，…，N。

2 投影數據訪問順序改進

2.1 最小相關性原理

將方程組（1）中的M個方程抽象化為N維空間內的M個超平面，若方程組的解只有一個，M個超平面的交點便是其解。若M>N，方程組便不會有唯一解，超平面的交點附近會產生震蕩；若M

初始值設為x0，將x0投影至第一個超平面，再將得到的值垂直投影至第二個超平面，反復操作直至達到兩個超平面的交點。從圖2得出結論：到達交點需要的迭代次數隨著夾角的增大而減少。若超平面的夾角為直角，迭代2次便可完成。同時，超平面的夾角較大時，可獲得較高精度的最終解。因此，為了達到收斂速度快、重建質量好的目的，需要對投影數據的訪問順序進行選擇，增大相鄰超平面間的夾角，減小數據相關性，這就是最小相關性原理。選擇高效的投影數據訪問順序可在迭代初期較快重建高頻成分[12-14]。這一高效的投影數據訪問順序除了需要滿足最小相關性原理外，還需要遵循以下原則：①在視角范圍內，投影角度應盡可能分布均勻；②不允許投影角度在某些角度被密集訪問。

2.2 改進的投影數據訪問順序

改進的投影數據訪問順序充分運用最小相關性原理，而且投影角度均勻分布，未出現投影角度密集訪問情況。假設投影角度范圍為[0，θ]，采用間隔為α的均勻采樣，投影數據訪問順序改進步驟如下：①設i=0為首個投影角度；②投影角度按以下循環體選擇：j=i，j<θ，j=j+90；③i=i+α，重復步驟②，直至遍歷所有投影角度。通過以上方式選擇的投影數據訪問順序，使得相鄰射線間的夾角盡量滿足90°，從而減小射線穿過像素的信息相關性。如θ=360，α=3時，按照以上步驟可得投影數據的訪問順序為：0，90，180，270，360，3，93，183，273，…。

3 實驗結果與分析

3.1 圖像質量評價方式

圖像質量的好壞可用主觀和客觀兩種方式評價。主觀評價就是通過人眼進行直觀判斷，這種方式雖然直觀可見，但容易受到人的主觀影響，缺乏科學依據。客觀評價就是通過計算機模擬計算出圖像的各項衡量指標，結合原始圖像和重建圖像的匹配程度來評價。為了直觀對比，本文綜合這兩種方式來評估圖像質量，并采用以下3項指標作為客觀評價標準：（1）密度曲線對比圖，即根據圖像的某一列或某一行的像素值繪制密度曲線圖，對比分析重建圖像與原始圖像像素的偏離大小，直觀判斷重建算法的優劣。（2）峰值信噪比PSNR，設xi 表示重建后圖像的像素值，其均方誤差MSE可表示為：

MSE = 1M×N∑0≤i < N（xi -xi ）2（3）則峰值信噪比PSNR可表示為： PSNR=10lg255×255MSE（4）由式（4）可知，PSNR的值越大，重建圖像中的噪聲越小，重建的質量越好。（3）歸一化平均絕對距離r，計算公式如下： r = ∑0≤i < Nxi -xi ∑0≤i < Nxi （5）由式（5）可知，r值越小，圖像失真越小，重建效果越好。

3.2 實驗結果與分析

為了驗證改進算法的正確性與可行性，選用集算法研究、建模仿真、可視化、數據分析等功能為一體的MATALB仿真軟件進行仿真實驗。實驗采用Sheep-Logan頭部模型切片來模擬分析不同投影數據訪問順序對CT圖像重建的影響。該頭部模型是CT圖像重建領域內的經典仿真模型，像素取值范圍為[0，255]，實驗采用的模型大小為128×128。

實驗采用4種算法分別在無噪聲條件下進行圖像重建，投影數據訪問順序各不相同，分別為：順序訪問（SAR）、隨機訪問（RAS）、固定角度訪問（FAS）、改進的投影數據訪問（PAS），其中FAS采用的固定角度為67°。實驗中的投影射線均為平行束投影，投影角度取值范圍為[0，180°]，均勻采樣間隔為1°，即投影射線每間隔1°掃描一次，185條投影射線平均分布在每個投影角度上，因此獲得的投影矩陣大小為185×180。圖像的初值為0，松弛參數為1，迭代次數為6次，圖3為4種訪問方式重建的圖像，圖4為在x=64即垂直方向上各個重建圖像的像素值圖。

由圖3可以看出4種不同訪問順序的ART算法均能成功重建圖像，但都有一定條狀偽影。從圖像可以看出，采用改進的投影數據訪問順序重建出的圖像較之其它3種方式重建的圖像愈發清晰，偽影較少。由圖4可以發現，基于這4種不同訪問順序得到的重建圖像像素值較之原始圖像像素值差別均不大，但是采用改進的投影數據訪問順序的密度曲線與原始圖像更為接近，誤差最小。

為了更加客觀地分析對比圖像重建的質量和速度，不斷增加迭代次數，觀察隨著迭代次數的增加，采用不同訪問順序重建圖像的值和值的變化。圖5所示為隨著迭代次數的增加兩種評價標準的變化曲線，其中（a）為PSNR值變化曲線、（b）為r值變化曲線。由圖5可知，改進的投影數據訪問順序重建的圖像PSNR值最大，r值最小，即圖像重建質量最好。隨著迭代次數增加，收斂速度也越快，最終迭代6次后，改進的CT圖像重建算法趨于穩定，且重建圖像質量始終優于其它訪問順序的算法。

使用MATLAB仿真軟件記錄程序運行時間t，進一步評價各種算法效率。表1所示為迭代6次時，不同投影數據訪問順序評價標準值。由表1可以看出，在迭代6次時，本文算法得到的重建圖像誤差最小，質量最好，同時運行時間最短。以上各項實驗數據和仿真圖像，充分證實了改進的算法不但提高了圖像的重建質量，還實現了重建算法的快速收斂。

4 結語

CT圖像重建影響因素眾多，本文對投影數據的訪問順序作了改進。對實驗結果分析發現，投影數據的訪問順序直接決定著圖像質量的優劣和重建速度的快慢。本文提出的投影數據訪問順序不僅在視角范圍內均勻分布且不密集出現，而且滿足了最小相關原理，使得連續投影射線間的投影數據相關性最小。改進算法不僅大大提高了重建質量，而且明顯加快了重建速度。本文的仿真實驗均在無噪聲條件下進行，但實際環境中CT掃描難免受到噪聲污染，今后將進一步研究在有噪聲情況下CT圖像重建算法的改進，使其更好地應用到實際場景中。

參考文獻：

[1] HERMAN G T.Image reconstruction from projections[M].London：ACADADEMIC P-RESS， 1980.

[2] 徐國良，陳沖，李明.圖像重構的數值方法[M].北京：科學出版社，2015.

[3] 莊天戈.CT原理與算法[M].上海：上海交通大學出版社，1992.

[4] NATTERER F.The mathematics of computerized tomography[M].Society for Industrial Mathematics， 2001.

[5] 楊文良，魏東波.一種改進投影系數計算的快速ART算法[J].CT理論與應用研究，2012， 21（2）： 187-195.

[6] GUAN H，GORDON R.A projection access order for speedy convergence of ART：a multilevel scheme for computed tomography[J]. Physics in Medicine and Biology， 1994， 39（11）： 2005-2022.

[7] HERMAN G， MEYER L. Algebraic reconstruction can be made computationally efficient[J]. IEEE Trans， Med Img，1993， 12（3）： 600-609.endprint

[8] VAN DIJKE.Iterative methods in image reconstruction[M]. The Netherlands， 1992.

[9] KESIDIS A L，PAPAMARKOS N.Exact image reconstruction from a limited number of projections[J]. Journal of Visual Communication and Image Reconstruction， 2008， 19（5）： 285-298.

[10] 李志鵬，叢鵬，鄔海峰.代數迭代算法進行CT圖像重建的研究[J].核電子學與探測技術，2005，25（2）： 184-186.

[11] MAN B D. Statistics methods for image reconstruction part 3： x-ray computed tomography[C]. IEEE Nuclear Science Symposium and Medical Imaging Conference， Honollu， 2007.

[12] 王宏鈞，路宏年，傅建.代數重建技術中投影序列選擇次序的研究[J].光學技術，2006，32（3）： 1005-1022.

[13] THIBAULT JB， SAUER KD， BOUMAN CA， et al. A three-dimensional statistical approach to improve image quality for multislice helical CT[J]. Medical Physics， 2007， 34（11）： 4526-4544.

[14] MUELLER K， YAGEL R， CORNHILL J F. The weighted distance scheme： a globally optimizing projection ordering method for the algebraic reconstruction technique （ART） [J]. IEEE Transaction on Medical Image， 1997， 16（2）： 223-230.endprint